* Copyright (c) 1993 Landon Curt Noll
* Permission is granted to use, distribute, or modify this source,
* provided that this copyright notice remains intact.
* chongo@toad.com -or- ...!{pyramid,sun,uunet}!hoptoad!chongo
* lucas - perform a Lucas primality test on h*2^n-1
* On 6 August 1989 at 00:53 PDT, the 'Amdahl 6', a team consisting of
* John Brown, Landon Curt Noll, Bodo Parady, Gene Smith, Joel Smith and
* Sergio Zarantonello proved the following 65087 digit number to be prime:
* At the time of discovery, this number was the largest known prime.
* The primality was demonstrated by a program implementing the test
* found in these routines. An Amdahl 1200 takes 1987 seconds to test
* the primality of this number. A Cray 2 took several hours to
* confirm this prime. As of 28 Aug 1993, this prime was the 2nd
* largest known prime and the largest known non-Mersenne prime.
* The same team also discovered the following twin prime pair:
* 1706595 * 2 -1 1706595 * 2 +1
* As of 28 Aug 1993, these primes was still the largest known twin prime pair.
* ON GAINING A WORLD RECORD:
* The routines in calc were designed to be portable, and to work on
* numbers of 'sane' size. The Amdahl 6 team used a 'ultra-high speed
* multi-precision' package that a machine dependent collection of routines
* tuned for a long trace vector processor to work with very large numbers.
* The heart of the package was a multiplication and square routine that
* was based on the PFA Fast Fourier Transform and on Winograd's radix FFTs.
* Having a fast computer, and a good multi-precision package are
* critical, but one also needs to know where to look in order to have
* a good chance at a record. Knowing what to test is beyond the scope
* of this routine. However the following observations are noted:
* test numbers of the form h*2^n-1
* fix a value of n and vary the value h
* h*2^n-1 is not divisible by any small prime < 2^40
* h*2^n+1 is not divisible by any small prime < 2^40
* The Mersenne test for '2^n-1' is the fastest known primality test
* for a given large numbers. However, it is faster to search for
* primes of the form 'h*2^n-1'. When n is around 20000, one can find
* a prime of the form 'h*2^n-1' in about 1/2 the time.
* Critical to understanding why 'h*2^n-1' is to observe that primes of
* the form '2^n-1' seem to bunch around "islands". Such "islands"
* seem to be getting fewer and farther in-between, forcing the time
* for each test to grow longer and longer (worse then O(n^2 log n)).
* On the other hand, when one tests 'h*2^n-1', fixes 'n' and varies
* 'h', the time to test each number remains relatively constant.
* It is clearly a win to eliminate potential test candidates by
* rejecting numbers that that are divisible by 'small' primes. We
* (the "Amdahl 6") rejected all numbers that were divisible by primes
* less than '2^40'. We stopped looking for small factors at '2^40'
* when the rate of candidates being eliminated was slowed down to
* The 'n mod 128 == 0' restriction allows one to test for divisibility
* of small primes more quickly. To test of 'q' is a factor of 'k*2^n-1',
* one check to see if 'k*2^n mod q' == 1, which is the same a checking
* if 'h*(2^n mod q) mod q' == 1. One can compute '2^n mod q' by making
* 2^(2x) mod q == y^2 mod q 0 bit
* 2^(2x+1) mod q == 2*y^2 mod q 1 bit
* The choice of which expression depends on the binary pattern of 'n'.
* Since '1' bits require an extra step (multiply by 2), one should
* select value of 'n' that contain mostly '0' bits. The restriction
* of 'n mod 128 == 0' ensures that the bottom 7 bits of 'n' are 0.
* By limiting 'h' to '2^39' and eliminating all values divisible by
* small primes < twice the 'h' limit (2^40), one knows that all
* remaining candidates are relatively prime. Thus, when a candidate
* is proven to be composite (not prime) by the big test, one knows
* that the factors for that number (whatever they may be) will not
* be the factors of another candidate.
* Finally, one should eliminate all values of 'h*2^n-1' where
* 'h*2^n+1' is divisible by a small primes. The ideas behind this
* point is beyond the scope of this program.
global pprod256; /* product of "primes up to 256" / "primes up to 46" */
global lib_debug; /* 1 => print debug statements */
* lucas - lucas primality test on h*2^n-1
* This routine will perform a primality test on h*2^n-1 based on
* the mathematics of Lucas, Lehmer and Riesel. One should read
* "Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel,
* Mathematics of Computation, Vol 23 #108, pp. 869-875, Oct 1969
* The following book is also useful:
* "Prime numbers and Computer Methods for Factorization", by Hans Riesel,
* Birkhauser, 1985, pp 131-134, 278-285, 438-444
* A few useful Legendre identities may be found in:
* "Introduction to Analytic Number Theory", by Tom A. Apostol,
* Springer-Verlag, 1984, p 188.
* This test is performed as follows: (see Ref1, Theorem 5)
* a) generate u(0) (see the function gen_u0() below)
* b) generate u(n-2) according to the rule:
* u(i+1) = u(i)^2-2 mod h*2^n-1
* c) h*2^n-1 is prime if and only if u(n-2) == 0 Q.E.D. :-)
* Now the following conditions must be true for the test to work:
* In order to reduce the number of tests, as attempt to eliminate
* any number that is divisible by a prime less than 257. Valid prime
* candidates less than 257 are declared prime as a special case.
* The condition 'h mod 2 == 1' is not a problem. Say one is testing
* 'j*2^m-1', where j is even. If we note that:
* j mod 2^x == 0 for x>0 implies j*2^m-1 == ((j/2^x)*2^(m+x))-1,
* then we can let h=j/2^x and n=m+x and test 'h*2^n-1' which is the value.
* We need only consider odd values of h because we can rewrite our numbers
* 0 => h*2^n-1 is not prime
* -1 => a test could not be formed, or h >= 2^n, h <= 0, n <= 0
local testval; /* h*2^n-1 */
local shiftdown; /* the power of 2 that divides h */
local u; /* the u(i) sequence value */
local v1; /* the v(1) generator of u(0) */
local i; /* u sequence cycle number */
local oldh; /* pre-reduced h */
local oldn; /* pre-reduced n */
local bits; /* highbit of h*2^n-1 */
ldebug("lucas", "h is non-int");
quit "FATAL: bad args: h must be an integer";
ldebug("lucas", "n is non-int");
quit "FATAL: bad args: n must be an integer";
* we will force h to be odd by moving powers of two over to 2^n
shiftdown = fcnt(h,2); /* h % 2^shiftdown == 0, max shiftdown */
* enforce the 0 < h < 2^n rule
print "ERROR: reduced args violate the rule: 0 < h < 2^n";
print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
ldebug("lucas", "unknown: h <= 0 || n <= 0");
print "ERROR: reduced args violate the rule: h < 2^n";
print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
ldebug("lucas", "unknown: highbit(h) >= n");
* catch the degenerate case of h*2^n-1 == 1
ldebug("lucas", "not prime: h == 1 && n == 1");
return 0; /* 1*2^1-1 == 1 is not prime */
* catch the degenerate case of n==2
* n==2 and 0<h<2^n ==> 0<h<4
* Since h is now odd ==> h==1 or h==3
ldebug("lucas", "prime: h == 1 && n == 2");
return 1; /* 1*2^2-1 == 3 is prime */
ldebug("lucas", "prime: h == 3 && n == 2");
return 1; /* 3*2^2-1 == 11 is prime */
* catch small primes < 257
* We check for only a few primes because the other primes < 257
* violate the checks above.
if (n == 3 || n == 5 || n == 7) {
ldebug("lucas", "prime: 3, 7, 31, 127 are prime");
return 1; /* 3, 7, 31, 127 are prime */
if (n == 2 || n == 3 || n == 4 || n == 6) {
ldebug("lucas", "prime: 11, 23, 47, 191 are prime");
return 1; /* 11, 23, 47, 191 are prime */
ldebug("lucas", "prime: 79 is prime");
return 1; /* 79 is prime */
ldebug("lucas", "prime: 223 is prime");
return 1; /* 223 is prime */
ldebug("lucas", "prime: 239 is prime");
return 1; /* 239 is prime */
* Avoid any numbers divisible by small primes
* check for 3 <= prime factors < 29
* pfact(28)/2 = 111546435
if (gcd(testval, 111546435) > 1) {
/* a small 3 <= prime < 29 divides h*2^n-1 */
ldebug("lucas","not-prime: 3<=prime<29 divides h*2^n-1");
* check for 29 <= prime factors < 47
* pfact(46)/pfact(28) = 5864229
if (gcd(testval, 58642669) > 1) {
/* a small 29 <= prime < 47 divides h*2^n-1 */
ldebug("lucas","not-prime: 29<=prime<47 divides h*2^n-1");
* check for prime 47 <= factors < 257, if h*2^n-1 is large
* 2^282 > pfact(256)/pfact(46) > 2^281
pprod256 = pfact(256)/pfact(46);
if (gcd(testval, pprod256) > 1) {
/* a small 47 <= prime < 257 divides h*2^n-1 */
"not-prime: 47<=prime<257 divides h*2^n-1");
* We will use gen_v1() to give us a v(1) using the values
* of 'h' and 'n'. We will then use gen_u0() to convert
* If gen_v1() returns a negative value, then we failed to
* generate a test for h*2^n-1. This is because h mod 3 == 0
* is hard to do, and in rare cases, exceed the tables found
* in this program. We will generate an message and assume
* the number is not prime, even though if we had a larger
* table, we might have been able to show that it is prime.
v1 = gen_v1(h, n, testval);
/* failure to test number */
print "unable to compute v(1) for", h : "*2^" : n : "-1";
ldebug("lucas", "unknown: no v(1)");
u = gen_u0(h, n, testval, v1);
* return 1 if prime, 0 is not prime
ldebug("lucas", "prime: end of test");
ldebug("lucas", "not-prime: end of test");
* gen_u0 - determine the initial Lucas sequence for h*2^n-1
* According to Ref1, Theorem 5:
* u(0) = alpha^h + alpha^(-h)
* v(x) = alpha^x + alpha^(-x) (Ref1, bottom of page 872)
* We calculate v(h) as follows: (Ref1, top of page 873)
* v(0) = alpha^0 + alpha^(-0) = 2
* v(1) = alpha^1 + alpha^(-1) = gen_v1(h,n)
* v(n+2) = v(1)*v(n+1) - v(n)
* This function does not concern itself with the value of 'alpha'.
* The gen_v1() function is used to compute v(1), and identity
* functions take it from there.
* It can be shown that the following are true:
* v(2*n+1) = v(n+1)*v(n) - v(1)
* To prevent v(x) from growing too large, one may replace v(x) with
* `v(x) mod h*2^n-1' at any time.
* See the function gen_v1() for details on the value of v(1).
* h - h as in h*2^n-1 (h mod 2 != 0)
* v1 - gen_v1(h,n) (see function below)
* u(0) - initial value for Lucas test on h*2^n-1
* -1 - failed to generate u(0)
gen_u0(h, n, testval, v1)
local shiftdown; /* the power of 2 that divides h */
local r; /* low value: v(n) */
local s; /* high value: v(n+1) */
local hbits; /* highest bit set in h */
quit "bad args: h must be an integer";
quit "bad args: n must be an integer";
quit "bad args: testval must be an integer";
quit "bad args: v1 must be an integer";
quit "bogus arg: testval is <= 0";
quit "bogus arg: v1 is <= 0";
* enforce the h mod rules
quit "h must not be even";
* enforce the h > 0 and n >= 2 rules
quit "reduced args violate the rule: 0 < h < 2^n";
quit "reduced args violate the rule: 0 < h < 2^n";
* build up u2 based on the reversed bits of h
* deal with small h as a special case
* The h value is odd > 0, and it needs to be
* at least 2 bits long for the loop below to work.
ldebug("gen_u0", "quick h == 1 case");
/* cycle from second highest bit to second lowest bit of h */
for (i=hbits-1; i > 0; --i) {
/* compute v(2n+1) = v(r+1)*v(r)-v1 */
r = (r*s - v1) % testval;
/* compute v(2n+2) = v(r+1)^2-2 */
/* compute v(2n+1) = v(r+1)*v(r)-v1 */
s = (r*s - v1) % testval;
/* compute v(2n) = v(r)^-2 */
/* we know that h is odd, so the final bit(0) is 1 */
r = (r*s - v1) % testval;
/* compute the final u2 return value */
* Trial tables used by gen_v1()
* When h mod 3 == 0, one needs particular values of D, a and b (see gen_v1
* documentation) in order to find a value of v(1).
* This table defines 'quickmax' possible tests to be taken in ascending
* order. The v1_qval[x] refers to a v(1) value from Ref1, Table 1. A
* related D value is found in d_qval[x]. All D values expect d_qval[1]
* are also taken from Ref1, Table 1. The case of D == 21 as listed in
* Ref1, Table 1 can be changed to D == 7 for the sake of the test because
* It should be noted that the D values all satisfy the selection values
* as outlined in the gen_v1() function comments. That is:
* where f == 0 and g == 0, P == D. So we simply need to check that
* one of the following two cases are true:
* P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
* P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
* In all cases, the value of r is:
* r == Q*(2^j)*(3^k)*(z^2)
* where Q == 1. No further processing is needed to compute v(1) when r
d_qval[0] = 5; v1_qval[0] = 3; /* a=1 b=1 r=4 */
d_qval[1] = 7; v1_qval[1] = 5; /* a=3 b=1 r=12 D=21 */
d_qval[2] = 13; v1_qval[2] = 11; /* a=3 b=1 r=4 */
d_qval[3] = 11; v1_qval[3] = 20; /* a=3 b=1 r=2 */
d_qval[4] = 29; v1_qval[4] = 27; /* a=5 b=1 r=4 */
d_qval[5] = 53; v1_qval[5] = 51; /* a=53 b=1 r=4 */
d_qval[6] = 17; v1_qval[6] = 66; /* a=17 b=1 r=1 */
d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
* gen_v1 - compute the v(1) for a given h*2^n-1 if we can
* n > 2 (n==2 has already been eliminated)
* h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3)
* The generation of v(1) depends on the value of h. There are two cases
* to consider, h mod 3 != 0, and h mod 3 == 0.
* This case is easy and always finds v(1).
* In Ref1, page 869, one finds that if: (or see Ref2, page 131-132)
* which translates, gives the functions assumptions, into the condition:
* If this case condition is true, then:
* u(0) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
* = (2+sqrt(3))^h + (2+sqrt(3))^(-h)
* and since Ref1, Theorem 5 states:
* u(0) = alpha^h + alpha^(-h)
* r = abs(2^2 - 1^2*3) = 1
* and the bottom of Ref1, page 872 states:
* v(x) = alpha^x + alpha^(-x)
* v(1) = alpha^1 + alpha^(-1)
* = (2+sqrt(3)) + (2-sqrt(3))
* This case is not so easy and finds v(1) in most all cases. In this
* version of this program, we will simply return -1 (failure) if we
* hit one of the cases that fall thru the cracks. This does not happen
* often, so this is not too bad.
* Ref1, Theorem 5 contains the following definitions:
* alpha = (a + b*sqrt(D))^2/r
* where D is 'square free', and 'alpha = epsilon^s' (for some s>0) are units
* in the quadratic field K(sqrt(D)).
* One can find possible values for a, b and D in Ref1, Table 1 (page 872).
* (see the file lucas_tbl.cal)
* Now Ref1, Theorem 5 states that if:
* L(D, h*2^n-1) = -1 [condition 1]
* L(r, h*2^n-1) * (a^2 - b^2*D)/r = -1 [condition 2]
* where L(x,y) is the Legendre symbol (see below), then:
* u(0) = alpha^h + alpha^(-h)
* The bottom of Ref1, page 872 states:
* v(x) = alpha^x + alpha^(-x)
* v(1) = alpha^1 + alpha^(-1)
* Therefore we need to take a given (D,a,b), determine if the two conditions
* are true, and return the related v(1).
* Before we address the two conditions, we need some background information
* on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find
* the following definitions of J(a,p) and L(a,n):
* The Legendre symbol L(a,p) takes the value:
* L(a,p) == 1 => a is a quadratic residue of p
* L(a,p) == -1 => a is NOT a quadratic residue of p
* The value x is a quadratic residue of y if there exists some integer z
* The Jacobi symbol J(x,y) takes the value:
* J(x,y) == 1 => y is not prime, or x is a quadratic residue of y
* J(x,y) == -1 => x is NOT a quadratic residue of y
* In the following comments on Legendre and Jacobi identities, we shall
* assume that the arguments to the symbolic are valid over the symbol
* definitions as stated above.
* In Ref2, pp 280-284, we find that:
* L(a,p)*L(b,p) == L(a*b,p) {A3.5}
* J(x,y)*J(z,y) == J(x*z,y) {A3.14}
* L(a,p) == L(p,a) * (-1)^((a-1)*(p-1)/4) {A3.8}
* J(x,y) == J(y,x) * (-1)^((x-1)*(y-1)/4) {A3.17}
* The equality L(a,p) == J(a,p) when: {note 0}
* It can be shown that (see Ref3):
* L(a,p) == L(a mod p, p) {note 1}
* L(z^2, p) == 1 {note 2}
* p mod 8 == +/-1 implies L(2,p) == 1 {note 3}
* p mod 12 == +/-1 implies L(3,p) == 1 {note 4}
* Since h*2^n-1 mod 8 == -1, for n>2, note 3 implies:
* L(2, h*2^n-1) == 1 (n>2) {note 5}
* Since h=3*A, h*2^n-1 mod 12 == -1, for A>0, note 4 implies:
* L(3, h*2^n-1) == 1 {note 6}
* By use of {A3.5}, {note 2}, {note 5} and {note 6}, one can show:
* L((2^g)*(3^l)*(z^2), h*2^n-1) == 1 (g>=0,l>=0,z>0,n>2) {note 7}
* Returning to the testing of conditions, take condition 1:
* L(D, h*2^n-1) == -1 [condition 1]
* In order for J(D, h*2^n-1) to be defined, we must ensure that D
* is not a factor of h*2^n-1. This is done by pre-screening h*2^n-1 to
* not have small factors and selecting D less than that factor check limit.
* By use of {note 7}, we can show that when we choose D to be:
* D = P*(2^f)*(3^g) (P is prime>2)
* The square free condition implies f = 0 or 1, g = 0 or 1. If f and g
* are both 1, P must be a prime > 3.
* So given such a D value:
* L(D, h*2^n-1) == L(P*(2^g)*(3^l), h*2^n-1)
* == L(P, h*2^n-1) * L((2^g)*(3^l), h*2^n-1) {A3.5}
* == L(P, h*2^n-1) * 1 {note 7}
* == L(h*2^n-1, P)*(-1)^((h*2^n-2)*(P-1)/4) {A3.8}
* == L(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 1}
* == J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 0}
* When does J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) take the value of -1,
* thus satisfy [condition 1]? The answer depends on P. Now P is a prime>2,
* thus P mod 4 == 1 or -1.
* P mod 4 == 1 implies (-1)^((h*2^n-2)*(P-1)/4) == 1
* L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
* P mod 4 == -1 implies (-1)^((h*2^n-2)*(P-1)/4) == -1
* L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
* == L(h*2^n-1 mod P, P) * -1
* == -J(h*2^n-1 mod P, P)
* Therefore [condition 1] is met if, and only if, one of the following
* P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
* P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
* Now consider [condition 2]:
* L(r, h*2^n-1) * (a^2 - b^2*D)/r == -1 [condition 2]
* We select only a, b, r and D values where:
* Therefore in order for [condition 2] to be met, we must show that:
* If we select r to be of the form:
* r == Q*(2^j)*(3^k)*(z^2) (Q == 1, j>=0, k>=0, z>0)
* then by use of {note 7}:
* L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
* == L((2^j)*(3^k)*(z^2), h*2^n-1)
* and thus, [condition 2] is met.
* If we select r to be of the form:
* r == Q*(2^j)*(3^k)*(z^2) (Q is prime>2, j>=0, k>=0, z>0)
* then by use of {note 7}:
* L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
* == L(Q, h*2^n-1) * L((2^j)*(3^k)*(z^2), h*2^n-1) {A3.5}
* == L(Q, h*2^n-1) * 1 {note 2}
* == L(h*2^n-1, Q) * (-1)^((h*2^n-2)*(Q-1)/4) {A3.8}
* == L(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 1}
* == J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 0}
* When does J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) take the value of 1,
* thus satisfy [condition 2]? The answer depends on Q. Now Q is a prime>2,
* thus Q mod 4 == 1 or -1.
* Q mod 4 == 1 implies (-1)^((h*2^n-2)*(Q-1)/4) == 1
* L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
* Q mod 4 == -1 implies (-1)^((h*2^n-2)*(Q-1)/4) == -1
* L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
* == L(h*2^n-1 mod Q, Q) * -1
* == -J(h*2^n-1 mod Q, Q)
* Therefore [condition 2] is met by selecting D = Q*(2^j)*(3^k)*(z^2),
* where Q is prime>2, j>=0, k>=0, z>0; if and only if one of the following
* Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
* Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
* In conclusion, we can compute v(1) by attempting to do the following:
* alpha == (a + b*sqrt(D))^2/r
* v(1) = alpha^1 + alpha^(-1)
* if and only if we can find a given a, b, D that obey all the
* following selection rules:
* D == P*(2^f)*(3^g) (P is prime>2, f,g == 0 or 1)
* r == Q*(2^j)*(3^k)*(z^2) (Q==1 or Q is prime>2, j>=0, k>=0, z>0)
* one of the following is true:
* P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
* P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
* if Q is prime, then one of the following is true:
* Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
* Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
* If we cannot find a v(1) quickly enough, then we will give up
* testing h*2^n-1. This does not happen too often, so this hack
* returns v(1), or -1 is there is no quick way
local d; /* the 'D' value to try */
local val_mod; /* h*2^n-1 mod 'D' */
/* v(1) is easy to compute */
* We will try all 'D' values until we find a proper v(1)
* or run out of 'D' values.
for (i=0; i < quickmax; ++i) {
/* compute h*2^n-1 mod 'D' quickly */
val_mod = (h*pmod(2,n%(d-1),d)-1) % d;
* if 'D' mod 4 == 1, then
* (h*2^n-1) mod 'D' can not be a quadratic residue of 'D'
* (h*2^n-1) mod 'D' must be a quadratic residue of 'D'
/* D mod 4 == 1, so check for J(D, h*2^n-1) == -1 */
if (jacobi(val_mod, d) == -1) {
/* it worked, return the related v(1) value */
/* D mod 4 == -1, so check for J(D, h*2^n-1) == 1 */
if (jacobi(val_mod, d) == 1) {
/* it worked, return the related v(1) value */
* This is an example of a more complex proof construction.
* The code above will not be able to find the v(1) for:
* v(1)=81 D=6557 a=79 b=1 r=316
* Now, D==79*83 and r=79*2^2. If we show that:
* J(h*2^n-1 mod 79, 79) == -1
* J(h*2^n-1 mod 83, 83) == 1
* then we will satisfy [condition 1]. Observe:
* 79 mod 4 == -1 implies (-1)^((h*2^n-2)*(79-1)/4) == -1
* 83 mod 4 == -1 implies (-1)^((h*2^n-2)*(83-1)/4) == -1
* J(D, h*2^n-1) == J(83, h*2^n-1) * J(79, h*2^n-1)
* == J(h*2^n-1, 83) * (-1)^((h*2^n-2)*(83-1)/4) *
* J(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
* == J(h*2^n-1 mod 83, 83) * -1 *
* J(h*2^n-1 mod 79, 79) * -1
* We will also satisfy [condition 2]. Observe:
* (a^2 - b^2*D)/r == (79^2 - 1^1*6557)/316
* L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
* == L(79, h*2^n-1) * L(2^2, h*2^n-1)
* == L(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
* == L(h*2^n-1 mod 79, 79) * -1
* == J(h*2^n-1 mod 79, 79) * -1
if (jacobi( ((h*pmod(2,n%(79-1),79)-1)%79), 79 ) == -1 &&
jacobi( ((h*pmod(2,n%(83-1),83)-1)%83), 83 ) == 1) {
/* return the associated v(1)=81 */
/* no quick and dirty v(1), so return -1 */
* ldebug - print a debug statement
* funct name of calling function
print "DEBUG:", funct:":", str;
print "lucas(h, n) defined";
projects / unix-history / blob