apcalc-common 2.12.5.0-1build1 / usr / share / calc / lucas.cal

/*
 * lucas - perform a Lucas primality test on h*2^n-1
 *
 * Copyright (C) 1999  Landon Curt Noll
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 * @(#) $Revision: 30.2 $
 * @(#) $Id: lucas.cal,v 30.2 2013/09/27 08:58:46 chongo Exp $
 * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lucas.cal,v $
 *
 * Under source code control:	1990/05/03 16:49:51
 * File existed as early as:	1990
 *
 * chongo <was here> /\oo/\	http://www.isthe.com/chongo/
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * NOTE: This is a standard calc resource file.  For information on calc see:
 *
 *	http://www.isthe.com/chongo/tech/comp/calc/index.html
 *
 * To obtain your own copy of calc, see:
 *
 *	http://www.isthe.com/chongo/tech/comp/calc/calc-download.html
 */

/*
 * HISTORICAL NOTE:
 *
 * 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:
 *
 *				  216193
 *			391581 * 2	-1
 *
 * 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 31 Dec 1995, this prime was the 3rd
 * largest known prime and the largest known non-Mersenne prime.
 *
 * The same team also discovered the following twin prime pair:
 *
 *			   11235		   11235
 *		1706595 * 2	-1	1706595 * 2	+1
 *
 * At the time of discovery, this was the largest known twin prime pair.
 *
 * See:
 *
 *	http://www.isthe.com/chongo/tech/math/prime/amdahl6.html
 *
 * for more information on the Amdahl 6 group.
 *
 * NOTE: Both largest known and largest known twin prime records have been
 *	 broken.  Rather than update this file each time, I'll just
 *	 congratulate the finders and encourage others to try for
 *	 larger finds.	Records were made to be broken afterall!
 */

/* 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
 *	n mod 2^x == 0  for some value of x, say > 7 or more
 *	h*2^n-1 is not divisible by any small prime < 2^40
 *	0 < h < 2^39
 *	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 200000, 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
 * just a trickle.
 *
 * 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
 * use of the following:
 *
 *	if
 *		y = 2^x mod q
 *	then
 *		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" */

/*
 * lucas - lucas primality test on h*2^n-1
 *
 * ABOUT THE TEST:
 *
 * This routine will perform a primality test on h*2^n-1 based on
 * the mathematics of Lucas, Lehmer and Riesel.	 One should read
 * the following article:
 *
 * Ref1:
 *	"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:
 *
 * Ref2:
 *	"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:
 *
 * Ref3:
 *	"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:
 *
 *	n >= 2
 *	h >= 1
 *	h < 2^n
 *	h mod 2 == 1
 *
 * A few misc notes:
 *
 * 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.
 *
 * In real life, you would eliminate candidates by checking for
 * divisibility by a prime much larger than 257 (perhaps as high
 * as 2^39).
 *
 * 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
 * do make this so.
 *
 * input:
 *	h    the h as in h*2^n-1
 *	n    the n as in h*2^n-1
 *
 * returns:
 *	1 => h*2^n-1 is prime
 *	0 => h*2^n-1 is not prime
 *     -1 => a test could not be formed, or h >= 2^n, h <= 0, n <= 0
 */
