AuroraMiddleware/libtommath
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

simplifications: prime functions

Browse Source
This commit is contained in:
Daniel Mendler 2019-10-29 20:07:29 +01:00
parent 8ac493512c
commit 448f35e2e1
No known key found for this signature in database
GPG Key ID: D88ADB2A2693CA43
9 changed files with 85 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
mp_prime_fermat.c
View File

@ -16,9 +16,6 @@ mp_err mp_prime_fermat(const mp_int *a, const mp_int *b, bool *result)
mp_int t; mp_int t;
mp_err err; mp_err err;
/* default to composite */
*result = false;
/* ensure b > 1 */ /* ensure b > 1 */
if (mp_cmp_d(b, 1uL) != MP_GT) { if (mp_cmp_d(b, 1uL) != MP_GT) {
return MP_VAL; return MP_VAL;
@ -31,16 +28,13 @@ mp_err mp_prime_fermat(const mp_int *a, const mp_int *b, bool *result)
/* compute t = b**a mod a */ /* compute t = b**a mod a */
if ((err = mp_exptmod(b, a, a, &t)) != MP_OKAY) { if ((err = mp_exptmod(b, a, a, &t)) != MP_OKAY) {
goto LBL_T; goto LBL_ERR;
} }
/* is it equal to b? */ /* is it equal to b? */
if (mp_cmp(&t, b) == MP_EQ) { *result = mp_cmp(&t, b) == MP_EQ;
*result = true;
}
err = MP_OKAY; LBL_ERR:
LBL_T:
mp_clear(&t); mp_clear(&t);
return err; return err;
} }

66
mp_prime_frobenius_underwood.c
View File

