diff options
Diffstat (limited to 'libtommath/bn_mp_prime_strong_lucas_selfridge.c')
| -rw-r--r-- | libtommath/bn_mp_prime_strong_lucas_selfridge.c | 258 |
1 files changed, 68 insertions, 190 deletions
diff --git a/libtommath/bn_mp_prime_strong_lucas_selfridge.c b/libtommath/bn_mp_prime_strong_lucas_selfridge.c index ae4fc8f..b50bbcd 100644 --- a/libtommath/bn_mp_prime_strong_lucas_selfridge.c +++ b/libtommath/bn_mp_prime_strong_lucas_selfridge.c @@ -1,22 +1,13 @@ #include "tommath_private.h" #ifdef BN_MP_PRIME_STRONG_LUCAS_SELFRIDGE_C -/* LibTomMath, multiple-precision integer library -- Tom St Denis - * - * LibTomMath is a library that provides multiple-precision - * integer arithmetic as well as number theoretic functionality. - * - * The library was designed directly after the MPI library by - * Michael Fromberger but has been written from scratch with - * additional optimizations in place. - * - * SPDX-License-Identifier: Unlicense - */ +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ /* * See file bn_mp_prime_is_prime.c or the documentation in doc/bn.tex for the details */ -#ifndef LTM_USE_FIPS_ONLY +#ifndef LTM_USE_ONLY_MR /* * 8-bit is just too small. You can try the Frobenius test @@ -28,33 +19,21 @@ * multiply bigint a with int d and put the result in c * Like mp_mul_d() but with a signed long as the small input */ -static int s_mp_mul_si(const mp_int *a, long d, mp_int *c) +static mp_err s_mp_mul_si(const mp_int *a, int32_t d, mp_int *c) { mp_int t; - int err, neg = 0; + mp_err err; if ((err = mp_init(&t)) != MP_OKAY) { return err; } - if (d < 0) { - neg = 1; - d = -d; - } /* * mp_digit might be smaller than a long, which excludes * the use of mp_mul_d() here. */ - if ((err = mp_set_long(&t, (unsigned long) d)) != MP_OKAY) { - goto LBL_MPMULSI_ERR; - } - if ((err = mp_mul(a, &t, c)) != MP_OKAY) { - goto LBL_MPMULSI_ERR; - } - if (neg == 1) { - c->sign = (a->sign == MP_NEG) ? MP_ZPOS: MP_NEG; - } -LBL_MPMULSI_ERR: + mp_set_i32(&t, d); + err = mp_mul(a, &t, c); mp_clear(&t); return err; } @@ -75,14 +54,14 @@ LBL_MPMULSI_ERR: (If that name sounds familiar, he is the guy who found the fdiv bug in the Pentium (P5x, I think) Intel processor) */ -int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) +mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, mp_bool *result) { /* CZ TODO: choose better variable names! */ mp_int Dz, gcd, Np1, Uz, Vz, U2mz, V2mz, Qmz, Q2mz, Qkdz, T1z, T2z, T3z, T4z, Q2kdz; /* CZ TODO: Some of them need the full 32 bit, hence the (temporary) exclusion of MP_8BIT */ int32_t D, Ds, J, sign, P, Q, r, s, u, Nbits; - int e; - int isset, oddness; + mp_err err; + mp_bool oddness; *result = MP_NO; /* @@ -93,9 +72,9 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) included. */ - if ((e = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz, - NULL)) != MP_OKAY) { - return e; + if ((err = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz, + NULL)) != MP_OKAY) { + return err; } D = 5; @@ -104,12 +83,9 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) for (;;) { Ds = sign * D; sign = -sign; - if ((e = mp_set_long(&Dz, (unsigned long)D)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_gcd(a, &Dz, &gcd)) != MP_OKAY) { - goto LBL_LS_ERR; - } + mp_set_u32(&Dz, (uint32_t)D); + if ((err = mp_gcd(a, &Dz, &gcd)) != MP_OKAY) goto LBL_LS_ERR; + /* if 1 < GCD < N then N is composite with factor "D", and Jacobi(D,N) is technically undefined (but often returned as zero). */ @@ -119,9 +95,7 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) if (Ds < 0) { Dz.sign = MP_NEG; } - if ((e = mp_kronecker(&Dz, a, &J)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_kronecker(&Dz, a, &J)) != MP_OKAY) goto LBL_LS_ERR; if (J == -1) { break; @@ -129,7 +103,7 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) D += 2; if (D > (INT_MAX - 2)) { - e = MP_VAL; + err = MP_VAL; goto LBL_LS_ERR; } } @@ -169,9 +143,7 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) Baillie-PSW test based on the strong Lucas-Selfridge test should be more reliable. */ - if ((e = mp_add_d(a, 1uL, &Np1)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_add_d(a, 1uL, &Np1)) != MP_OKAY) goto LBL_LS_ERR; s = mp_cnt_lsb(&Np1); /* CZ @@ -181,9 +153,7 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) * dividing an even number by two does not produce * any leftovers. */ - if ((e = mp_div_2d(&Np1, s, &Dz, NULL)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_div_2d(&Np1, s, &Dz, NULL)) != MP_OKAY) goto LBL_LS_ERR; /* We must now compute U_d and V_d. Since d is odd, the accumulated values U and V are initialized to U_1 and V_1 (if the target index were even, U and V would be initialized instead to U_0=0 @@ -200,34 +170,10 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) mp_set(&U2mz, 1uL); /* U_1 */ mp_set(&V2mz, (mp_digit)P); /* V_1 */ - if (Q < 0) { - Q = -Q; - if ((e = mp_set_long(&Qmz, (unsigned long)Q)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - /* Initializes calculation of Q^d */ - if ((e = mp_set_long(&Qkdz, (unsigned long)Q)) != MP_OKAY) { - goto LBL_LS_ERR; - } - Qmz.sign = MP_NEG; - Q2mz.sign = MP_NEG; - Qkdz.sign = MP_NEG; - Q = -Q; - } else { - if ((e = mp_set_long(&Qmz, (unsigned long)Q)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - /* Initializes calculation of Q^d */ - if ((e = mp_set_long(&Qkdz, (unsigned long)Q)) != MP_OKAY) { - goto LBL_LS_ERR; - } - } + mp_set_i32(&Qmz, Q); + if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) goto LBL_LS_ERR; + /* Initializes calculation of Q^d */ + mp_set_i32(&Qkdz, Q); Nbits = mp_count_bits(&Dz); @@ -240,37 +186,20 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) * V_2m = V_m*V_m - 2*Q^m */ - if ((e = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mod(&U2mz, a, &U2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_sqr(&V2mz, &V2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_sub(&V2mz, &Q2mz, &V2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mod(&V2mz, a, &V2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mod(&U2mz, a, &U2mz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_sqr(&V2mz, &V2mz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_sub(&V2mz, &Q2mz, &V2mz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mod(&V2mz, a, &V2mz)) != MP_OKAY) goto LBL_LS_ERR; + /* Must calculate powers of Q for use in V_2m, also for Q^d later */ - if ((e = mp_sqr(&Qmz, &Qmz)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_sqr(&Qmz, &Qmz)) != MP_OKAY) goto LBL_LS_ERR; + /* prevents overflow */ /* CZ still necessary without a fixed prealloc'd mem.? */ - if ((e = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((isset = mp_get_bit(&Dz, u)) == MP_VAL) { - e = isset; - goto LBL_LS_ERR; - } - if (isset == MP_YES) { + 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 (s_mp_get_bit(&Dz, (unsigned int)u) == MP_YES) { /* Formulas for addition of indices (carried out mod N); * * U_(m+n) = (U_m*V_n + U_n*V_m)/2 @@ -278,79 +207,46 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) * * Be careful with division by 2 (mod N)! */ - if ((e = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mul(&Uz, &V2mz, &T2z)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mul(&V2mz, &Vz, &T3z)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mul(&U2mz, &Uz, &T4z)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = s_mp_mul_si(&T4z, (long)Ds, &T4z)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_add(&T1z, &T2z, &Uz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if (mp_isodd(&Uz) != MP_NO) { - if ((e = mp_add(&Uz, a, &Uz)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mul(&Uz, &V2mz, &T2z)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mul(&V2mz, &Vz, &T3z)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mul(&U2mz, &Uz, &T4z)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = s_mp_mul_si(&T4z, Ds, &T4z)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_add(&T1z, &T2z, &Uz)) != MP_OKAY) goto LBL_LS_ERR; + if (MP_IS_ODD(&Uz)) { + if ((err = mp_add(&Uz, a, &Uz)) != MP_OKAY) goto LBL_LS_ERR; } /* CZ * This should round towards negative infinity because * Thomas R. Nicely used GMP's mpz_fdiv_q_2exp(). * But mp_div_2() does not do so, it is truncating instead. */ - oddness = mp_isodd(&Uz); - if ((e = mp_div_2(&Uz, &Uz)) != MP_OKAY) { - goto LBL_LS_ERR; - } + oddness = MP_IS_ODD(&Uz) ? MP_YES : MP_NO; + if ((err = mp_div_2(&Uz, &Uz)) != MP_OKAY) goto LBL_LS_ERR; if ((Uz.sign == MP_NEG) && (oddness != MP_NO)) { - if ((e = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - } - if ((e = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; + if ((err = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY) goto LBL_LS_ERR; } - if (mp_isodd(&Vz) != MP_NO) { - if ((e = mp_add(&Vz, a, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - } - oddness = mp_isodd(&Vz); - if ((e = mp_div_2(&Vz, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; + if ((err = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY) goto LBL_LS_ERR; + if (MP_IS_ODD(&Vz)) { + if ((err = mp_add(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR; } + oddness = MP_IS_ODD(&Vz) ? MP_YES : MP_NO; + if ((err = mp_div_2(&Vz, &Vz)) != MP_OKAY) goto LBL_LS_ERR; if ((Vz.sign == MP_NEG) && (oddness != MP_NO)) { - if ((e = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - } - if ((e = mp_mod(&Uz, a, &Uz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mod(&Vz, a, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; + if ((err = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY) goto LBL_LS_ERR; } + if ((err = mp_mod(&Uz, a, &Uz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR; + /* Calculating Q^d for later use */ - if ((e = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR; } } /* If U_d or V_d is congruent to 0 mod N, then N is a prime or a strong Lucas pseudoprime. */ - if ((mp_iszero(&Uz) != MP_NO) || (mp_iszero(&Vz) != MP_NO)) { + if (MP_IS_ZERO(&Uz) || MP_IS_ZERO(&Vz)) { *result = MP_YES; goto LBL_LS_ERR; } @@ -367,45 +263,27 @@ int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result) Lucas pseudoprime. */ /* Initialize 2*Q^(d*2^r) for V_2m */ - if ((e = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) goto LBL_LS_ERR; for (r = 1; r < s; r++) { - if ((e = mp_sqr(&Vz, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_sub(&Vz, &Q2kdz, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mod(&Vz, a, &Vz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if (mp_iszero(&Vz) != MP_NO) { + if ((err = mp_sqr(&Vz, &Vz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_sub(&Vz, &Q2kdz, &Vz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR; + if (MP_IS_ZERO(&Vz)) { *result = MP_YES; goto LBL_LS_ERR; } /* Calculate Q^{d*2^r} for next r (final iteration irrelevant). */ if (r < (s - 1)) { - if ((e = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) { - goto LBL_LS_ERR; - } - if ((e = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) { - goto LBL_LS_ERR; - } + if ((err = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR; + if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) goto LBL_LS_ERR; } } LBL_LS_ERR: mp_clear_multi(&Q2kdz, &T4z, &T3z, &T2z, &T1z, &Qkdz, &Q2mz, &Qmz, &V2mz, &U2mz, &Vz, &Uz, &Np1, &gcd, &Dz, NULL); - return e; + return err; } #endif #endif #endif - -/* ref: HEAD -> master, tag: v1.1.0 */ -/* git commit: 08549ad6bc8b0cede0b357a9c341c5c6473a9c55 */ -/* commit time: 2019-01-28 20:32:32 +0100 */ |