CryptX
view release on metacpan or search on metacpan
src/ltm/mp_prime_strong_lucas_selfridge.c view on Meta::CPAN
err = MP_VAL;
goto LBL_LS_ERR;
}
}
P = 1; /* Selfridge's choice */
Q = (1 - Ds) / 4; /* Required so D = P*P - 4*Q */
/* NOTE: The conditions (a) N does not divide Q, and
(b) D is square-free or not a perfect square, are included by
some authors; e.g., "Prime numbers and computer methods for
factorization," Hans Riesel (2nd ed., 1994, Birkhauser, Boston),
p. 130. For this particular application of Lucas sequences,
these conditions were found to be immaterial. */
/* Now calculate N - Jacobi(D,N) = N + 1 (even), and calculate the
odd positive integer d and positive integer s for which
N + 1 = 2^s*d (similar to the step for N - 1 in Miller's test).
The strong Lucas-Selfridge test then returns N as a strong
Lucas probable prime (slprp) if any of the following
conditions is met: U_d=0, V_d=0, V_2d=0, V_4d=0, V_8d=0,
V_16d=0, ..., etc., ending with V_{2^(s-1)*d}=V_{(N+1)/2}=0
(all equalities mod N). Thus d is the highest index of U that
must be computed (since V_2m is independent of U), compared
to U_{N+1} for the standard Lucas-Selfridge test; and no
index of V beyond (N+1)/2 is required, just as in the
standard Lucas-Selfridge test. However, the quantity Q^d must
be computed for use (if necessary) in the latter stages of
the test. The result is that the strong Lucas-Selfridge test
has a running time only slightly greater (order of 10 %) than
that of the standard Lucas-Selfridge test, while producing
only (roughly) 30 % as many pseudoprimes (and every strong
Lucas pseudoprime is also a standard Lucas pseudoprime). Thus
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 ((err = mp_add_d(a, 1uL, &Np1)) != MP_OKAY) 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 ((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
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,
U_4 and V_4, U_8 and V_8, etc. If the corresponding bits
(1, 2, 3, ...) of t are on (the zero bit having been accounted
for in the initialization of U and V), these values are then
combined with the previous totals for U and V, using the
composition formulas for addition of indices. */
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 */
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);
for (u = 1; u < Nbits; u++) { /* zero bit off, already accounted for */
/* Formulas for doubling of indices (carried out mod N). Note that
* the indices denoted as "2m" are actually powers of 2, specifically
* 2^(ul-1) beginning each loop and 2^ul ending each loop.
*
* U_2m = U_m*V_m
* V_2m = V_m*V_m - 2*Q^m
*/
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 ((err = mp_sqr(&Qmz, &Qmz)) != MP_OKAY) goto LBL_LS_ERR;
/* prevents overflow */ /* CZ still necessary without a fixed prealloc'd mem.? */
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, u)) {
/* 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 ((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_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_isodd(&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 ((err = mp_div_2(&Uz, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
if (mp_isneg(&Uz) && oddness) {
if ((err = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
}
if ((err = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
if (mp_isodd(&Vz)) {
if ((err = mp_add(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
}
oddness = mp_isodd(&Vz);
if ((err = mp_div_2(&Vz, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
if (mp_isneg(&Vz) && oddness) {
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 ((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_iszero(&Vz)) {
*result = MP_YES;
goto LBL_LS_ERR;
}
/* NOTE: Ribenboim ("The new book of prime number records," 3rd ed.,
1995/6) omits the condition V0 on p.142, but includes it on
p. 130. The condition is NECESSARY; otherwise the test will
return MP_NO negatives---e.g., the primes 29 and 2000029 will be
returned as composite. */
/* Otherwise, we must compute V_2d, V_4d, V_8d, ..., V_{2^(s-1)*d}
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 ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) goto LBL_LS_ERR;
for (r = 1; r < s; r++) {
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_iszero(&Vz)) {
*result = MP_YES;
goto LBL_LS_ERR;
}
/* Calculate Q^{d*2^r} for next r (final iteration irrelevant). */
if (r < (s - 1)) {
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 err;
}
#endif
#endif
( run in 0.941 second using v1.01-cache-2.11-cpan-364913b4093 )