@ -23,17 +23,16 @@
mp_err mp_prime_frobenius_underwood(const mp_int *N, bool *result) mp_err mp_prime_frobenius_underwood(const mp_int *N, bool *result)
{ {
mp_int T1z, T2z, Np1z, sz, tz; mp_int T1z, T2z, Np1z, sz, tz;
int a, ap2, i;
int a, ap2, length, i, j;
mp_err err; mp_err err;
*result = false;
if ((err = mp_init_multi(&T1z, &T2z, &Np1z, &sz, &tz, NULL)) != MP_OKAY) { if ((err = mp_init_multi(&T1z, &T2z, &Np1z, &sz, &tz, NULL)) != MP_OKAY) {
return err; return err;
} }
for (a = 0; a < LTM_FROBENIUS_UNDERWOOD_A; a++) { for (a = 0; a < LTM_FROBENIUS_UNDERWOOD_A; a++) {
int j;
/* TODO: That's ugly! No, really, it is! */ /* TODO: That's ugly! No, really, it is! */
if ((a==2) || (a==4) || (a==7) || (a==8) || (a==10) || if ((a==2) || (a==4) || (a==7) || (a==8) || (a==10) ||
(a==14) || (a==18) || (a==23) || (a==26) || (a==28)) { (a==14) || (a==18) || (a==23) || (a==26) || (a==28)) {
@ -42,7 +41,7 @@ mp_err mp_prime_frobenius_underwood(const mp_int *N, bool *result)
mp_set_i32(&T1z, (int32_t)((a * a) - 4)); mp_set_i32(&T1z, (int32_t)((a * a) - 4));
if ((err = mp_kronecker(&T1z, N, &j)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_kronecker(&T1z, N, &j)) != MP_OKAY) goto LBL_END;
if (j == -1) { if (j == -1) {
break; break;
@ -50,73 +49,76 @@ mp_err mp_prime_frobenius_underwood(const mp_int *N, bool *result)
if (j == 0) { if (j == 0) {
/* composite */ /* composite */
goto LBL_FU_ERR; *result = false;
goto LBL_END;
} }
} }
/* Tell it a composite and set return value accordingly */ /* Tell it a composite and set return value accordingly */
if (a >= LTM_FROBENIUS_UNDERWOOD_A) { if (a >= LTM_FROBENIUS_UNDERWOOD_A) {
err = MP_ITER; err = MP_ITER;
goto LBL_FU_ERR; goto LBL_END;
} }
/* Composite if N and (a+4)*(2*a+5) are not coprime */ /* Composite if N and (a+4)*(2*a+5) are not coprime */
mp_set_u32(&T1z, (uint32_t)((a+4)*((2*a)+5))); mp_set_u32(&T1z, (uint32_t)((a+4)*((2*a)+5)));
if ((err = mp_gcd(N, &T1z, &T1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_gcd(N, &T1z, &T1z)) != MP_OKAY) goto LBL_END;
if (!((T1z.used == 1) && (T1z.dp[0] == 1u))) goto LBL_FU_ERR; if (!((T1z.used == 1) && (T1z.dp[0] == 1u))) {
/* composite */
*result = false;
goto LBL_END;
}
ap2 = a + 2; ap2 = a + 2;
if ((err = mp_add_d(N, 1uL, &Np1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_add_d(N, 1uL, &Np1z)) != MP_OKAY) goto LBL_END;
mp_set(&sz, 1uL); mp_set(&sz, 1uL);
mp_set(&tz, 2uL); mp_set(&tz, 2uL);
length = mp_count_bits(&Np1z);
for (i = length - 2; i >= 0; i--) { for (i = mp_count_bits(&Np1z) - 2; i >= 0; i--) {
/* /*
* temp = (sz*(a*sz+2*tz))%N; * temp = (sz*(a*sz+2*tz))%N;
* tz = ((tz-sz)*(tz+sz))%N; * tz = ((tz-sz)*(tz+sz))%N;
* sz = temp; * sz = temp;
*/ */
if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY) goto LBL_END;
/* a = 0 at about 50% of the cases (non-square and odd input) */ /* a = 0 at about 50% of the cases (non-square and odd input) */
if (a != 0) { if (a != 0) {
if ((err = mp_mul_d(&sz, (mp_digit)a, &T1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mul_d(&sz, (mp_digit)a, &T1z)) != MP_OKAY) goto LBL_END;
if ((err = mp_add(&T1z, &T2z, &T2z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_add(&T1z, &T2z, &T2z)) != MP_OKAY) goto LBL_END;
} }
if ((err = mp_mul(&T2z, &sz, &T1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mul(&T2z, &sz, &T1z)) != MP_OKAY) goto LBL_END;
if ((err = mp_sub(&tz, &sz, &T2z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_sub(&tz, &sz, &T2z)) != MP_OKAY) goto LBL_END;
if ((err = mp_add(&sz, &tz, &sz)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_add(&sz, &tz, &sz)) != MP_OKAY) goto LBL_END;
if ((err = mp_mul(&sz, &T2z, &tz)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mul(&sz, &T2z, &tz)) != MP_OKAY) goto LBL_END;
if ((err = mp_mod(&tz, N, &tz)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mod(&tz, N, &tz)) != MP_OKAY) goto LBL_END;
if ((err = mp_mod(&T1z, N, &sz)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mod(&T1z, N, &sz)) != MP_OKAY) goto LBL_END;
if (s_mp_get_bit(&Np1z, (unsigned int)i)) { if (s_mp_get_bit(&Np1z, i)) {
/* /*
* temp = (a+2) * sz + tz * temp = (a+2) * sz + tz
* tz = 2 * tz - sz * tz = 2 * tz - sz
* sz = temp * sz = temp
*/ */
if (a == 0) { if (a == 0) {
if ((err = mp_mul_2(&sz, &T1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mul_2(&sz, &T1z)) != MP_OKAY) goto LBL_END;
} else { } else {
if ((err = mp_mul_d(&sz, (mp_digit)ap2, &T1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mul_d(&sz, (mp_digit)ap2, &T1z)) != MP_OKAY) goto LBL_END;
} }
if ((err = mp_add(&T1z, &tz, &T1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_add(&T1z, &tz, &T1z)) != MP_OKAY) goto LBL_END;
if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mul_2(&tz, &T2z)) != MP_OKAY) goto LBL_END;
if ((err = mp_sub(&T2z, &sz, &tz)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_sub(&T2z, &sz, &tz)) != MP_OKAY) goto LBL_END;
mp_exch(&sz, &T1z); mp_exch(&sz, &T1z);
} }
} }
mp_set_u32(&T1z, (uint32_t)((2 * a) + 5)); mp_set_u32(&T1z, (uint32_t)((2 * a) + 5));
if ((err = mp_mod(&T1z, N, &T1z)) != MP_OKAY) goto LBL_FU_ERR; if ((err = mp_mod(&T1z, N, &T1z)) != MP_OKAY) goto LBL_END;
if (mp_iszero(&sz) && (mp_cmp(&tz, &T1z) == MP_EQ)) {
*result = true;
}
LBL_FU_ERR: *result = mp_iszero(&sz) && (mp_cmp(&tz, &T1z) == MP_EQ);
LBL_END:
mp_clear_multi(&tz, &sz, &Np1z, &T2z, &T1z, NULL); mp_clear_multi(&tz, &sz, &Np1z, &T2z, &T1z, NULL);
return err; return err;
} }

14
mp_prime_is_prime.c
View File

@ -13,14 +13,12 @@ static unsigned int s_floor_ilog2(int value)
return r; return r;
} }
mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result) mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
{ {
mp_int b; mp_int b;
int ix, p_max = 0, size_a, len; int ix;
bool res; bool res;
mp_err err; mp_err err;
unsigned int fips_rand, mask;
/* default to no */ /* default to no */
*result = false; *result = false;
@ -133,6 +131,8 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
TODO: can be made a bit finer grained but comparing is not free. TODO: can be made a bit finer grained but comparing is not free.
*/ */
if (t < 0) { if (t < 0) {
int p_max = 0;
/* /*
Sorenson, Jonathan; Webster, Jonathan (2015). Sorenson, Jonathan; Webster, Jonathan (2015).
"Strong Pseudoprimes to Twelve Prime Bases". "Strong Pseudoprimes to Twelve Prime Bases".
@ -174,6 +174,9 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
See Fips 186.4 p. 126ff See Fips 186.4 p. 126ff
*/ */
else if (t > 0) { else if (t > 0) {
unsigned int mask;
int size_a;
/* /*
* The mp_digit's have a defined bit-size but the size of the * The mp_digit's have a defined bit-size but the size of the
* array a.dp is a simple 'int' and this library can not assume full * array a.dp is a simple 'int' and this library can not assume full
@ -219,6 +222,9 @@ mp_err mp_prime_is_prime(const mp_int *a, int t, bool *result)
need to be prime. need to be prime.
*/ */
for (ix = 0; ix < t; ix++) { for (ix = 0; ix < t; ix++) {
unsigned int fips_rand;
int len;
/* mp_rand() guarantees the first digit to be non-zero */ /* mp_rand() guarantees the first digit to be non-zero */
if ((err = mp_rand(&b, 1)) != MP_OKAY) { if ((err = mp_rand(&b, 1)) != MP_OKAY) {
goto LBL_B; goto LBL_B;

28
mp_prime_miller_rabin.c
View File

@ -16,9 +16,6 @@ mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result)
mp_err err; mp_err err;
int s, j; int s, j;
/* default */
*result = false;
/* ensure b > 1 */ /* ensure b > 1 */
if (mp_cmp_d(b, 1uL) != MP_GT) { if (mp_cmp_d(b, 1uL) != MP_GT) {
return MP_VAL; return MP_VAL;
@ -29,12 +26,12 @@ mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result)
return err; return err;
} }
if ((err = mp_sub_d(&n1, 1uL, &n1)) != MP_OKAY) { if ((err = mp_sub_d(&n1, 1uL, &n1)) != MP_OKAY) {
goto LBL_N1; goto LBL_ERR1;
} }
/* set 2**s * r = n1 */ /* set 2**s * r = n1 */
if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) { if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) {
goto LBL_N1; goto LBL_ERR1;
} }
/* count the number of least significant bits /* count the number of least significant bits
@ -44,15 +41,15 @@ mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result)
/* now divide n - 1 by 2**s */ /* now divide n - 1 by 2**s */
if ((err = mp_div_2d(&r, s, &r, NULL)) != MP_OKAY) { if ((err = mp_div_2d(&r, s, &r, NULL)) != MP_OKAY) {
goto LBL_R; goto LBL_ERR2;
} }
/* compute y = b**r mod a */ /* compute y = b**r mod a */
if ((err = mp_init(&y)) != MP_OKAY) { if ((err = mp_init(&y)) != MP_OKAY) {
goto LBL_R; goto LBL_ERR2;
} }
if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) { if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) {
goto LBL_Y; goto LBL_END;
} }
/* if y != 1 and y != n1 do */ /* if y != 1 and y != n1 do */
@ -61,12 +58,13 @@ mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result)
/* while j <= s-1 and y != n1 */ /* while j <= s-1 and y != n1 */
while ((j <= (s - 1)) && (mp_cmp(&y, &n1) != MP_EQ)) { while ((j <= (s - 1)) && (mp_cmp(&y, &n1) != MP_EQ)) {
if ((err = mp_sqrmod(&y, a, &y)) != MP_OKAY) { if ((err = mp_sqrmod(&y, a, &y)) != MP_OKAY) {
goto LBL_Y; goto LBL_END;
} }
/* if y == 1 then composite */ /* if y == 1 then composite */
if (mp_cmp_d(&y, 1uL) == MP_EQ) { if (mp_cmp_d(&y, 1uL) == MP_EQ) {
goto LBL_Y; *result = false;
goto LBL_END;
} }
++j; ++j;
@ -74,17 +72,19 @@ mp_err mp_prime_miller_rabin(const mp_int *a, const mp_int *b, bool *result)
/* if y != n1 then composite */ /* if y != n1 then composite */
if (mp_cmp(&y, &n1) != MP_EQ) { if (mp_cmp(&y, &n1) != MP_EQ) {
goto LBL_Y; *result = false;
goto LBL_END;
} }
} }
/* probably prime now */ /* probably prime now */
*result = true; *result = true;
LBL_Y:
LBL_END:
mp_clear(&y); mp_clear(&y);
LBL_R: LBL_ERR2:
mp_clear(&r); mp_clear(&r);
LBL_N1: LBL_ERR1:
mp_clear(&n1); mp_clear(&n1);
return err; return err;
} }

27
mp_prime_next_prime.c
View File

@ -10,11 +10,10 @@
*/ */
mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style) mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
{ {
int x, y; int x;
mp_ord cmp;
mp_err err; mp_err err;
bool res = false; bool res = false;
mp_digit res_tab[MP_PRIME_TAB_SIZE], step, kstep; mp_digit res_tab[MP_PRIME_TAB_SIZE], kstep;
mp_int b; mp_int b;
/* force positive */ /* force positive */
@ -24,7 +23,7 @@ mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
if (mp_cmp_d(a, s_mp_prime_tab[MP_PRIME_TAB_SIZE-1]) == MP_LT) { if (mp_cmp_d(a, s_mp_prime_tab[MP_PRIME_TAB_SIZE-1]) == MP_LT) {
/* find which prime it is bigger than "a" */ /* find which prime it is bigger than "a" */
for (x = 0; x < MP_PRIME_TAB_SIZE; x++) { for (x = 0; x < MP_PRIME_TAB_SIZE; x++) {
cmp = mp_cmp_d(a, s_mp_prime_tab[x]); mp_ord cmp = mp_cmp_d(a, s_mp_prime_tab[x]);
if (cmp == MP_EQ) { if (cmp == MP_EQ) {
continue; continue;
} }
@ -42,11 +41,7 @@ mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
} }
/* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */ /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
if (bbs_style) { kstep = bbs_style ? 4 : 2;
kstep = 4;
} else {
kstep = 2;
}
/* at this point we will use a combination of a sieve and Miller-Rabin */ /* at this point we will use a combination of a sieve and Miller-Rabin */
@ -79,11 +74,12 @@ mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
} }
for (;;) { for (;;) {
mp_digit step = 0;
bool y;
/* skip to the next non-trivially divisible candidate */ /* skip to the next non-trivially divisible candidate */
step = 0;
do { do {
/* y == 1 if any residue was zero [e.g. cannot be prime] */ /* y == true if any residue was zero [e.g. cannot be prime] */
y = 0; y = false;
/* increase step to next candidate */ /* increase step to next candidate */
step += kstep; step += kstep;
@ -100,10 +96,10 @@ mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
/* set flag if zero */ /* set flag if zero */
if (res_tab[x] == 0u) { if (res_tab[x] == 0u) {
y = 1; y = true;
} }
} }
} while ((y == 1) && (step < (((mp_digit)1 << MP_DIGIT_BIT) - kstep))); } while (y && (step < (((mp_digit)1 << MP_DIGIT_BIT) - kstep)));
/* add the step */ /* add the step */
if ((err = mp_add_d(a, step, a)) != MP_OKAY) { if ((err = mp_add_d(a, step, a)) != MP_OKAY) {
@ -111,7 +107,7 @@ mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
} }
/* if didn't pass sieve and step == MP_MAX then skip test */ /* if didn't pass sieve and step == MP_MAX then skip test */
if ((y == 1) && (step >= (((mp_digit)1 << MP_DIGIT_BIT) - kstep))) { if (y && (step >= (((mp_digit)1 << MP_DIGIT_BIT) - kstep))) {
continue; continue;
} }
@ -123,7 +119,6 @@ mp_err mp_prime_next_prime(mp_int *a, int t, bool bbs_style)
} }
} }
err = MP_OKAY;
LBL_ERR: LBL_ERR:
mp_clear(&b); mp_clear(&b);
return err; return err;

2
mp_prime_strong_lucas_selfridge.c
View File

@ -192,7 +192,7 @@ mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, bool *result)
if ((err = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY) goto LBL_LS_ERR; if ((err = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY) goto LBL_LS_ERR;
if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) goto LBL_LS_ERR; if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) goto LBL_LS_ERR;
if (s_mp_get_bit(&Dz, (unsigned int)u)) { if (s_mp_get_bit(&Dz, u)) {
/* Formulas for addition of indices (carried out mod N); /* Formulas for addition of indices (carried out mod N);
* *
* U_(m+n) = (U_m*V_n + U_n*V_m)/2 * U_(m+n) = (U_m*V_n + U_n*V_m)/2

6
s_mp_get_bit.c
View File

@ -5,12 +5,12 @@
/* SPDX-License-Identifier: Unlicense */ /* SPDX-License-Identifier: Unlicense */
/* Get bit at position b and return true if the bit is 1, false if it is 0 */ /* Get bit at position b and return true if the bit is 1, false if it is 0 */
bool s_mp_get_bit(const mp_int *a, unsigned int b) bool s_mp_get_bit(const mp_int *a, int b)
{ {
mp_digit bit; mp_digit bit;
int limb = (int)(b / MP_DIGIT_BIT); int limb = b / MP_DIGIT_BIT;
if (limb >= a->used) { if (limb < 0 || limb >= a->used) {
return false; return false;
} }

18
s_mp_prime_is_divisible.c
View File

@ -10,16 +10,12 @@
*/ */
mp_err s_mp_prime_is_divisible(const mp_int *a, bool *result) mp_err s_mp_prime_is_divisible(const mp_int *a, bool *result)
{ {
int ix; int i;
mp_err err; for (i = 0; i < MP_PRIME_TAB_SIZE; i++) {
mp_digit res; /* what is a mod LBL_prime_tab[i] */
mp_err err;
/* default to not */ mp_digit res;
*result = false; if ((err = mp_mod_d(a, s_mp_prime_tab[i], &res)) != MP_OKAY) {
for (ix = 0; ix < MP_PRIME_TAB_SIZE; ix++) {
/* what is a mod LBL_prime_tab[ix] */
if ((err = mp_mod_d(a, s_mp_prime_tab[ix], &res)) != MP_OKAY) {
return err; return err;
} }
@ -30,6 +26,8 @@ mp_err s_mp_prime_is_divisible(const mp_int *a, bool *result)
} }
} }
/* default to not */
*result = false;
return MP_OKAY; return MP_OKAY;
} }
#endif #endif

2
tommath_private.h
View File

@ -188,7 +188,7 @@ MP_STATIC_ASSERT(prec_geq_min_prec, MP_PREC >= MP_MIN_PREC)
extern MP_PRIVATE mp_err(*s_mp_rand_source)(void *out, size_t size); extern MP_PRIVATE mp_err(*s_mp_rand_source)(void *out, size_t size);
/* lowlevel functions, do not call! */ /* lowlevel functions, do not call! */
MP_PRIVATE bool s_mp_get_bit(const mp_int *a, unsigned int b); MP_PRIVATE bool s_mp_get_bit(const mp_int *a, int b);
MP_PRIVATE mp_err s_mp_add(const mp_int *a, const mp_int *b, mp_int *c) MP_WUR; MP_PRIVATE mp_err s_mp_add(const mp_int *a, const mp_int *b, mp_int *c) MP_WUR;
MP_PRIVATE mp_err s_mp_sub(const mp_int *a, const mp_int *b, mp_int *c) MP_WUR; MP_PRIVATE mp_err s_mp_sub(const mp_int *a, const mp_int *b, mp_int *c) MP_WUR;
MP_PRIVATE mp_err s_mp_mul_digs_fast(const mp_int *a, const mp_int *b, mp_int *c, int digs) MP_WUR; MP_PRIVATE mp_err s_mp_mul_digs_fast(const mp_int *a, const mp_int *b, mp_int *c, int digs) MP_WUR;