dropbear: libtommath/bn_mp_prime_strong_lucas

comparison libtommath/bn_mp_prime_strong_lucas_selfridge.c @ 1692:1051e4eea25a

Update LibTomMath to 1.2.0 (#84) * update C files * update other files * update headers * update makefiles * remove mp_set/get_double() * use ltm 1.2.0 API * update ltm_desc * use bundled tommath if system-tommath is too old * XMALLOC etc. were changed to MP_MALLOC etc.

author	Steffen Jaeckel <s@jaeckel.eu>
date	Tue, 26 May 2020 17:36:47 +0200
parents	f52919ffd3b1
children

comparison

equal deleted inserted replaced

-:2d3745d58843
+:1051e4eea25a
 #include "tommath_private.h"
 #ifdef BN_MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
-/* LibTomMath, multiple-precision integer library -- Tom St Denis
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
-*
+/* SPDX-License-Identifier: Unlicense */
-* 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
-*/
 /*
 *  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
 *  but that frobenius test can fail, too, for the same reason.
 */
 /*
 * 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) {
+mp_set_i32(&t, d);
-goto LBL_MPMULSI_ERR;
+err = mp_mul(a, &t, c);
-}
-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_clear(&t);
 return err;
 }
 /*
 Strong Lucas-Selfridge test.
 Additional comments marked "CZ" (without the quotes) are by the code-portist.
 (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;
+mp_err err;
-int isset, oddness;
+mp_bool oddness;
 *result = MP_NO;
 /*
 Find the first element D in the sequence {5, -7, 9, -11, 13, ...}
 such that Jacobi(D,N) = -1 (Selfridge's algorithm). Theory
 indicates that, if N is not a perfect square, D will "nearly
 always" be "small." Just in case, an overflow trap for D is
 included.
 */
-if ((e = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz,
+if ((err = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz,
 NULL)) != MP_OKAY) {
-return e;
+return err;
 }
 D = 5;
 sign = 1;
 for (;;) {
 Ds   = sign * D;
 sign = -sign;
-if ((e = mp_set_long(&Dz, (unsigned long)D)) != MP_OKAY) {
+mp_set_u32(&Dz, (uint32_t)D);
-goto LBL_LS_ERR;
+if ((err = mp_gcd(a, &Dz, &gcd)) != MP_OKAY)                goto LBL_LS_ERR;
-}
-if ((e = 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). */
 if ((mp_cmp_d(&gcd, 1uL) == MP_GT) && (mp_cmp(&gcd, a) == MP_LT)) {
 goto LBL_LS_ERR;
 }
 if (Ds < 0) {
 Dz.sign = MP_NEG;
 }
-if ((e = mp_kronecker(&Dz, a, &J)) != MP_OKAY) {
+if ((err = mp_kronecker(&Dz, a, &J)) != MP_OKAY)            goto LBL_LS_ERR;
-goto LBL_LS_ERR;
-}
 if (J == -1) {
 break;
 }
 D += 2;
 if (D > (INT_MAX - 2)) {
-e = MP_VAL;
+err = MP_VAL;
 goto LBL_LS_ERR;
 }
 }
 the evidence indicates that the strong Lucas-Selfridge test is
 more effective than the standard Lucas-Selfridge test, and a
 Baillie-PSW test based on the strong Lucas-Selfridge test
 should be more reliable. */
-if ((e = mp_add_d(a, 1uL, &Np1)) != MP_OKAY) {
+if ((err = mp_add_d(a, 1uL, &Np1)) != MP_OKAY)                 goto LBL_LS_ERR;
-goto LBL_LS_ERR;
-}
 s = mp_cnt_lsb(&Np1);
 /* CZ
 * This should round towards zero because
 * Thomas R. Nicely used GMP's mpz_tdiv_q_2exp()
 * and mp_div_2d() is equivalent. Additionally:
 * dividing an even number by two does not produce
 * any leftovers.
 */
-if ((e = mp_div_2d(&Np1, s, &Dz, NULL)) != MP_OKAY) {
+if ((err = mp_div_2d(&Np1, s, &Dz, NULL)) != MP_OKAY)          goto LBL_LS_ERR;
-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
 and V_0=2). The values of U_2m and V_2m are also initialized to
 U_1 and V_1; the FOR loop calculates in succession U_2 and V_2,
 mp_set(&Uz, 1uL);    /* U=U_1 */
 mp_set(&Vz, (mp_digit)P);    /* V=V_1 */
 mp_set(&U2mz, 1uL);  /* U_1 */
 mp_set(&V2mz, (mp_digit)P);  /* V_1 */
-if (Q < 0) {
+mp_set_i32(&Qmz, Q);
-Q = -Q;
+if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY)                  goto LBL_LS_ERR;
-if ((e = mp_set_long(&Qmz, (unsigned long)Q)) != MP_OKAY) {
+/* Initializes calculation of Q^d */
-goto LBL_LS_ERR;
+mp_set_i32(&Qkdz, Q);
-}
-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;
-}
-}
 Nbits = mp_count_bits(&Dz);
 for (u = 1; u < Nbits; u++) { /* zero bit off, already accounted for */
 /* Formulas for doubling of indices (carried out mod N). Note that
 *
 * U_2m = U_m*V_m
 * V_2m = V_m*V_m - 2*Q^m
 */
-if ((e = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY) {
+if ((err = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY)         goto LBL_LS_ERR;
-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 ((e = mp_mod(&U2mz, a, &U2mz)) != MP_OKAY) {
+if ((err = mp_sub(&V2mz, &Q2mz, &V2mz)) != MP_OKAY)         goto LBL_LS_ERR;
-goto LBL_LS_ERR;
+if ((err = mp_mod(&V2mz, a, &V2mz)) != 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;
-}
 /* Must calculate powers of Q for use in V_2m, also for Q^d later */
-if ((e = mp_sqr(&Qmz, &Qmz)) != MP_OKAY) {
+if ((err = mp_sqr(&Qmz, &Qmz)) != MP_OKAY)                  goto LBL_LS_ERR;
-goto LBL_LS_ERR;
-}
 /* prevents overflow */ /* CZ  still necessary without a fixed prealloc'd mem.? */
-if ((e = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY) {
+if ((err = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY)               goto LBL_LS_ERR;
-goto LBL_LS_ERR;
+if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY)               goto LBL_LS_ERR;
-}
-if ((e = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) {
+if (s_mp_get_bit(&Dz, (unsigned int)u) == MP_YES) {
-goto LBL_LS_ERR;
-}
-if ((isset = mp_get_bit(&Dz, u)) == MP_VAL) {
-e = isset;
-goto LBL_LS_ERR;
-}
-if (isset == MP_YES) {
 /* Formulas for addition of indices (carried out mod N);
 *
 * U_(m+n) = (U_m*V_n + U_n*V_m)/2
 * V_(m+n) = (V_m*V_n + D*U_m*U_n)/2
 *
 * Be careful with division by 2 (mod N)!
 */
-if ((e = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY) {
+if ((err = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY)         goto LBL_LS_ERR;
-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 ((e = mp_mul(&Uz, &V2mz, &T2z)) != MP_OKAY) {
+if ((err = mp_mul(&U2mz, &Uz, &T4z)) != MP_OKAY)         goto LBL_LS_ERR;
-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 ((e = mp_mul(&V2mz, &Vz, &T3z)) != MP_OKAY) {
+if (MP_IS_ODD(&Uz)) {
-goto LBL_LS_ERR;
+if ((err = mp_add(&Uz, a, &Uz)) != 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;
-}
 }
 /* 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);
+oddness = MP_IS_ODD(&Uz) ? MP_YES : MP_NO;
-if ((e = mp_div_2(&Uz, &Uz)) != MP_OKAY) {
+if ((err = mp_div_2(&Uz, &Uz)) != MP_OKAY)               goto LBL_LS_ERR;
-goto LBL_LS_ERR;
-}
 if ((Uz.sign == MP_NEG) && (oddness != MP_NO)) {
-if ((e = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY) {
+if ((err = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY)       goto LBL_LS_ERR;
-goto LBL_LS_ERR;
+}
-}
+if ((err = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY)          goto LBL_LS_ERR;
-}
+if (MP_IS_ODD(&Vz)) {
-if ((e = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY) {
+if ((err = mp_add(&Vz, a, &Vz)) != MP_OKAY)           goto LBL_LS_ERR;
-goto LBL_LS_ERR;
+}
-}
+oddness = MP_IS_ODD(&Vz) ? MP_YES : MP_NO;
-if (mp_isodd(&Vz) != MP_NO) {
+if ((err = mp_div_2(&Vz, &Vz)) != MP_OKAY)               goto LBL_LS_ERR;
-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 ((Vz.sign == MP_NEG) && (oddness != MP_NO)) {
-if ((e = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY) {
+if ((err = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY)       goto LBL_LS_ERR;
-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;
-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;
-}
 /* Calculating Q^d for later use */
-if ((e = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY) {
+if ((err = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY)       goto LBL_LS_ERR;
-goto LBL_LS_ERR;
+if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY)          goto LBL_LS_ERR;
-}
-if ((e = 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;
 }
 /* NOTE: Ribenboim ("The new book of prime number records," 3rd ed.,
 by repeated use of the formula V_2m = V_m*V_m - 2*Q^m. If any of
 these are congruent to 0 mod N, then N is a prime or a strong
 Lucas pseudoprime. */
 /* Initialize 2*Q^(d*2^r) for V_2m */
-if ((e = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) {
+if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY)                goto LBL_LS_ERR;
-goto LBL_LS_ERR;
-}
 for (r = 1; r < s; r++) {
-if ((e = mp_sqr(&Vz, &Vz)) != MP_OKAY) {
+if ((err = mp_sqr(&Vz, &Vz)) != MP_OKAY)                    goto LBL_LS_ERR;
-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 ((e = mp_sub(&Vz, &Q2kdz, &Vz)) != MP_OKAY) {
+if (MP_IS_ZERO(&Vz)) {
-goto LBL_LS_ERR;
-}
-if ((e = mp_mod(&Vz, a, &Vz)) != MP_OKAY) {
-goto LBL_LS_ERR;
-}
-if (mp_iszero(&Vz) != MP_NO) {
 *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) {
+if ((err = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY)             goto LBL_LS_ERR;
-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;
-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;
-}
 }
 }
 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 */

Mercurial > dropbear

comparison libtommath/bn_mp_prime_strong_lucas_selfridge.c @ 1692:1051e4eea25a