define
lucas(h, n)
{
	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 */

	/*
	 * check arg types
	 */
	if (!isint(h)) {
		ldebug("lucas", "h is non-int");
		quit "FATAL: bad args: h must be an integer";
	}
	if (!isint(n)) {
		ldebug("lucas", "n is non-int");
		quit "FATAL: bad args: n must be an integer";
	}

	/*
	 * reduce h if even
	 *
	 * we will force h to be odd by moving powers of two over to 2^n
	 */
	oldh = h;
	oldn = n;
	shiftdown = fcnt(h,2);	/* h % 2^shiftdown == 0, max shiftdown */
	if (shiftdown > 0) {
		h >>= shiftdown;
		n += shiftdown;
	}

	/*
	 * enforce the 0 < h < 2^n rule
	 */
	if (h <= 0 || n <= 0) {
		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");
		return -1;
	}
	if (highbit(h) >= n) {
		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");
		return -1;
	}

	/*
	 * catch the degenerate case of h*2^n-1 == 1
	 */
	if (h == 1 && n == 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
	 */
	if (h == 1 && n == 2) {
		ldebug("lucas", "prime: h == 1 && n == 2");
		return 1;	/* 1*2^2-1 == 3 is prime */
	}
	if (h == 3 && n == 2) {
		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 (h == 1) {
		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 (h == 3) {
		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 */
		}
	}
	if (h == 5 && n == 4) {
		ldebug("lucas", "prime: 79 is prime");
		return 1;		/* 79 is prime */
	}
	if (h == 7 && n == 5) {
		ldebug("lucas", "prime: 223 is prime");
		return 1;		/* 223 is prime */
	}
	if (h == 15 && n == 4) {
		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
	 */
	testval = h*2^n - 1;
	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");
		return 0;
	}
	/*
	 * 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");
		return 0;
	}
	/*
	 * check for prime 47 <= factors < 257, if h*2^n-1 is large
	 * 2^282 > pfact(256)/pfact(46) > 2^281
	 */
	bits = highbit(testval);
	if (bits >= 281) {
		if (pprod256 <= 0) {
			pprod256 = pfact(256)/pfact(46);
		}
		if (gcd(testval, pprod256) > 1) {
			/* a small 47 <= prime < 257 divides h*2^n-1 */
			ldebug("lucas",\
			    "not-prime: 47<=prime<257 divides h*2^n-1");
			return 0;
		}
	}

	/*
	 * try to compute u(0)
	 *
	 * 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
	 * the v(1) into u(0).
	 *
	 * 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);
	if (v1 < 0) {
		/* failure to test number */
		print "unable to compute v(1) for", h : "*2^" : n : "-1";
		ldebug("lucas", "unknown: no v(1)");
		return -1;
	}
	u = gen_u0(h, n, v1);

	/*
	 * compute u(n-2)
	 */
	for (i=3; i <= n; ++i) {
		/* u = (u^2 - 2) % testval; */
		u = hnrmod(u^2 - 2, h, n, -1);
	}

	/*
	 * return 1 if prime, 0 is not prime
	 */
	if (u == 0) {
		ldebug("lucas", "prime: end of test");
		return 1;
	} else {
		ldebug("lucas", "not-prime: end of test");
		return 0;
	}
}

/*
 * gen_u0 - determine the initial Lucas sequence for h*2^n-1
 *
 * According to Ref1, Theorem 5:
 *
 *	u(0) = alpha^h + alpha^(-h)
 *
 * Now:
 *
 *	v(x) = alpha^x + alpha^(-x)	(Ref1, bottom of page 872)
 *
 * Therefore:
 *
 *	u(0) = v(h)
 *
 * 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) = v(n)^2 - 2
 *	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).
 *
 * input:
 *	h	- h as in h*2^n-1
 *	n	- n as in h*2^n-1
 *	v1	- gen_v1(h,n)		(see function below)
 *
 * returns:
 *	u(0)	- initial value for Lucas test on h*2^n-1
 *	-1	- failed to generate u(0)
 */
define
gen_u0(h, n, 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 */
	local i;

	/*
	 * check arg types
	 */
	if (!isint(h)) {
		quit "bad args: h must be an integer";
	}
	if (!isint(n)) {
		quit "bad args: n must be an integer";
	}
	if (!isint(v1)) {
		quit "bad args: v1 must be an integer";
	}
	if (v1 <= 0) {
		quit "bogus arg: v1 is <= 0";
	}

	/*
	 * enforce the h > 0 and n >= 2 rules
	 */
	if (h <= 0 || n < 1) {
		quit "reduced args violate the rule: 0 < h < 2^n";
	}
	hbits = highbit(h);
	if (hbits >= n) {
		quit "reduced args violate the rule: 0 < h < 2^n";
	}

	/*
	 * build up u2 based on the reversed bits of h
	 */
	/* setup for bit loop */
	r = v1;
	s = (r^2 - 2);

	/*
	 * 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.
	 */
	if (h == 1) {
		ldebug("gen_u0", "quick h == 1 case");
		/* return r%(h*2^n-1); */
		return hnrmod(r, h, n, -1);
	}

	/* cycle from second highest bit to second lowest bit of h */
	for (i=hbits-1; i > 0; --i) {

		/* bit(i) is 1 */
		if (bit(h,i)) {

			/* compute v(2n+1) = v(r+1)*v(r)-v1 */
			/* r = (r*s - v1) % (h*2^n-1); */
			r = hnrmod((r*s - v1), h, n, -1);

			/* compute v(2n+2) = v(r+1)^2-2 */
			/* s = (s^2 - 2) % (h*2^n-1); */
			s = hnrmod((s^2 - 2), h, n, -1);

		/* bit(i) is 0 */
		} else {

			/* compute v(2n+1) = v(r+1)*v(r)-v1 */
			/* s = (r*s - v1) % (h*2^n-1); */
			s = hnrmod((r*s - v1), h, n, -1);

			/* compute v(2n) = v(r)^-2 */
			/* r = (r^2 - 2) % (h*2^n-1); */
			r = hnrmod((r^2 - 2), h, n, -1);
		}
	}

	/* we know that h is odd, so the final bit(0) is 1 */
	/* r = (r*s - v1) % (h*2^n-1); */
	r = hnrmod((r*s - v1), h, n, -1);

	/* compute the final u2 return value */
	return r;
}

/*
 * 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
 * of {note 6}.
 *
 * It should be noted that the D values all satisfy the selection values
 * as outlined in the gen_v1() function comments.  That is:
 *
 *	   D == P*(2^f)*(3^g)
 *
 * 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
 * is of this form.
 */
quickmax = 8;
mat d_qval[quickmax];
mat v1_qval[quickmax];
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
 *
 * This function assumes:
 *
 *	n > 2			(n==2 has already been eliminated)
 *	h mod 2 == 1
 *	h < 2^n
 *	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.
 *
 ***
 *
 * Case 1:	(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)
 *
 *	h mod 6 == +/-1
 *	h*2^n-1 mod 3 != 0
 *
 * which translates, gives the functions assumptions, into the condition:
 *
 *	h mod 3 != 0
 *
 * 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)
 *
 * If we let:
 *
 *	alpha = (2+sqrt(3))
 *
 * then
 *
 *	u(0) = v(h)
 *
 * so we simply return
 *
 *	v(1) = alpha^1 + alpha^(-1)
 *	     = (2+sqrt(3)) + (2-sqrt(3))
 *	     = 4
 *
 ***
 *
 * Case 2:	(h mod 3 == 0)
 *
 * 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:
 *
 *	    r = abs(a^2 - b^2*D)
 *	    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)
 *
 * thus since:
 *
 *	u(0) = v(h)
 *
 * so we want to return:
 *
 *	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
 *
 * when
 *
 *	p is prime
 *	p mod 2 == 1
 *	gcd(a,p) == 1
 *
 * The value x is a quadratic residue of y if there exists some integer z
 * such that:
 *
 *	z^2 mod y == x
 *
 * 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
 *
 * when
 *
 *	y mod 2 == 1
 *	gcd(x,y) == 1
 *
 * 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}
 *
 *	p is prime
 *	p mod 2 == 1
 *	gcd(a,p) == 1
 *
 * It can be shown that (see Ref3):
 *
 *	L(a,p) == L(a mod p, p)					{note 1}
 *	L(z^2, p) == 1						{note 2}
 *
 * From Ref2, table 32:
 *
 *	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 is square free
 *	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.
 *
 * Take P mod 4 == 1:
 *
 *	P mod 4 == 1  implies  (-1)^((h*2^n-2)*(P-1)/4) == 1
 *
 * Thus:
 *
 *	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)
 *		      == J(h*2^n-1 mod P, P)
 *
 * Take P mod 4 == -1:
 *
 *	P mod 4 == -1  implies	(-1)^((h*2^n-2)*(P-1)/4) == -1
 *
 * Thus:
 *
 *	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
 * to 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
 *
 * 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:
 *
 *	(a^2 - b^2*D)/r == -1
 *
 * Therefore in order for [condition 2] to be met, we must show that:
 *
 *	L(r, h*2^n-1) == 1
 *
 * 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)
 *		      == 1					       {note 2}
 *
 * 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.
 *
 * Take Q mod 4 == 1:
 *
 *	Q mod 4 == 1  implies  (-1)^((h*2^n-2)*(Q-1)/4) == 1
 *
 * Thus:
 *
 *	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)
 *		      == J(h*2^n-1 mod Q, Q)
 *
 * Take Q mod 4 == -1:
 *
 *	Q mod 4 == -1  implies	(-1)^((h*2^n-2)*(Q-1)/4) == -1
 *
 * Thus:
 *
 *	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
 * to cases are 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
 *
 ***
 *
 * In conclusion, we can compute v(1) by attempting to do the following:
 *
 * h mod 3 != 0
 *
 *     we return:
 *
 *	   v(1) == 4
 *
 * h mod 3 == 0
 *
 *     define:
 *
 *	   r == abs(a^2 - b^2*D)
 *	   alpha == (a + b*sqrt(D))^2/r
 *
 *     we return:
 *
 *	   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 is square free
 *
 *	   D == P*(2^f)*(3^g)		(P is prime>2, f,g == 0 or 1)
 *
 *	   (a^2 - b^2*D)/r == -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
 *     is not too bad.
 *
 ***
 *
 * input:
 *	h	h as in h*2^n-1
 *	n	n as in h*2^n-1
 *
 * output:
 *	returns v(1), or -1 is there is no quick way
 */
