diff --git a/bn_mp_next_small_prime.c b/bn_mp_next_small_prime.c new file mode 100644 index 000000000..4ba19114c --- /dev/null +++ b/bn_mp_next_small_prime.c @@ -0,0 +1,35 @@ +#include "tommath_private.h" +#ifdef BN_MP_NEXT_SMALL_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + +/* + * Mimics behaviour of Pari/GP's nextprime(n) + * If n is prime set *result to n else set *result to first prime > n + * and 0 in case of error + */ +mp_err mp_next_small_prime(mp_sieve_prime n, mp_sieve_prime *result, mp_sieve *sieve) +{ + mp_sieve_prime ret = 0; + int e = MP_OKAY; + + if (n < 2) { + *result = 2; + return e; + } + + for (; ret == 0 && n != 0; n++) { + if (n > MP_SIEVE_BIGGEST_PRIME) { + return MP_SIEVE_MAX_REACHED; + } + /* just call mp_is_small_prime(), it does all of the heavy work */ + if ((e = mp_is_small_prime(n, &ret, sieve)) != MP_OKAY) { + *result = 0; + return e; + } + } + *result = n - 1; + return e; +} +#endif + diff --git a/bn_mp_prec_small_prime.c b/bn_mp_prec_small_prime.c new file mode 100644 index 000000000..2797838c6 --- /dev/null +++ b/bn_mp_prec_small_prime.c @@ -0,0 +1,37 @@ +#include "tommath_private.h" +#ifdef BN_MP_PREC_SMALL_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + +/* + * Mimics behaviour of Pari/GP's precprime(n) + * If n is prime set *result to n else set *result to first prime < n + * and 0 in case of error + */ +mp_err mp_prec_small_prime(mp_sieve_prime n, mp_sieve_prime *result, mp_sieve *sieve) +{ + mp_sieve_prime ret = 0; + int e = MP_OKAY; + + if (n == 2) { + *result = 2; + return e; + } + + if (n < 2) { + *result = 0; + return e; + } + + for (; ret == 0; n--) { + if ((e = mp_is_small_prime(n, &ret, sieve)) != MP_OKAY) { + *result = 0; + return e; + } + } + *result = n + 1; + + return e; +} +#endif + diff --git a/bn_mp_sieve.c b/bn_mp_sieve.c new file mode 100644 index 000000000..ffc710d8f --- /dev/null +++ b/bn_mp_sieve.c @@ -0,0 +1,393 @@ +#include "tommath_private.h" +#ifdef BN_MP_SIEVE_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + +/* Segmented version of an Eratosthenes sieve */ +#define mp_sieve_prime_NUM_BITS (sizeof(mp_sieve_prime)*CHAR_BIT) + +static void s_mp_sieve_setall(mp_single_sieve *bst); +static void s_mp_sieve_clear(mp_single_sieve *bst, mp_sieve_prime n); +static mp_sieve_prime s_mp_sieve_get(mp_single_sieve *bst, mp_sieve_prime n); +static mp_sieve_prime s_mp_sieve_nextset(mp_single_sieve *bst, mp_sieve_prime n); + +static mp_sieve_prime s_isqrt(mp_sieve_prime n); + +static void s_mp_eratosthenes(mp_single_sieve *bst); +static int s_mp_eratosthenes_init(mp_sieve_prime n, mp_single_sieve *bst); + +static void s_mp_eratosthenes_segment(mp_sieve_prime a, mp_sieve_prime b, mp_single_sieve *base, + mp_single_sieve *segment); +static int s_mp_eratosthenes_segment_init(mp_sieve_prime a, mp_sieve_prime b, mp_single_sieve *base, + mp_single_sieve *segment); +static void s_mp_eratosthenes_segment_clear(mp_single_sieve *segment, mp_sieve_prime *single_segment_a); + +static int s_init_base_sieve(mp_single_sieve *base_sieve); +static int s_init_single_segment_with_start(mp_sieve_prime a, mp_single_sieve *base_sieve, + mp_single_sieve *single_segment, mp_sieve_prime *single_segment_a); + +static mp_sieve_prime s_nextprime(mp_sieve_prime p, mp_single_sieve *bst); + + +/* Manual memset to avoid the inclusion of string.h */ +static void s_mp_sieve_setall(mp_single_sieve *bst) +{ + mp_sieve_prime i, bs_size; + bs_size = bst->alloc / sizeof(mp_sieve_prime); + for (i = 0; i < bs_size; i++) { + (bst)->content[i] = MP_SIEVE_PRIME_MAX; + } +} + +static void s_mp_sieve_clear(mp_single_sieve *bst, mp_sieve_prime n) +{ + ((*((bst)->content+(n/mp_sieve_prime_NUM_BITS)) + &= ~(1lu<<(n % mp_sieve_prime_NUM_BITS)))); +} + +static mp_sieve_prime s_mp_sieve_get(mp_single_sieve *bst, mp_sieve_prime n) +{ + return (((*((bst)->content+(n/mp_sieve_prime_NUM_BITS)) + & (1lu<<(n % mp_sieve_prime_NUM_BITS))) != 0)); +} + +static mp_sieve_prime s_mp_sieve_nextset(mp_single_sieve *bst, mp_sieve_prime n) +{ + while ((n < ((bst)->size)) && (!s_mp_sieve_get(bst, n))) { + n++; + } + return n; +} + + +/* + * Integer square root, hardware style + * Wikipedia calls it "Digit-by-digit calculation" + * https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation + * This is the base 2 method described at + * https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2) + */ +static mp_sieve_prime s_isqrt(mp_sieve_prime n) +{ + mp_sieve_prime s, rem, root; + + if (n < 1) { + return 0; + } + /* highest power of four <= n */ + s = ((mp_sieve_prime) 1) << (mp_sieve_prime_NUM_BITS - 2); + + rem = n; + root = 0; + while (s > n) { + s >>= 2; + } + while (s != 0) { + if (rem >= (s | root)) { + rem -= (s | root); + root >>= 1; + root |= s; + } else { + root >>= 1; + } + s >>= 2; + } + return root; +} + +/* + * Set range_a_b to sqrt(MP_SIEVE_PRIME_MAX) + * TODO: Make it const or put it in bncore.c because it is said to be faster + * if the size of range_a_b fits into the L2-cache. + * Not much difference on the author's machine for 32 bit but quite + * a large one for 64 bit and large limits. YMMV, as always. + * Please be aware that range_a_b is in bits, not bytes and memory + * allocation rounds up and adds CHAR_BIT*sizeof(mp_sieve_prime) bits. + */ +#ifndef MP_SIEVE_RANGE_A_B +#if ( (defined MP_64BIT) && (defined MP_SIEVE_USE_LARGE_SIEVE) ) +#define MP_SIEVE_RANGE_A_B 0x400000uL +#else +#define MP_SIEVE_RANGE_A_B ((mp_sieve_prime) MP_SIEVE_PRIME_MAX_SQRT) +#endif +#endif +#define MP_SIEVE_BASE_SIEVE_SIZE ((mp_sieve_prime)MP_SIEVE_PRIME_MAX_SQRT) + +/* TODO: Some redundant code below, needs a clean-up */ + + +/* + * Initiate a sieve that stores the odd numbers only: + * allocate memory, set actual size and allocated size and fill it + */ +static int s_mp_eratosthenes_init(mp_sieve_prime n, mp_single_sieve *bst) +{ + n = (n - 1) / 2; + bst->content = + MP_MALLOC(((size_t)n + sizeof(mp_sieve_prime)) / sizeof(mp_sieve_prime) + sizeof(mp_sieve_prime)); + if (bst->content == NULL) { + return MP_MEM; + } + + bst->alloc = + (n + sizeof(mp_sieve_prime)) / sizeof(mp_sieve_prime) + sizeof(mp_sieve_prime); + bst->size = 2 * n + 1; + + s_mp_eratosthenes(bst); + return MP_OKAY; +} + +/* + * Simple Eratosthenes' sieve, starting at zero + * Keeping odd primes only as the single optimization + */ +static void s_mp_eratosthenes(mp_single_sieve *bst) +{ + mp_sieve_prime n, k, r, j; + + n = bst->size; + r = s_isqrt(n); + s_mp_sieve_setall(bst); + for (k = 1; k < ((r - 1) / 2); k += 1) { + if (s_mp_sieve_get(bst, k)) { + for (j = k * 2 * (k + 1); j < (n - 1) / 2; j += 2 * k + 1) { + s_mp_sieve_clear(bst, j); + } + } + } +} + +/* + * For a sieve that has only the odd numbers stored. + * It returns next prime even if p is prime in contrast to + * e.g.: Pari/GP's nextprime() function + * Used internally for the segment to get the primes from + * the base sieve. + */ +static mp_sieve_prime s_nextprime(mp_sieve_prime p, mp_single_sieve *bst) +{ + mp_sieve_prime ret; + if (p <= 1) + return 2; + ret = s_mp_sieve_nextset(bst, ((p - 1) / 2) + 1); + return 2 * ret + 1; +} + +/* + * Init sieve "segment" of the range [a,b]: + * allocate memory, set actual size and allocated size and fill it from sieve "base" + * TODO: merge with s_mp_eratosthenes_init() + */ +static int s_mp_eratosthenes_segment_init(mp_sieve_prime a, mp_sieve_prime b, + mp_single_sieve *segment, mp_single_sieve *base) +{ + mp_sieve_prime n; + + n = b - a; + + segment->content = + MP_MALLOC(((size_t)n + sizeof(mp_sieve_prime)) / sizeof(mp_sieve_prime) + + sizeof(mp_sieve_prime)); + if (segment->content == NULL) { + return MP_MEM; + } + segment->alloc = + (n + sizeof(mp_sieve_prime)) / sizeof(mp_sieve_prime) + + sizeof(mp_sieve_prime); + segment->size = n; + + s_mp_eratosthenes_segment(a, b, base, segment); + return MP_OKAY; +} + +/* + * Clear sieve "segment" + * free memory and reset "single_segment_a" + */ +static void s_mp_eratosthenes_segment_clear(mp_single_sieve *segment, + mp_sieve_prime *single_segment_a) +{ + if (segment->content != NULL) { + MP_FREE(segment->content, ((size_t)segment->size + sizeof(mp_sieve_prime)) / + sizeof(mp_sieve_prime) + sizeof(mp_sieve_prime)); + } + *single_segment_a = 0; + +} + +/* + * TODO: Not memory optimized, it stores the even numbers, too. + * Should be ok if the segments are small but needs to + * be done at some time in the near future. + * + * Fill sieve "segment" of the range [a,b] from the basic sieve "base" + */ +static void s_mp_eratosthenes_segment(mp_sieve_prime a, mp_sieve_prime b, + mp_single_sieve *base, mp_single_sieve *segment) +{ + mp_sieve_prime r, j, i; + mp_sieve_prime p; + r = s_isqrt(b); + s_mp_sieve_setall(segment); + p = 0; + for (i = 0; p <= r; i++) { + p = s_nextprime(p, base); + j = p * p; + if (j < a) { + j = ((a + p - 1) / p) * p; + } + for (; j <= b; j += p) { + /* j+=p can overflow */ + if (j >= a) { + s_mp_sieve_clear(segment, j - a); + } else { + break; + } + } + } +} + +/* Build a basic sieve with the largest reasonable size */ +static int s_init_base_sieve(mp_single_sieve *base_sieve) +{ + mp_sieve_prime n; + /* + * That is a risky idea. If range_a_b is smaller than sqrt(n_max) + * the whole thing stops working! + */ + /* n = range_a_b; */ + n = MP_SIEVE_BASE_SIEVE_SIZE; + return s_mp_eratosthenes_init(n, base_sieve); +} + +/* + * Build a segment sieve with the largest reasonable size. "a" is the start of + * the sieve Size is MIN(range_a_b,MP_SIEVE_PRIME_MAX-a) + */ +static int s_init_single_segment_with_start(mp_sieve_prime a, + mp_single_sieve *base_sieve, + mp_single_sieve *single_segment, + mp_sieve_prime *single_segment_a) +{ + mp_sieve_prime b; + int e = MP_OKAY; + + /* last segment might not fit, depending on size of range_a_b */ + if (a > (MP_SIEVE_PRIME_MAX - MP_SIEVE_RANGE_A_B)) { + b = MP_SIEVE_PRIME_MAX - 1; + } else { + b = a + MP_SIEVE_RANGE_A_B; + } + if ((e = + s_mp_eratosthenes_segment_init(a, b, single_segment, + base_sieve)) != MP_OKAY) { + return e; + } + *single_segment_a = a; + return e; +} + +/* + * Sets "result" to one if n is prime or zero respectively. + * Also sets "result" to zero in case of error. + * Worst case runtime is: building a base sieve and a segment and + * search the segment + */ +mp_err mp_is_small_prime(mp_sieve_prime n, mp_sieve_prime *result, mp_sieve *sieve) +{ + int e = MP_OKAY; + mp_sieve_prime a = 0, b = 0; + + if (n < 2) { + *result = 0; + return e; + } + + if ((n & 1) == 0) { + *result = 0; + return e; + } + + /* neither of 2^16-x, 2^32-x, or 2^64-x are prime for 0<=x<=4 */ + if (n >= (MP_SIEVE_PRIME_MAX - 3)) { + *result = 0; + return e; + } + + if (sieve->base.content == NULL) { + if ((e = s_init_base_sieve(&(sieve->base))) != MP_OKAY) { + *result = 0; + return e; + } + } + + /* No need to generate a segment if n is in the base sieve */ + if (n < MP_SIEVE_BASE_SIEVE_SIZE) { + /* might have been a small sieve, so check size of sieve first */ + if (n < sieve->base.size) { + *result = s_mp_sieve_get(&(sieve->base), (n - 1) / 2); + return e; + } + } + + /* no further shortcuts to apply, build and search a segment */ + + /* + * TODO: if base_sieve is < sqrt(MP_SIEVE_PRIME_MAX) it is not possible to get + * all primes <= MP_SIEVE_PRIME_MAX so add check if n > size(base_sieve)^2 and either + * a. make size(base_sieve) = sqrt(MP_SIEVE_PRIME_MAX) + * b. make size(base_sieve) = 2 * size(base_sieve) + * with 2 * size(base_sieve) <= sqrt(MP_SIEVE_PRIME_MAX) + * c. give up and return MP_VAL + */ + /* we have a segment and may be able to use it */ + if (sieve->segment.content != NULL) { + a = sieve->single_segment_a; + /* last segment may not fit into range_a_b */ + if (a > (MP_SIEVE_PRIME_MAX - MP_SIEVE_RANGE_A_B)) { + b = MP_SIEVE_PRIME_MAX - 1; + } else { + b = a + MP_SIEVE_RANGE_A_B; + } + /* check if n is inside the bounds of the segment */ + if (n >= a && n <= b) { + *result = s_mp_sieve_get(&(sieve->segment), n - a); + return e; + } + /* get a clean slate */ + else { + s_mp_eratosthenes_segment_clear(&(sieve->segment), &(sieve->single_segment_a)); + } + } + + /* + * A bit of heuristics ( "heuristics" is a more pretentious word for the + * commonly known expression "wild guess") + * Based on the vague idea of the assumption that most sieves get used for + * sequential series of primes or a single test here and there, but not + * for massive amounts of random requests. + */ + if (n > a) { + if (n > (MP_SIEVE_PRIME_MAX - MP_SIEVE_RANGE_A_B)) { + a = MP_SIEVE_PRIME_MAX - MP_SIEVE_RANGE_A_B; + } else { + a = n; + } + } else { + if (n < MP_SIEVE_RANGE_A_B) { + a = MP_SIEVE_BASE_SIEVE_SIZE - 1; + } else { + /* TODO: there should be no need for an overlap, check */ + a = n - (MP_SIEVE_RANGE_A_B - 2); + } + } + if ((e = s_init_single_segment_with_start(a, + &(sieve->base), &(sieve->segment), &(sieve->single_segment_a))) != MP_OKAY) { + *result = 0; + return e; + } + /* finally, check for primality */ + *result = s_mp_sieve_get(&(sieve->segment), n - a); + return e; +} +#endif diff --git a/bn_mp_sieve_clear.c b/bn_mp_sieve_clear.c new file mode 100644 index 000000000..b311f95f9 --- /dev/null +++ b/bn_mp_sieve_clear.c @@ -0,0 +1,21 @@ +#include "tommath_private.h" +#ifdef BN_MP_SIEVE_CLEAR_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + +/* Free memory of one sieve */ +static void s_clear_one(mp_single_sieve *sieve) +{ + if (sieve->content != NULL) { + MP_FREE(sieve->content, ((size_t)segment->size + sizeof(mp_sieve_prime)) / + sizeof(mp_sieve_prime) + sizeof(mp_sieve_prime)); + } +} + +void mp_sieve_clear(mp_sieve *sieve) +{ + s_clear_one(&(sieve->base)); + s_clear_one(&(sieve->segment)); +} +#endif + diff --git a/bn_mp_sieve_init.c b/bn_mp_sieve_init.c new file mode 100644 index 000000000..55fc73b66 --- /dev/null +++ b/bn_mp_sieve_init.c @@ -0,0 +1,14 @@ +#include "tommath_private.h" +#ifdef BN_MP_SIEVE_INIT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + +/* Set the default values for the sieve */ +void mp_sieve_init(mp_sieve *sieve) +{ + sieve->base.content = NULL; + sieve->segment.content = NULL; + sieve->single_segment_a = 0; +} +#endif + diff --git a/callgraph.txt b/callgraph.txt index 4eb61ea68..b5b5495e5 100644 --- a/callgraph.txt +++ b/callgraph.txt @@ -3879,6 +3879,9 @@ BN_MP_NEG_C | +--->BN_MP_GROW_C +BN_MP_NEXT_SMALL_PRIME_C + + BN_MP_N_ROOT_C +--->BN_MP_N_ROOT_EX_C | +--->BN_MP_INIT_C @@ -4501,6 +4504,9 @@ BN_MP_OR_C +--->BN_MP_CLEAR_C +BN_MP_PREC_SMALL_PRIME_C + + BN_MP_PRIME_FERMAT_C +--->BN_MP_CMP_D_C +--->BN_MP_INIT_C @@ -12593,6 +12599,15 @@ BN_MP_SET_LONG_LONG_C BN_MP_SHRINK_C +BN_MP_SIEVE_C + + +BN_MP_SIEVE_CLEAR_C + + +BN_MP_SIEVE_INIT_C + + BN_MP_SIGNED_BIN_SIZE_C +--->BN_MP_UNSIGNED_BIN_SIZE_C | +--->BN_MP_COUNT_BITS_C diff --git a/demo/test.c b/demo/test.c index 69c8fbc80..54adb8f9f 100644 --- a/demo/test.c +++ b/demo/test.c @@ -1560,6 +1560,161 @@ static int test_mp_decr(void) mp_clear_multi(&a, &b, NULL); return EXIT_FAILURE; } +static int test_mp_is_small_prime(void) +{ + mp_sieve sieve; + int e; + int i, test_size; + + mp_sieve_prime to_test[] = { + 52, 137, 153, 179, 6, 153, 53, 132, 150, 65, + 27414, 36339, 36155, 11067, 52060, 5741, + 29755, 2698, 52572, 13053, 9375, 47241 +#ifndef MP_8BIT + ,39626, 207423, 128857, 37419, 141696, 189465, + 41503, 127370, 91673, 8473, 479142414, 465566339, + 961126169, 1057886067, 1222702060, 1017450741, + 1019879755, 72282698, 2048787577, 2058368053 +#endif + }; + mp_sieve_prime tested[] = { + 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 +#ifndef MP_8BIT + ,0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0 +#endif + }; + mp_sieve_prime result; + + mp_sieve_init(&sieve); + + test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime)); + + for (i = 0; i < test_size; i++) { + if ((e = mp_is_small_prime(to_test[i], &result, &sieve)) != MP_OKAY) { + fprintf(stderr,"mp_is_small_prime failed: \"%s\"\n",mp_error_to_string(e)); + goto LTM_ERR; + } + if (result != tested[i]) { + fprintf(stderr,"mp_is_small_prime failed for %lu. Said %lu but is %lu \n", + (unsigned long)to_test[i], (unsigned long)result, (unsigned long)tested[i]); + goto LTM_ERR; + } + } + mp_sieve_clear(&sieve); + return EXIT_SUCCESS; +LTM_ERR: + mp_sieve_clear(&sieve); + return EXIT_FAILURE; +} + +static int test_mp_next_small_prime(void) +{ + mp_sieve sieve; + mp_sieve_prime ret; + int e; + int i, test_size; + + mp_sieve_prime to_test[] = { + 52, 137, 153, 179, 6, 153, 53, 132, 150, 65, + 27414, 36339, 36155, 11067, 52060, 5741, + 29755, 2698, 52572, 13053, 9375, 47241 +#ifndef MP_8BIT + ,39626, 207423, 128857, 37419, 141696, 189465, + 41503, 127370, 91673, 8473, 479142414, 465566339, + 961126169, 1057886067, 1222702060, 1017450741, + 1019879755, 72282698, 2048787577, 2058368053 +#endif + }; + mp_sieve_prime tested[] = { + 53, 137, 157, 179, 7, 157, 53, 137, 151, 67, + 27427, 36341, 36161, 11069, 52067, 5741, + 29759, 2699, 52579, 13063, 9377, 47251 +#ifndef MP_8BIT + ,39631, 207433, 128857, 37423, 141697, 189467, + 41507, 127373, 91673, 8501, 479142427, 465566393, + 961126169, 1057886083, 1222702081, 1017450823, + 1019879761, 72282701, 2048787577, 2058368113 +#endif + }; + + mp_sieve_init(&sieve); + + test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime)); + + for (i = 0; i < test_size; i++) { + if ((e = mp_next_small_prime(to_test[i], &ret, &sieve)) != MP_OKAY) { + fprintf(stderr,"mp_next_small_prime failed with \"%s\" at index %d\n", + mp_error_to_string(e), i); + goto LTM_ERR; + } + if (ret != tested[i]) { + fprintf(stderr,"mp_next_small_prime failed for %lu. Said %lu but is %lu \n", + (unsigned long)to_test[i], (unsigned long)ret, (unsigned long)tested[i]); + goto LTM_ERR; + } + } + mp_sieve_clear(&sieve); + return EXIT_SUCCESS; +LTM_ERR: + mp_sieve_clear(&sieve); + return EXIT_FAILURE; +} + +static int test_mp_prec_small_prime(void) +{ + mp_sieve sieve; + mp_sieve_prime ret; + int e; + int i, test_size; + + mp_sieve_prime to_test[] = { + 52, 137, 153, 179, 6, 153, 53, 132, 150, 65, + 27414, 36339, 36155, 11067, 52060, 5741, + 29755, 2698, 52572, 13053, 9375, 47241 +#ifndef MP_8BIT + ,39626, 207423, 128857, 37419, 141696, 189465, + 41503, 127370, 91673, 8473, 479142414, 465566339, + 961126169, 1057886067, 1222702060, 1017450741, + 1019879755, 72282698, 2048787577, 2058368053 +#endif + }; + mp_sieve_prime tested[] = { + 47, 137, 151, 179, 5, 151, 53, 131, 149, 61, + 27409, 36319, 36151, 11059, 52057, 5741, + 29753, 2693, 52571, 13049, 9371, 47237 +#ifndef MP_8BIT + ,39623, 207409, 128857, 37409, 141689, 189463, + 41491, 127363, 91673, 8467, 479142413, 465566323, + 961126169, 1057886029, 1222702051, 1017450739, + 1019879717, 72282697, 2048787577, 2058368051 +#endif + }; + + mp_sieve_init(&sieve); + + test_size = (int)(sizeof(to_test)/sizeof(mp_sieve_prime)); + + for (i = 0; i < test_size; i++) { + if ((e = mp_prec_small_prime(to_test[i], &ret, &sieve)) != MP_OKAY) { + fprintf(stderr,"mp_prec_small_prime failed with \"%s\" at index %d\n", + mp_error_to_string(e), i); + goto LTM_ERR; + } + if (ret != tested[i]) { + fprintf(stderr,"mp_prec_small_prime failed for %lu. Said %lu but is %lu \n", + (unsigned long)to_test[i], (unsigned long)ret, (unsigned long)tested[i]); + goto LTM_ERR; + } + } + + mp_sieve_clear(&sieve); + return EXIT_SUCCESS; +LTM_ERR: + mp_sieve_clear(&sieve); + return EXIT_FAILURE; +} + + /* Cannot test mp_exp(_d) without mp_n_root and vice versa. @@ -1879,7 +2034,11 @@ int unit_tests(int argc, char **argv) T(mp_tc_or), T(mp_tc_xor), T(s_mp_balance_mul), - T(s_mp_jacobi) + T(s_mp_jacobi), + T(mp_ilogb), + T(mp_is_small_prime), + T(mp_next_small_prime), + T(mp_prec_small_prime) #undef T }; unsigned long i; diff --git a/doc/bn.tex b/doc/bn.tex index f4bfda571..1a1125ac7 100644 --- a/doc/bn.tex +++ b/doc/bn.tex @@ -2147,6 +2147,245 @@ \subsection{To ASCII} int mp_fwrite(const mp_int *a, int radix, FILE *stream); \end{alltt} +\chapter{Factorization} + +\section{Prime Sieve} +The prime sieve is implemented as a simple segmented Sieve of Eratosthenes. It is only moderately optimized for space and runtime but should be small enough and also fast enough for almost all use-cases; quite quick for sequential access but relatively slow for random access. + +\subsection{Initialization and Clearing} +Initializing. It cannot fail because it only sets some default values. Memory is allocated later according to needs. +\index{mp\_sieve\_init} +\begin{alltt} +void mp_sieve_init(mp_sieve *sieve); +\end{alltt} +The function \texttt{mp\_sieve\_init} is equivalent to +\begin{alltt} +mp_sieve sieve = {NULL, NULL, 0}; +\end{alltt} + +Free the memory used by the sieve with +\index{mp\_sieve\_clear} +\begin{alltt} +void mp_sieve_clear(mp_sieve *sieve); +\end{alltt} + +\subsection{Primality Test of Small Numbers} +Individual small numbers can be tested for primality with +\index{mp\_is\_small\_prime} +\begin{alltt} +int mp_is_small_prime(mp_sieve_prime n, mp_sieve_prime *result, + mp_sieve *sieve); +\end{alltt} +The implementation of this function does all the heavy lifting--the building of the base sieve and the segment if one is necessary. The base sieve stays, so this function can be used to ``warm up'' the sieve and, if \texttt{n} is slightly larger than the upper limit of the base sieve, ``warm up'' the first segment, too. +It will return \texttt{MP\_SIEVE\_MAX\_REACHED} to flag the content of \texttt{result} as the last valid one. + +\subsection{Find Adjacent Primes} +To find the prime bigger than a number \texttt{n} use +\index{mp\_next\_small\_prime} +\begin{alltt} +int mp_next_small_prime(mp_sieve_prime n, mp_sieve_prime *result, + mp_sieve *sieve); +\end{alltt} +and to find the one smaller than \texttt{n} +\begin{alltt} +int mp_prec_small_prime(mp_sieve_prime n, mp_sieve_prime *result, + mp_sieve *sieve); +\end{alltt} + +\subsection{Useful Constants} +\begin{description} +\item[\texttt{MP\_SIEVE\_BIGGEST\_PRIME}] \texttt{read-only} The biggest prime the sieve can offer. It is be $65\,521$ for \texttt{MP\_8BIT}, + $4\,294\,967\,291$ for \texttt{MP\_16BIT}, \texttt{MP\_32BIT} and \texttt{MP\_64BIT}; and + $18\,446\,744\,073\,709\,551\,557$ for \texttt{MP\_64BIT} if the macro\\ + \texttt{LTM\_SIEVE\_USE\_LARGE\_SIEVE} is defined. + +\item[\texttt{mp\_sieve\_prime}] \texttt{read-only} The basic type for the numbers in the sieve. It is be \texttt{uint16\_t} for \texttt{MP\_8BIT}, \texttt{uint32\_t} for \texttt{MP\_16BIT}, \texttt{MP\_32BIT} and \texttt{MP\_64BIT}; and \texttt{uint64\_t} for \texttt{MP\_64BIT} if the macro \texttt{LTM\_SIEVE\_USE\_LARGE\_SIEVE} is defined. + +\item[\texttt{LTM\_SIEVE\_USE\_LARGE\_SIEVE}] \texttt{read-only} A flag to make a large sieve. No advantage has been seen in using 64-bit integers if available except the ability to get a sieve up to $2^64$. But in this case the base sieve gets 0.25 Gibibytes large and the segments 0.5 Gibibytes (although you can change \texttt{LTM\_SIEVE\_RANGE\_A\_B} in \texttt{bn\_mp\_sieve.c} to get smaller segments) and need a long time to fill. + + +\item[\texttt{MP\_SIEVE\_PR\_UINT}] Choses the correct macro from \texttt{inttypes.h} to print a\\ + \texttt{mp\_sieve\_prime}. The header \texttt{inttypes.h} must be included before\\ + \texttt{tommath.h} for it to work. +\end{description} + + +\subsection{Examples}\label{sec:spnexamples} +\subsubsection{Initialization and Clearing} +Using a sieve follows the same procedure as using a bigint: +\begin{alltt} +/* Declaration */ +mp_sieve sieve; + +/* + Initialization. + Cannot fail, only sets a handful of default values. + Memory allocation is done in the library itself + according to needs. + */ +mp_sieve_init(&sieve); + +/* use the sieve */ + +/* Clean up */ +mp_sieve_clear(&sieve); +\end{alltt} +\subsubsection{Primality Test} +A small program to test the input of a small number for primality. +\begin{alltt} +#include +#include +#include +/*inttypes.h is needed for printing and must be included before tommath.h*/ +#include +#include "tommath.h" +int main(int argc, char **argv) +{ + mp_sieve_prime number; + mp_sieve *base = NULL; + mp_sieve *segment = NULL; + mp_sieve_prime single_segment_a = 0; + int e; + + /* variable holding the result of mp_is_small_prime */ + mp_sieve_prime result; + + if (argc != 2) { + fprintf(stderr,"Usage %s number\textbackslash{}n", argv[0]); + exit(EXIT_FAILURE); + } + + number = (mp_sieve_prime) strtoul(argv[1],NULL, 10); + if (errno == ERANGE) { + fprintf(stderr,"strtoul(l) failed: input out of range\textbackslash{}n"); + goto LTM_ERR; + } + + mp_sieve_init(&sieve); + + if ((e = mp_is_small_prime(number, &result, &sieve)) != MP_OKAY) { + fprintf(stderr,"mp_is_small_prime failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } + + printf("The number %" LTM_SIEVE_PR_UINT " is %s prime\textbackslash{}n", + number,(result)?"":"not"); + + + mp_sieve_clear(&sieve); + exit(EXIT_SUCCESS); +LTM_ERR: + mp_sieve_clear(&sieve); + exit(EXIT_FAILURE); +} +\end{alltt} +\subsubsection{Find Adjacent Primes} +To sum up all primes up to and including \texttt{MP\_SIEVE\_BIGGEST\_PRIME} you might do something like: +\begin{alltt} +#include +#include +#include +#include +int main(int argc, char **argv) +{ + mp_sieve_prime number; + mp_sieve sieve; + mp_sieve_prime k, ret; + mp_int total, t; + int e; + + if (argc != 2) { + fprintf(stderr,"Usage %s integer\textbackslash{}n", argv[0]); + exit(EXIT_FAILURE); + } + + if ((e = mp_init_multi(&total, &t, NULL)) != MP_OKAY) { + fprintf(stderr,"mp_init_multi(segment): \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR_1; + } + errno = 0; +#if ( (defined MP_64BIT) && (defined LTM_SIEVE_USE_LARGE_SIEVE) ) + number = (mp_sieve_prime) strtoull(argv[1],NULL, 10); +#else + number = (mp_sieve_prime) strtoul(argv[1],NULL, 10); +#endif + if (errno == ERANGE) { + fprintf(stderr,"strtoul(l) failed: input out of range\textbackslash{}n"); + return EXIT_FAILURE + } + + mp_sieve_init(&sieve); + + for (k = 0, ret = 0; ret < number; k = ret) { + if ((e = mp_next_small_prime(k + 1, &ret, &sieve)) != MP_OKAY) { + if (e == LTM_SIEVE_MAX_REACHED) { +#ifdef MP_64BIT + if ((e = mp_add_d(&total, (mp_digit) k, &total)) != MP_OKAY) { + fprintf(stderr,"mp_add_d (1) failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } +#else + if ((e = mp_set_long(&t, k)) != MP_OKAY) { + fprintf(stderr,"mp_set_long (1) failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } + if ((e = mp_add(&total, &t, &total)) != MP_OKAY) { + fprintf(stderr,"mp_add (1) failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } +#endif + break; + } + fprintf(stderr,"mp_next_small_prime failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } + /* The check if the prime is below the given limit + * cannot be done in the for-loop conditions because if + * it could we wouldn't need the sieve in the first place. + */ + if (ret <= number) { +#ifdef MP_64BIT + if ((e = mp_add_d(&total, (mp_digit) k, &total)) != MP_OKAY) { + fprintf(stderr,"mp_add_d failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } +#else + if ((e = mp_set_long(&t, k)) != MP_OKAY) { + fprintf(stderr,"mp_set_long failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } + if ((e = mp_add(&total, &t, &total)) != MP_OKAY) { + fprintf(stderr,"mp_add failed: \textbackslash{}"%s\textbackslash{}"\textbackslash{}n", + mp_error_to_string(e)); + goto LTM_ERR; + } +#endif + } + } + printf("total: "); + mp_fwrite(&total,10,stdout); + putchar('\textbackslash{}n'); + + mp_clear_multi(&total, &t, NULL); + mp_sieve_clear(&sieve); + exit(EXIT_SUCCESS); +LTM_ERR: + mp_clear_multi(&total, &t, NULL); + mp_sieve_clear(&sieve); + exit(EXIT_FAILURE); +} +\end{alltt} +It took about a minute on the authors machine from 2015 to get the expected $425\,649\,736\,193\,687\,430$ for the sum of all primes up to $2^{32}$, about the same runtime as Pari/GP version 2.9.4 (with a GMP-5.1.3 kernel). + + \subsection{From ASCII} \index{mp\_read\_radix} diff --git a/libtommath_VS2008.vcproj b/libtommath_VS2008.vcproj index 7b054169c..3ab591d8b 100644 --- a/libtommath_VS2008.vcproj +++ b/libtommath_VS2008.vcproj @@ -600,10 +600,18 @@ RelativePath="bn_mp_neg.c" > + + + + @@ -724,6 +732,18 @@ RelativePath="bn_mp_shrink.c" > + + + + + + diff --git a/makefile b/makefile index 16232cbb2..9c419f114 100644 --- a/makefile +++ b/makefile @@ -37,24 +37,24 @@ bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o bn_mp bn_mp_invmod.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \ bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \ bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \ -bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \ -bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \ -bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \ -bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \ -bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o bn_mp_read_unsigned_bin.o bn_mp_reduce.o \ -bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o bn_mp_reduce_2k_setup_l.o \ -bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o bn_mp_set.o \ -bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ -bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o \ -bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o \ -bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o \ -bn_mp_toradix.o bn_mp_toradix_n.o bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o \ -bn_s_mp_add.o bn_s_mp_balance_mul.o bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o \ -bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o \ -bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o \ -bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o bn_s_mp_rand_platform.o bn_s_mp_reverse.o \ -bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o - +bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_next_small_prime.o bn_mp_or.o bn_mp_prec_small_prime.o \ +bn_mp_prime_fermat.o bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o \ +bn_mp_prime_is_prime.o bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o \ +bn_mp_prime_rabin_miller_trials.o bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o \ +bn_mp_radix_size.o bn_mp_radix_smap.o bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o \ +bn_mp_read_unsigned_bin.o bn_mp_reduce.o bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o \ +bn_mp_reduce_2k_setup_l.o bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o \ +bn_mp_set.o bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ +bn_mp_sieve.o bn_mp_sieve_clear.o bn_mp_sieve_init.o bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o \ +bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o \ +bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o \ +bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o bn_mp_toradix.o bn_mp_toradix_n.o \ +bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o bn_s_mp_add.o bn_s_mp_balance_mul.o \ +bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o \ +bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o \ +bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o \ +bn_s_mp_rand_platform.o bn_s_mp_reverse.o bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o \ +bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o #END_INS $(OBJECTS): $(HEADERS) diff --git a/makefile.mingw b/makefile.mingw index f3304912a..ac8244bd4 100644 --- a/makefile.mingw +++ b/makefile.mingw @@ -40,23 +40,24 @@ bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o bn_mp bn_mp_invmod.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \ bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \ bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \ -bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \ -bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \ -bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \ -bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \ -bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o bn_mp_read_unsigned_bin.o bn_mp_reduce.o \ -bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o bn_mp_reduce_2k_setup_l.o \ -bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o bn_mp_set.o \ -bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ -bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o \ -bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o \ -bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o \ -bn_mp_toradix.o bn_mp_toradix_n.o bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o \ -bn_s_mp_add.o bn_s_mp_balance_mul.o bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o \ -bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o \ -bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o \ -bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o bn_s_mp_rand_platform.o bn_s_mp_reverse.o \ -bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o +bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_next_small_prime.o bn_mp_or.o bn_mp_prec_small_prime.o \ +bn_mp_prime_fermat.o bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o \ +bn_mp_prime_is_prime.o bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o \ +bn_mp_prime_rabin_miller_trials.o bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o \ +bn_mp_radix_size.o bn_mp_radix_smap.o bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o \ +bn_mp_read_unsigned_bin.o bn_mp_reduce.o bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o \ +bn_mp_reduce_2k_setup_l.o bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o \ +bn_mp_set.o bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ +bn_mp_sieve.o bn_mp_sieve_clear.o bn_mp_sieve_init.o bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o \ +bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o \ +bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o \ +bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o bn_mp_toradix.o bn_mp_toradix_n.o \ +bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o bn_s_mp_add.o bn_s_mp_balance_mul.o \ +bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o \ +bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o \ +bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o \ +bn_s_mp_rand_platform.o bn_s_mp_reverse.o bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o \ +bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o HEADERS_PUB=tommath.h tommath_class.h tommath_superclass.h diff --git a/makefile.msvc b/makefile.msvc index 78030a471..ca0888351 100644 --- a/makefile.msvc +++ b/makefile.msvc @@ -32,23 +32,24 @@ bn_mp_init_copy.obj bn_mp_init_multi.obj bn_mp_init_set.obj bn_mp_init_set_int.o bn_mp_invmod.obj bn_mp_is_square.obj bn_mp_iseven.obj bn_mp_isodd.obj bn_mp_kronecker.obj bn_mp_lcm.obj bn_mp_lshd.obj \ bn_mp_mod.obj bn_mp_mod_2d.obj bn_mp_mod_d.obj bn_mp_montgomery_calc_normalization.obj bn_mp_montgomery_reduce.obj \ bn_mp_montgomery_setup.obj bn_mp_mul.obj bn_mp_mul_2.obj bn_mp_mul_2d.obj bn_mp_mul_d.obj bn_mp_mulmod.obj \ -bn_mp_n_root.obj bn_mp_n_root_ex.obj bn_mp_neg.obj bn_mp_or.obj bn_mp_prime_fermat.obj \ -bn_mp_prime_frobenius_underwood.obj bn_mp_prime_is_divisible.obj bn_mp_prime_is_prime.obj \ -bn_mp_prime_miller_rabin.obj bn_mp_prime_next_prime.obj bn_mp_prime_rabin_miller_trials.obj \ -bn_mp_prime_rand.obj bn_mp_prime_strong_lucas_selfridge.obj bn_mp_radix_size.obj bn_mp_radix_smap.obj \ -bn_mp_rand.obj bn_mp_read_radix.obj bn_mp_read_signed_bin.obj bn_mp_read_unsigned_bin.obj bn_mp_reduce.obj \ -bn_mp_reduce_2k.obj bn_mp_reduce_2k_l.obj bn_mp_reduce_2k_setup.obj bn_mp_reduce_2k_setup_l.obj \ -bn_mp_reduce_is_2k.obj bn_mp_reduce_is_2k_l.obj bn_mp_reduce_setup.obj bn_mp_rshd.obj bn_mp_set.obj \ -bn_mp_set_double.obj bn_mp_set_int.obj bn_mp_set_long.obj bn_mp_set_long_long.obj bn_mp_shrink.obj \ -bn_mp_signed_bin_size.obj bn_mp_sqr.obj bn_mp_sqrmod.obj bn_mp_sqrt.obj bn_mp_sqrtmod_prime.obj bn_mp_sub.obj \ -bn_mp_sub_d.obj bn_mp_submod.obj bn_mp_tc_and.obj bn_mp_tc_div_2d.obj bn_mp_tc_or.obj bn_mp_tc_xor.obj \ -bn_mp_to_signed_bin.obj bn_mp_to_signed_bin_n.obj bn_mp_to_unsigned_bin.obj bn_mp_to_unsigned_bin_n.obj \ -bn_mp_toradix.obj bn_mp_toradix_n.obj bn_mp_unsigned_bin_size.obj bn_mp_xor.obj bn_mp_zero.obj bn_prime_tab.obj \ -bn_s_mp_add.obj bn_s_mp_balance_mul.obj bn_s_mp_exptmod.obj bn_s_mp_exptmod_fast.obj bn_s_mp_get_bit.obj \ -bn_s_mp_invmod_fast.obj bn_s_mp_invmod_slow.obj bn_s_mp_karatsuba_mul.obj bn_s_mp_karatsuba_sqr.obj \ -bn_s_mp_montgomery_reduce_fast.obj bn_s_mp_mul_digs.obj bn_s_mp_mul_digs_fast.obj bn_s_mp_mul_high_digs.obj \ -bn_s_mp_mul_high_digs_fast.obj bn_s_mp_rand_jenkins.obj bn_s_mp_rand_platform.obj bn_s_mp_reverse.obj \ -bn_s_mp_sqr.obj bn_s_mp_sqr_fast.obj bn_s_mp_sub.obj bn_s_mp_toom_mul.obj bn_s_mp_toom_sqr.obj +bn_mp_n_root.obj bn_mp_n_root_ex.obj bn_mp_neg.obj bn_mp_next_small_prime.obj bn_mp_or.obj bn_mp_prec_small_prime.obj \ +bn_mp_prime_fermat.obj bn_mp_prime_frobenius_underwood.obj bn_mp_prime_is_divisible.obj \ +bn_mp_prime_is_prime.obj bn_mp_prime_miller_rabin.obj bn_mp_prime_next_prime.obj \ +bn_mp_prime_rabin_miller_trials.obj bn_mp_prime_rand.obj bn_mp_prime_strong_lucas_selfridge.obj \ +bn_mp_radix_size.obj bn_mp_radix_smap.obj bn_mp_rand.obj bn_mp_read_radix.obj bn_mp_read_signed_bin.obj \ +bn_mp_read_unsigned_bin.obj bn_mp_reduce.obj bn_mp_reduce_2k.obj bn_mp_reduce_2k_l.obj bn_mp_reduce_2k_setup.obj \ +bn_mp_reduce_2k_setup_l.obj bn_mp_reduce_is_2k.obj bn_mp_reduce_is_2k_l.obj bn_mp_reduce_setup.obj bn_mp_rshd.obj \ +bn_mp_set.obj bn_mp_set_double.obj bn_mp_set_int.obj bn_mp_set_long.obj bn_mp_set_long_long.obj bn_mp_shrink.obj \ +bn_mp_sieve.obj bn_mp_sieve_clear.obj bn_mp_sieve_init.obj bn_mp_signed_bin_size.obj bn_mp_sqr.obj bn_mp_sqrmod.obj \ +bn_mp_sqrt.obj bn_mp_sqrtmod_prime.obj bn_mp_sub.obj bn_mp_sub_d.obj bn_mp_submod.obj bn_mp_tc_and.obj \ +bn_mp_tc_div_2d.obj bn_mp_tc_or.obj bn_mp_tc_xor.obj bn_mp_to_signed_bin.obj bn_mp_to_signed_bin_n.obj \ +bn_mp_to_unsigned_bin.obj bn_mp_to_unsigned_bin_n.obj bn_mp_toradix.obj bn_mp_toradix_n.obj \ +bn_mp_unsigned_bin_size.obj bn_mp_xor.obj bn_mp_zero.obj bn_prime_tab.obj bn_s_mp_add.obj bn_s_mp_balance_mul.obj \ +bn_s_mp_exptmod.obj bn_s_mp_exptmod_fast.obj bn_s_mp_get_bit.obj bn_s_mp_invmod_fast.obj bn_s_mp_invmod_slow.obj \ +bn_s_mp_karatsuba_mul.obj bn_s_mp_karatsuba_sqr.obj bn_s_mp_montgomery_reduce_fast.obj bn_s_mp_mul_digs.obj \ +bn_s_mp_mul_digs_fast.obj bn_s_mp_mul_high_digs.obj bn_s_mp_mul_high_digs_fast.obj bn_s_mp_rand_jenkins.obj \ +bn_s_mp_rand_platform.obj bn_s_mp_reverse.obj bn_s_mp_sqr.obj bn_s_mp_sqr_fast.obj bn_s_mp_sub.obj \ +bn_s_mp_toom_mul.obj bn_s_mp_toom_sqr.obj HEADERS_PUB=tommath.h tommath_class.h tommath_superclass.h diff --git a/makefile.shared b/makefile.shared index b4be47c33..7d5059747 100644 --- a/makefile.shared +++ b/makefile.shared @@ -34,23 +34,24 @@ bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o bn_mp bn_mp_invmod.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \ bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \ bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \ -bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \ -bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \ -bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \ -bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \ -bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o bn_mp_read_unsigned_bin.o bn_mp_reduce.o \ -bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o bn_mp_reduce_2k_setup_l.o \ -bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o bn_mp_set.o \ -bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ -bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o \ -bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o \ -bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o \ -bn_mp_toradix.o bn_mp_toradix_n.o bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o \ -bn_s_mp_add.o bn_s_mp_balance_mul.o bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o \ -bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o \ -bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o \ -bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o bn_s_mp_rand_platform.o bn_s_mp_reverse.o \ -bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o +bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_next_small_prime.o bn_mp_or.o bn_mp_prec_small_prime.o \ +bn_mp_prime_fermat.o bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o \ +bn_mp_prime_is_prime.o bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o \ +bn_mp_prime_rabin_miller_trials.o bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o \ +bn_mp_radix_size.o bn_mp_radix_smap.o bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o \ +bn_mp_read_unsigned_bin.o bn_mp_reduce.o bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o \ +bn_mp_reduce_2k_setup_l.o bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o \ +bn_mp_set.o bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ +bn_mp_sieve.o bn_mp_sieve_clear.o bn_mp_sieve_init.o bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o \ +bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o \ +bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o \ +bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o bn_mp_toradix.o bn_mp_toradix_n.o \ +bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o bn_s_mp_add.o bn_s_mp_balance_mul.o \ +bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o \ +bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o \ +bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o \ +bn_s_mp_rand_platform.o bn_s_mp_reverse.o bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o \ +bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o #END_INS diff --git a/makefile.unix b/makefile.unix index 6e097cf59..6016fdb5b 100644 --- a/makefile.unix +++ b/makefile.unix @@ -41,23 +41,24 @@ bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o bn_mp_init_set_int.o bn_mp bn_mp_invmod.o bn_mp_is_square.o bn_mp_iseven.o bn_mp_isodd.o bn_mp_kronecker.o bn_mp_lcm.o bn_mp_lshd.o \ bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mod_d.o bn_mp_montgomery_calc_normalization.o bn_mp_montgomery_reduce.o \ bn_mp_montgomery_setup.o bn_mp_mul.o bn_mp_mul_2.o bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_mulmod.o \ -bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_or.o bn_mp_prime_fermat.o \ -bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o bn_mp_prime_is_prime.o \ -bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o bn_mp_prime_rabin_miller_trials.o \ -bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o bn_mp_radix_size.o bn_mp_radix_smap.o \ -bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o bn_mp_read_unsigned_bin.o bn_mp_reduce.o \ -bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o bn_mp_reduce_2k_setup_l.o \ -bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o bn_mp_set.o \ -bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ -bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o \ -bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o \ -bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o \ -bn_mp_toradix.o bn_mp_toradix_n.o bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o \ -bn_s_mp_add.o bn_s_mp_balance_mul.o bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o \ -bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o \ -bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o \ -bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o bn_s_mp_rand_platform.o bn_s_mp_reverse.o \ -bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o +bn_mp_n_root.o bn_mp_n_root_ex.o bn_mp_neg.o bn_mp_next_small_prime.o bn_mp_or.o bn_mp_prec_small_prime.o \ +bn_mp_prime_fermat.o bn_mp_prime_frobenius_underwood.o bn_mp_prime_is_divisible.o \ +bn_mp_prime_is_prime.o bn_mp_prime_miller_rabin.o bn_mp_prime_next_prime.o \ +bn_mp_prime_rabin_miller_trials.o bn_mp_prime_rand.o bn_mp_prime_strong_lucas_selfridge.o \ +bn_mp_radix_size.o bn_mp_radix_smap.o bn_mp_rand.o bn_mp_read_radix.o bn_mp_read_signed_bin.o \ +bn_mp_read_unsigned_bin.o bn_mp_reduce.o bn_mp_reduce_2k.o bn_mp_reduce_2k_l.o bn_mp_reduce_2k_setup.o \ +bn_mp_reduce_2k_setup_l.o bn_mp_reduce_is_2k.o bn_mp_reduce_is_2k_l.o bn_mp_reduce_setup.o bn_mp_rshd.o \ +bn_mp_set.o bn_mp_set_double.o bn_mp_set_int.o bn_mp_set_long.o bn_mp_set_long_long.o bn_mp_shrink.o \ +bn_mp_sieve.o bn_mp_sieve_clear.o bn_mp_sieve_init.o bn_mp_signed_bin_size.o bn_mp_sqr.o bn_mp_sqrmod.o \ +bn_mp_sqrt.o bn_mp_sqrtmod_prime.o bn_mp_sub.o bn_mp_sub_d.o bn_mp_submod.o bn_mp_tc_and.o \ +bn_mp_tc_div_2d.o bn_mp_tc_or.o bn_mp_tc_xor.o bn_mp_to_signed_bin.o bn_mp_to_signed_bin_n.o \ +bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o bn_mp_toradix.o bn_mp_toradix_n.o \ +bn_mp_unsigned_bin_size.o bn_mp_xor.o bn_mp_zero.o bn_prime_tab.o bn_s_mp_add.o bn_s_mp_balance_mul.o \ +bn_s_mp_exptmod.o bn_s_mp_exptmod_fast.o bn_s_mp_get_bit.o bn_s_mp_invmod_fast.o bn_s_mp_invmod_slow.o \ +bn_s_mp_karatsuba_mul.o bn_s_mp_karatsuba_sqr.o bn_s_mp_montgomery_reduce_fast.o bn_s_mp_mul_digs.o \ +bn_s_mp_mul_digs_fast.o bn_s_mp_mul_high_digs.o bn_s_mp_mul_high_digs_fast.o bn_s_mp_rand_jenkins.o \ +bn_s_mp_rand_platform.o bn_s_mp_reverse.o bn_s_mp_sqr.o bn_s_mp_sqr_fast.o bn_s_mp_sub.o \ +bn_s_mp_toom_mul.o bn_s_mp_toom_sqr.o HEADERS_PUB=tommath.h tommath_class.h tommath_superclass.h diff --git a/tommath.h b/tommath.h index d70eef659..174500020 100644 --- a/tommath.h +++ b/tommath.h @@ -125,7 +125,8 @@ typedef enum { MP_ERR = -1, MP_MEM = -2, MP_VAL = -3, - MP_ITER = -4 + MP_ITER = -4, + MP_SIEVE_MAX_REACHED = -5 } mp_err; #else typedef int mp_sign; @@ -145,6 +146,7 @@ typedef int mp_err; #define MP_VAL -3 /* invalid input */ #define MP_RANGE (MP_DEPRECATED_PRAGMA("MP_RANGE has been deprecated in favor of MP_VAL") MP_VAL) #define MP_ITER -4 /* Max. iterations reached */ +#define MP_SIEVE_MAX_REACHED -5 /* Ret. for. max. poss. prime found in mp_next_small_prime */ #endif /* tunable cutoffs */ @@ -562,6 +564,37 @@ mp_err mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y) /* ---> Primes <--- */ +/* + * All MP_xBIT sizes need one data type that has twice the size of "x", + * that means that the type uint16_t must be available, even for MP_8BIT. + */ +#ifdef MP_8BIT +typedef uint16_t mp_sieve_prime; +# define MP_SIEVE_BIGGEST_PRIME 65521lu +# define MP_SIEVE_PR_UINT PRIu16 +#elif ( (defined MP_64BIT) && (defined MP_SIEVE_USE_LARGE_SIEVE) ) +typedef uint64_t mp_sieve_prime; +# define MP_SIEVE_BIGGEST_PRIME 18446744073709551557llu +# define MP_SIEVE_PR_UINT PRIu64 +#else +typedef uint32_t mp_sieve_prime; +# define MP_SIEVE_BIGGEST_PRIME 4294967291lu +# define MP_SIEVE_PR_UINT PRIu32 +#endif + +typedef struct mp_single_sieve_t { + mp_sieve_prime *content; /* bitset holding the sieve */ + mp_sieve_prime size; /* number of entries (which is a slightly misleading description) */ + mp_sieve_prime alloc; /* size in bytes */ +} mp_single_sieve; + +typedef struct mp_sieve_t { + mp_single_sieve base; /* base sieve (0..MP_SIEVE_PRIME_MAX_SQRT) */ + mp_single_sieve segment; /* segment (range_a_b) */ + mp_sieve_prime single_segment_a; /* startpoint of segment */ +} mp_sieve; + + /* number of primes */ #ifdef MP_8BIT # define MP_PRIME_SIZE 31 @@ -617,6 +650,7 @@ mp_err mp_prime_frobenius_underwood(const mp_int *N, mp_bool *result) MP_WUR; */ mp_err mp_prime_is_prime(const mp_int *a, int t, mp_bool *result) MP_WUR; + /* finds the next prime after the number "a" using "t" trials * of Miller-Rabin. * @@ -652,9 +686,34 @@ MP_DEPRECATED(mp_prime_rand) mp_err mp_prime_random_ex(mp_int *a, int t, int siz private_mp_prime_callback cb, void *dat) MP_WUR; mp_err mp_prime_rand(mp_int *a, int t, int size, int flags) MP_WUR; + /* Integer logarithm to integer base */ mp_err mp_ilogb(const mp_int *a, mp_digit base, mp_int *c) MP_WUR; +/* Init a sieve. Sets the necessary defaults */ +void mp_sieve_init(mp_sieve *sieve); + +/* int mp_sieve_init(mp_sieve *sieve); */ +/* Clear a sieve. Frees the memory for the contents of base and segment */ +void mp_sieve_clear(mp_sieve *sieve); + +/* + Deterministicaly checks if a small prime (< MP_SIEVE_PRIME_MAX) is prime. + Uses the sieve so is not the fastest for random checks but quick at + linear access. + */ +mp_err mp_is_small_prime(mp_sieve_prime n, mp_sieve_prime *result, mp_sieve *sieve) MP_WUR; +/* + Puts the next prime >= n in "result". May return MP_SIEVE_MAX_REACHED to flag the content + of "result" as the last valid one. +*/ +mp_err mp_next_small_prime(mp_sieve_prime n, mp_sieve_prime *result, mp_sieve *sieve) MP_WUR; +/* + Puts the prime preceeding n in "result" +*/ +mp_err mp_prec_small_prime(mp_sieve_prime n, mp_sieve_prime *result, mp_sieve *sieve) MP_WUR; + + /* ---> radix conversion <--- */ int mp_count_bits(const mp_int *a) MP_WUR; diff --git a/tommath_class.h b/tommath_class.h index 880ea8bd2..9e8f23497 100644 --- a/tommath_class.h +++ b/tommath_class.h @@ -82,7 +82,9 @@ # define BN_MP_N_ROOT_C # define BN_MP_N_ROOT_EX_C # define BN_MP_NEG_C +# define BN_MP_NEXT_SMALL_PRIME_C # define BN_MP_OR_C +# define BN_MP_PREC_SMALL_PRIME_C # define BN_MP_PRIME_FERMAT_C # define BN_MP_PRIME_FROBENIUS_UNDERWOOD_C # define BN_MP_PRIME_IS_DIVISIBLE_C @@ -113,6 +115,9 @@ # define BN_MP_SET_LONG_C # define BN_MP_SET_LONG_LONG_C # define BN_MP_SHRINK_C +# define BN_MP_SIEVE_C +# define BN_MP_SIEVE_CLEAR_C +# define BN_MP_SIEVE_INIT_C # define BN_MP_SIGNED_BIN_SIZE_C # define BN_MP_SQR_C # define BN_MP_SQRMOD_C @@ -631,6 +636,10 @@ # define BN_MP_COPY_C #endif +#if defined(BN_MP_NEXT_SMALL_PRIME_C) +# define BN_MP_IS_SMALL_PRIME_C +#endif + #if defined(BN_MP_OR_C) # define BN_MP_INIT_COPY_C # define BN_MP_CLAMP_C @@ -638,6 +647,10 @@ # define BN_MP_CLEAR_C #endif +#if defined(BN_MP_PREC_SMALL_PRIME_C) +# define BN_MP_IS_SMALL_PRIME_C +#endif + #if defined(BN_MP_PRIME_FERMAT_C) # define BN_MP_CMP_D_C # define BN_MP_INIT_C @@ -883,6 +896,25 @@ #if defined(BN_MP_SHRINK_C) #endif +#if defined(BN_MP_SIEVE_C) +# define BN_S_MP_SIEVE_SETALL_C +# define BN_S_MP_SIEVE_CLEAR_C +# define BN_S_MP_SIEVE_GET_C +# define BN_S_MP_SIEVE_NEXTSET_C +# define BN_S_MP_ERATOSTHENES_C +# define BN_S_MP_ERATOSTHENES_INIT_C +# define BN_S_MP_ERATOSTHENES_SEGMENT_C +# define BN_S_MP_ERATOSTHENES_SEGMENT_INIT_C +# define BN_S_MP_ERATOSTHENES_SEGMENT_CLEAR_C +# define BN_MP_IS_SMALL_PRIME_C +#endif + +#if defined(BN_MP_SIEVE_CLEAR_C) +#endif + +#if defined(BN_MP_SIEVE_INIT_C) +#endif + #if defined(BN_MP_SIGNED_BIN_SIZE_C) # define BN_MP_UNSIGNED_BIN_SIZE_C #endif diff --git a/tommath_private.h b/tommath_private.h index 3bf4e7c43..48ff4f2fa 100644 --- a/tommath_private.h +++ b/tommath_private.h @@ -127,6 +127,18 @@ extern void *MP_CALLOC(size_t nmemb, size_t size); extern void MP_FREE(void *mem, size_t size); #endif +/* Size of the base sieve of mp_sieve*/ +#ifdef MP_8BIT +# define MP_SIEVE_PRIME_MAX 0xFFFFlu +# define MP_SIEVE_PRIME_MAX_SQRT 0xFFlu +#elif ( (defined MP_64BIT) && (defined MP_SIEVE_USE_LARGE_SIEVE) ) +# define MP_SIEVE_PRIME_MAX 0xFFFFFFFFFFFFFFFFllu +# define MP_SIEVE_PRIME_MAX_SQRT 0xFFFFFFFFllu +#else +# define MP_SIEVE_PRIME_MAX 0xFFFFFFFFlu +# define MP_SIEVE_PRIME_MAX_SQRT 0xFFFFlu +#endif + /* TODO: Remove PRIVATE_MP_WARRAY as soon as deprecated MP_WARRAY is removed from tommath.h */ #undef MP_WARRAY #define MP_WARRAY PRIVATE_MP_WARRAY