define
gen_v1(h, n)
{
	local d;	/* the 'D' value to try */
	local val_mod;	/* h*2^n-1 mod 'D' */
	local i;

	/*
	 * check for case 1
	 */
	if (h % 3 != 0) {
		/* v(1) is easy to compute */
		return 4;
	}

	/*
	 * 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) {

		/* grab our 'D' value */
		d = d_qval[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'
		 * else
		 *	(h*2^n-1) mod 'D' must be a quadratic residue of 'D'
		 */
		if (d%4 == 1) {
			/* 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 */
				return v1_qval[i];
			}
		} else {
			/* 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 */
				return v1_qval[i];
			}
		}
	}

	/*
	 * This is an example of a more complex proof construction.
	 * The code above will not be able to find the v(1) for:
	 *
	 *			81*2^81-1
	 *
	 * We will check with:
	 *
	 *	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
	 *		      ==  1 * -1 *
	 *			 -1 * -1
	 *		      == -1
	 *
	 * We will also satisfy [condition 2].	Observe:
	 *
	 *	(a^2 - b^2*D)/r == (79^2 - 1^1*6557)/316
	 *			== -1
	 *
	 *	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(79, h*2^n-1) * 1
	 *		      == L(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
	 *		      == L(h*2^n-1, 79) * -1
	 *		      == L(h*2^n-1 mod 79, 79) * -1
	 *		      == J(h*2^n-1 mod 79, 79) * -1
	 *		      == -1 * -1
	 *		      == 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 */
		return 81;
	}

	/* no quick and dirty v(1), so return -1 */
	return -1;
}

/*
 * ldebug - print a debug statement
 *
 * input:
 *	funct	name of calling function
 *	str	string to print
 */
define
ldebug(funct, str)
{
	if (config("resource_debug") & 8) {
		print "DEBUG:", funct:":", str;
	}
	return;
}