/*
 * lucas - perform a Lucas primality test on h*2^n-1
 *
 * Copyright (C) 1999,2017,2018,2021  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.
 *
 * 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/
 */

/*
 * For a general tutorial on how to find a new largest known prime, see:
 *
 *	http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf
 *
 * Also see the reference code, available both C and go:
 *
 *	https://github.com/arcetri/goprime
 */

/*
 * 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 after all!
 */

/*
 * ON GAINING A WORLD RECORD:
 *
 * For a general tutorial on how to find a new largest known prime, see:
 *
 *	http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf
 *
 * 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.
 *
 * NOTE: While the PFA Fast Fourier Transform and Winograd's radix FFTs
 *	 might have been optimal for the Amdahl 6 team at the time,
 *	 they might not be optimal for your CPU architecture.  See
 *	 the above mentioned tutorial for information on better
 *	 methods of performing multiplications and squares of very
 *	 large numbers.
 *
 * 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.
 *
 * NOTE: Today, for world record sized h*2^n-1 primes, one might
 *	 search for factors < 2^46 or more.  By excluding h*2^n-1
 *	 with prime factors < 2^46, where h*2^n-1 is a bit larger
 *	 than the largest known prime, one may exclude about 96.5%
 *	 of candidates that have "small" prime factors.
 */

pprod256 = 0;	/* 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
 *
 *	http://www.ams.org/journals/mcom/1969-23-108/S0025-5718-1969-0262163-1/
 *			   S0025-5718-1969-0262163-1.pdf
 *
 *	NOTE: Join the above two lines for the complete URL of the paper.
 *
 * 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.
 *
 * An excellent 5-page paper by Oystein J. Rodseth (we apologize that the
 * ASCII character set does not allow us to spell his name with the
 * umlaut marks on the O's):
 *
 * NOTE: The original Amdahl 6 method predates the publication of Ref4.
 *	 The gen_v1() function used by lucas() uses the Ref4 method.
 *	 See the 'Amdahl 6 legacy code' section below for the original
 *	 method of generating v(1).
 *
 * Ref4:
 *
 * 	"A note on primality tests for N = h*2^n-1", by Oystein J. Rodseth,
 *	Department of Mathematics, University of Bergen, BIT Numerical
 *	Mathematics. 34 (3): pp 451-454.
 *
 *	http://folk.uib.no/nmaoy/papers/luc.pdf
 *
 * This test is performed as follows:	(see Ref1, Theorem 5)
 *
 *	a) generate u(2)		(see the function gen_u2() below)
 *					(NOTE: some call this u(0))
 *
 *	b) generate u(n) 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) == 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 miscellaneous 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    h as in h*2^n-1 (must be >= 1)
 *	n    n as in h*2^n-1 (must be >= 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(2) */
	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)) {
		quit "FATAL: bad args: h must be an integer";
	}
	if (h < 1) {
		quit "FATAL: bad args: h must be an integer >= 1";
	}
	if (!isint(n)) {
		quit "FATAL: bad args: n must be an integer";
	}
	if (n < 1) {
		quit "FATAL: bad args: n must be an integer >= 1";
	}

	/*
	 * 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 */
	}

	/*
	 * Verify that h*2^n-1 is not a multiple of 3
	 *
	 * The case for h*2^n-1 == 3 is handled above.
	 */
	if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) {
		/* no need to test h*2^n-1, it is a multiple of 3 */
		ldebug("lucas","not-prime: != 3 and is a multiple of 3");
		return 0;
	}

	/*
	 * Avoid any numbers divisible by small primes
	 */
	/*
	 * check for 5 <= prime factors < 31
	 * pfact(30)/6 = 1078282205
	 */
	testval = h*2^n - 1;
	if (gcd(testval, 1078282205) > 1) {
		/* a small 5 <= prime < 31 divides h*2^n-1 */
		ldebug("lucas",\
		    "not-prime: a small 5<=prime<31 divides h*2^n-1");
		return 0;
	}
	/*
	 * check for 31 <= prime factors < 53
	 * pfact(52)/pfact(30) = 95041567
	 */
	if (gcd(testval, 95041567) > 1) {
		/* a small 31 <= prime < 53 divides h*2^n-1 */
		ldebug("lucas","not-prime: 31<=prime<53 divides h*2^n-1");
		return 0;
	}
	/*
	 * check for prime 53 <= factors < 257, if h*2^n-1 is large
	 * 2^276 > pfact(256)/pfact(52) > 2^275
	 */
	bits = highbit(testval);
	if (bits >= 275) {
		if (pprod256 <= 0) {
			pprod256 = pfact(256)/pfact(52);
		}
		if (gcd(testval, pprod256) > 1) {
			/* a small 53 <= prime < 257 divides h*2^n-1 */
			ldebug("lucas",\
			    "not-prime: 53<=prime<257 divides h*2^n-1");
			return 0;
		}
	}

	/*
	 * try to compute u(2)		(NOTE: some call this 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_u2() to convert
	 * the v(1) into u(2).
	 *
	 * If gen_v1() returns a negative value, then we failed to
	 * generate a test for h*2^n-1.	 The legacy function,
	 * legacy_gen_v1() used by the Amdahl 6 could have returned
	 * -1. The new gen_v1() based on the method outlined in Ref4
	 * will never return -1 if h*2^n-1 is not a multiple of 3.
	 * Because the "multiple of 3" case is handled above, the
	 * call below to gen_v1() will never return -1.
	 */
	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_u2(h, n, v1);

	/*
	 * compute u(n)			(NOTE: some call this 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_u2 - determine the initial Lucas sequence for h*2^n-1
 *
 * Historically many start the Lucas sequence with u(0).
 * Some, like the author of this code, prefer to start
 * with U(2).  This is so one may say:
 *
 *	2^p-1 is prime if u(p) = 0 mod 2^p-1
 * or:
 *	h*2^p-1 is prime if u(p) = 0 mod h*2^p-1
 *
 * According to Ref1, Theorem 5:
 *
 *	u(2) = alpha^h + alpha^(-h)	(NOTE: Ref1 calls it u(0))
 *
 * Now:
 *
 *	v(x) = alpha^x + alpha^(-x)	(Ref1, bottom of page 872)
 *
 * Therefore:
 *
 *	u(2) = v(h)			(NOTE: Ref1 calls it u(0))
 *
 * 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	(must be >= 1)
 *	n	- n as in h*2^n-1	(must be >= 1)
 *	v1	- gen_v1(h,n)		(must be >= 3) (see function below)
 *
 * returns:
 *	u(2)	- initial value for Lucas test on h*2^n-1
 *	-1	- failed to generate u(2)
 */
define
gen_u2(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 oldh;		/* pre-reduced h */
	local oldn;		/* pre-reduced n */
	local i;

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

	/*
	 * 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 h > 0 and n >= 2 rules
	 */
	if (h <= 0 || n < 1) {
		print "	   ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
		quit "reduced args violate the rule: 0 < h < 2^n";
	}
	hbits = highbit(h);
	if (hbits >= n) {
		print "	   ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":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_u2", "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;
}

/*
 * gen_u0 - determine the initial Lucas sequence for h*2^n-1
 *
 * Historically many start the Lucas sequence with u(0).
 * Some, like the author of this code, prefer to start
 * with u(2).  This is so one may say:
 *
 *	2^p-1 is prime if u(p) = 0 mod 2^p-1
 * or:
 *	h*2^n-1 is prime if U(n) = 0 mod h*2^n-1
 *
 * For those using the old code with gen_u0(), we
 * simply call gen_u2() instead.
 *
 * See the function gen_u2() for details.
 *
 * input:
 *	h	- h as in h*2^n-1	(must be >= 1)
 *	n	- n as in h*2^n-1	(must be >= 1)
 *	v1	- gen_v1(h,n)		(see function below)
 *
 * returns:
 *	u(2)	- initial value for Lucas test on h*2^n-1
 *	-1	- failed to generate u(2)
 */
define
gen_u0(h, n, v1)
{
	return gen_u2(h, n, v1);
}

/*
 * rodseth_xhn - determine if v(1) == x for h*2^n-1
 *
 * For a given h*2^n-1, v(1) == x if:
 *
 *	jacobi(x-2, h*2^n-1) == 1		(Ref4, condition 1) part 1
 *	jacobi(x+2, h*2^n-1) == -1		(Ref4, condition 1) part 2
 *
 * Now when x-2 <= 0:
 *
 *	jacobi(x-2, h*2^n-1) == 0
 *
 * because:
 *
 *	jacobi(x,y) == 0                        if x <= 0
 *
 * So for (Ref4, condition 1) part 1 to be true:
 *
 *	x-2 > 0
 *
 * And therefore:
 *
 *	x > 2
 *
 * input:
 *	x	potential v(1) value
 *	h	h as in h*2^n-1 (h must be odd >= 1)
 *	n	n as in h*2^n-1	(must be >= 1)
 *
 * returns:
 *	1	if v(1) == x for h*2^n-1
 *	0	otherwise
 */
define
rodseth_xhn(x, h, n)
{
	local testval;		/* h*2^n-1 */

	/*
	 * check arg types
	 */
	if (!isint(h)) {
		quit "bad args: h must be an integer";
	}
	if (iseven(h)) {
		quit "bad args: h must be an odd integer";
	}
	if (h < 1) {
		quit "bad args: h must be an integer >= 1";
	}
	if (!isint(n)) {
		quit "bad args: n must be an integer";
	}
	if (n < 1) {
		quit "bad args: n must be an integer >= 1";
	}
	if (!isint(x)) {
		quit "bad args: x must be an integer";
	}

	/*
	 * firewall
	 */
	if (x <= 2) {
		return 0;
	}

	/*
	 * Check for jacobi(x-2, h*2^n-1) == 1	(Ref4, condition 1) part 1
	 */
	testval = h*2^n-1;
	if (jacobi(x-2, testval) != 1) {
		return 0;
	}

	/*
	 * Check for jacobi(x+2, h*2^n-1) == -1	(Ref4, condition 1) part 2
	 */
	if (jacobi(x+2, testval) != -1) {
		return 0;
	}

	/*
	 * v(1) == x for this h*2^n-1
	 */
	return 1;
}

/*
 * Trial tables used by gen_v1()
 *
 * When h mod 3 == 0, according to Ref4 we need to find the first value X where:
 *
 *	jacobi(X-2, h*2^n-1) == 1		(Ref4, condition 1) part 1
 *	jacobi(X+2, h*2^n-1) == -1		(Ref4, condition 1) part 2
 *
 * We can show that X > 2.  See the comments in the rodseth_xhn(x,h,n) above.
 *
 * Some values of X satisfy more often than others. For example a large sample
 * of h*2^n-1, h odd multiple of 3, and large n (some around 1e4, some near 1e6,
 * others near 3e7) where the sample size was 66 973 365, here is the count of
 * the smallest value of X that satisfies conditions in Ref4, condition 1:
 *
 *      count  X
 *    ----------
 *    26791345 3
 *    17223016 5
 *     7829600 9
 *     6988774 11
 *     3301093 15
 *     1517149 17
 *      910346 21
 *      711791 29
 *      573403 20
 *      390395 27
 *      288637 35
 *      149751 36
 *      107733 39
 *       58743 41
 *       35619 45
 *       25052 32
 *       17775 51
 *       13031 44
 *        7563 56
 *        7540 49
 *        7060 59
 *        4407 57
 *        2948 65
 *        2502 55
 *        2388 69
 *        2094 71
 *         689 77
 *         626 81
 *         491 66
 *         426 95
 *         219 80
 *         203 67
 *         185 84
 *         152 99
 *         127 72
 *         102 74
 *          98 87
 *          67 90
 *          55 104
 *          48 101
 *          32 105
 *          17 109
 *          16 116
 *          15 111
 *          13 92
 *          12 125
 *           7 129
 *           3 146
 *           2 140
 *           2 120
 *           1 165
 *           1 161
 *           1 155
 *
 * It is important that we select the smallest possible v(1).  While testing
 * various values of X for V(1) is fast, using larger than necessary values
 * of V(1) of can slow down calculating V(h).
 *
 * The above distribution was found to hold fairly well over many values of
 * odd h that are also a multiple of 3 and for many values of n where h < 2^n.
 *
 * For example for in a sample size of 1000000 numbers of the form h*2^n-1
 * where h is an odd multiple of 3, 12996351 <= h <= 13002351,
 * 4331116 <= n <= 4332116, these are the smallest v(1) values that were found:
 *
 *  smallest  percentage
 *	v(1)     used
 *  --------  ---------
 *	  3   40.0000 %
 *	  5   25.6833 %
 *	  9   11.6924 %
 *	 11   10.4528 %
 *	 15    4.8048 %
 *	 17    2.3458 %
 *	 21    1.3734 %
 *	 29    1.0527 %
 *	 20    0.8595 %
 *	 27    0.5758 %
 *	 35    0.4420 %
 *	 36    0.2433 %
 *	 39    0.1779 %
 *	 41    0.0885 %
 *	 45    0.0571 %
 *	 32    0.0337 %
 *	 51    0.0289 %
 *	 44    0.0205 %
 *	 49    0.0176 %
 *	 56    0.0137 %
 *	 59    0.0108 %
 *	 57    0.0053 %
 *	 65    0.0047 %
 *	 55    0.0045 %
 *	 69    0.0031 %
 *	 71    0.0024 %
 *	 66    0.0011 %
 *	 95    0.0008 %
 *	 81    0.0008 %
 *	 77    0.0006 %
 *	 72    0.0005 %
 *	 99    0.0004 %
 *	 80    0.0003 %
 *	 74    0.0003 %
 *	 84    0.0002 %
 *	 67    0.0002 %
 *	 87    0.0001 %
 *	104    0.0001 %
 *	129    0.0001 %
 *
 * However, a case can be made for considering only odd values for v(1)
 * candidates.  When h * 2^n-1 is prime and h is an odd multiple of 3,
 * a smallest v(1) that is even is extremely rate.  Of the list of 146553
 * known primes of the form h*2^n-1 when h is an odd a multiple of 3,
 * none has an smallest v(1) that was even.
 *
 * See:
 *
 *	https://github.com/arcetri/verified-prime
 *
 * for that list of 146553 known primes of the form h*2^n-1.
 *
 * That same example for in a sample size of 1000000 numbers of the
 * form h*2^n-1 where h is an odd multiple of 3, 12996351 <= h <= 13002351,
 * 4331116 <= n <= 4332116, these are the smallest odd v(1) values that were
 * found:
 *
 *  smallest    percentage
 *  odd v(1)       used
 *  --------    ---------
 *	  3	40.0000 %
 *	  5	25.6833 %
 *	  9	11.6924 %
 *	 11	10.4528 %
 *	 15	 4.8048 %
 *	 17	 2.3458 %
 *	 21	 1.6568 %
 *	 29	 1.6174 %
 *	 35	 0.4529 %
 *	 27	 0.3546 %
 *	 39	 0.3470 %
 *	 41	 0.2159 %
 *	 45	 0.1173 %
 *	 31	 0.0661 %
 *	 51	 0.0619 %
 *	 55	 0.0419 %
 *	 59	 0.0250 %
 *	 49	 0.0170 %
 *	 69	 0.0110 %
 *	 65	 0.0098 %
 *	 71	 0.0078 %
 *	 85	 0.0048 %
 *	 81	 0.0044 %
 *	 95	 0.0038 %
 *	 99	 0.0021 %
 *	125	 0.0009 %
 *	 57	 0.0007 %
 *	111	 0.0005 %
 *	 77	 0.0003 %
 *	165	 0.0003 %
 *	155	 0.0002 %
 *	129	 0.0002 %
 *	101	 0.0002 %
 *	 53	 0.0001 %
 *
 * Moreover when evaluating odd candidates for v(1), one may cache Jacobi
 * symbol evaluations to reduce the number of Jacobi symbol evaluations to
 * a minimum. For example, if one tests 5 and finds that the 2nd case fails:
 *
 *	jacobi(5+2, h*2^n-1) != -1
 *
 * Then if one is later testing 9, the Jacobi symbol value for the first
 * 1st case:
 *
 *	jacobi(7-2, h*2^n-1)
 *
 * is already known.
 *
 * Without Jacobi symbol value caching, it requires on average
 * 4.851377 Jacobi symbol evaluations.  With Jacobi symbol value caching
 * cacheing, an average of 4.348820 Jacobi symbol evaluations is needed.
 *
 * Given this information, when odd h is a multiple of 3 we try, in order,
 * these odd values of X:
 *
 *	3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59,
 *	69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129,
 *	101, 83, 165, 155, 149, 141, 121, 109
 *
 * And stop on the first value of X where:
 *
 *      jacobi(X-2, h*2^n-1) == 1
 *      jacobi(X+2, h*2^n-1) == -1
 *
 * Less than 1 case out of 1000000 will not be satisfied by the above list.
 * If no value in that list works, we start simple search starting with X = 167
 * and incrementing by 2 until a value of X is found.
 *
 * The x_tbl[] matrix contains those values of X to try in order.
 * If all x_tbl_len fail to satisfy Ref4 condition 1 (this happens less than
 * 1 in 1000000 cases), then we begin a linear search of odd values starting at
 * next_x until we find a proper X value.
 */
x_tbl_len = 42;
mat x_tbl[x_tbl_len];
x_tbl = {
    3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59,
    69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129,
    101, 83, 165, 155, 149, 141, 121, 109
};
next_x = 167;	/* must be 2 more than the largest value in x_tbl[] */

/*
 * 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.
 *
 * 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(2) = (2+sqrt(3))^h + (2-sqrt(3))^h	 (see Ref1, page 869)
 *	     = (2+sqrt(3))^h + (2+sqrt(3))^(-h)	 (NOTE: some call this u(2))
 *
 * and since Ref1, Theorem 5 states:
 *
 *	u(2) = alpha^h + alpha^(-h)		 (NOTE: some call this u(2))
 *	r = abs(2^2 - 1^2*3) = 1
 *
 * where these values work for Case 1:		 (h mod 3 != 0)
 *
 *	a = 1
 *	b = 2
 *	D = 1
 *
 * Now at the bottom of Ref1, page 872 states:
 *
 *	v(x) = alpha^x + alpha^(-x)
 *
 * If we let:
 *
 *	alpha = (2+sqrt(3))
 *
 * then
 *
 *	u(2) = v(h)				 (NOTE: some call this u(2))
 *
 * so we can always return
 *
 *	v(1) = alpha^1 + alpha^(-1)
 *	     = (2+sqrt(3)) + (2-sqrt(3))
 *	     = 4
 *
 * In 40% of the cases when h is not a multiple of 3, 3 is a valid value
 * for v(1).  We can test if 3 is a valid value for v(1) in this case:
 *
 *	if jacobi(1, h*2^n-1) == 1 and jacobi(5, h*2^n-1) == -1, then
 *	    v(1) = 3
 *	else
 *	    v(1) = 4
 *
 * NOTE: The above "if then else" works only of h is not a multiple of 3.
 *
 ***
 *
 * Case 2:	(h mod 3 == 0)
 *
 * For the case where h is a multiple of 3, we turn to Ref4.
 *
 * The central theorem on page 3 of that paper states that
 * we may set v(1) to the first value X that satisfies:
 *
 *	jacobi(X-2, h*2^n-1) == 1		(Ref4, condition 1)
 *	jacobi(X+2, h*2^n-1) == -1		(Ref4, condition 1)
 *
 *	NOTE: Ref4 uses P, which we shall refer to as X.
 *	      Ref4 uses N, which we shall refer to as h*2^n-1.
 *
 *	NOTE: Ref4 uses the term Legendre-Jacobi symbol, which
 *	      we shall refer to as the Jacobi symbol.
 *
 * 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 jacobi(a,b) and L(a,p):
 *
 * 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 a is a quadratic residue of b if there exists some integer z
 * such that:
 *
 *	z^2 mod b == a
 *
 * The Jacobi symbol jacobi(a,b) takes the value:
 *
 *	jacobi(a,b) == 1	=> b is not prime,
 *				   or a is a quadratic residue of b
 *	jacobi(a,b) == -1	=> a is NOT a quadratic residue of b
 *
 * when
 *
 *	b mod 2 == 1
 *	gcd(a,b) == 1
 *
 * It is worth noting for the Legendre symbol, in order for L(X+/-2,
 * h*2^n-1) to be defined, we must ensure that neither X-2 nor X+2 are
 * factors of h*2^n-1.  This is done by pre-screening h*2^n-1 to not
 * have small factors and keeping X+2 less than that small factor
 * limit.  It is worth noting that in lucas(h, n), we first verify
 * that h*2^n-1 does not have a factor < 257 before performing the
 * primality test.  So while X+/-2 < 257, we know that
 * gcd(X+/-2, h*2^n-1) == 1.
 *
 * Returning to the testing of conditions in Ref4, condition 1:
 *
 *	jacobi(X-2, h*2^n-1) == 1
 *	jacobi(X+2, h*2^n-1) == -1
 *
 * When such an X is found, we set:
 *
 *	v(1) = X
 *
 ***
 *
 * 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
 *
 *     we return:
 *
 *	    v(1) = X
 *
 *     where X > 2 in a integer such that:
 *
 *	    jacobi(X-2, h*2^n-1) == 1
 *	    jacobi(X+2, h*2^n-1) == -1
 *
 ***
 *
 * input:
 *	h	h as in h*2^n-1 (h must be odd >= 1)
 *	n	n as in h*2^n-1	(must be >= 1)
 *
 * output:
 *	returns v(1), or
 *		-1 when h*2^n-1 is a multiple of 3
 */
define
gen_v1(h, n)
{
	local x;	/* potential v(1) to test */
	local i;	/* x_tbl index */
	local v1m2;	/* X-2 1st case */
	local v1p2;	/* X+2 2nd case */
	local testval;	/* h*2^n-1 - value we are testing if prime */
	local mat cached_v1[next_x];	/* cached Jacobi symbol values or 0 */

	/*
	 * check arg types
	 */
	if (!isint(h)) {
		quit "bad args: h must be an integer";
	}
	if (iseven(h)) {
		quit "bad args: h must be an odd integer";
	}
	if (h < 1) {
		quit "bad args: h must be an integer >= 1";
	}
	if (!isint(n)) {
		quit "bad args: n must be an integer";
	}
	if (n < 1) {
		quit "bad args: n must be an integer >= 1";
	}

	/*
	 * pretest: Verify that h*2^n-1 is not a multiple of 3
	 */
	if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) {
		/* no need to test h*2^n-1, it is not prime */
		return -1;
	}

	/*
	 * Common Mersenne number case:
	 *
	 * For Mersenne numbers:
	 *
	 *	2^n-1
	 *
	 * we can use, 40% of the time, v(1) == 3.  However nearly all code that
	 * implements the Lucas-Lehmer test uses v(1) == 4.  Whenever for
	 * h != 0 mod 3, and particular the Mersenne number case of when h == 1:
	 *
	 *	1*2^n-1
	 *
	 * v(1) == 4 always works.  For this reason, we return 4 when h == 1.
	 */
	if (h == 1) {
		/* v(1) == 4 always works for the Mersenne number case */
		return 4;
	}

	/*
	 * check for Case 1:      (h mod 3 != 0)
	 */
	if (h % 3 != 0) {
		if (rodseth_xhn(3, h, n) == 1) {
			/* 40% of the time, 3 works when h mod 3 != 0 */
			return 3;
		} else {
			/* otherwise 4 always works when h mod 3 != 0 */
			return 4;
		}
	}

	/*
	 * What follow is Case 2:      (h mod 3 == 0)
	 */

	/*
	 * clear cache
	 */
	matfill(cached_v1, 0);

	/*
	 * We will look for x that satisfies conditions in Ref4, condition 1:
	 *
	 *	jacobi(X-2, h*2^n-1) == 1		part 1
	 *	jacobi(X+2, h*2^n-1) == -1		part 2
	 *
	 * NOTE: If we wanted to be super optimal, we would cache
	 *	 jacobi(X+2, h*2^n-1) that that when we increment X
	 *	 to the next odd value, the now jacobi(X-2, h*2^n-1)
	 *	 does not need to be re-evaluated.
	 */
	testval = h*2^n-1;
	for (i=0; i < x_tbl_len; ++i) {

		/*
		 * obtain the next test candidate
		 */
		x = x_tbl[i];

		/*
		 * Check x for condition 1 part 1
		 *
		 *	jacobi(x-2, h*2^n-1) == 1
		 */
		v1m2 = x-2;
		if (cached_v1[v1m2] == 0) {
			cached_v1[v1m2] = jacobi(v1m2, testval);
		}
		if (cached_v1[v1m2] != 1) {
			continue;
		}

		/*
		 * Check x for condition 1 part 2
		 *
		 *	jacobi(x+2, h*2^n-1) == -1
		 */
		v1p2 = x+2;
		if (cached_v1[v1p2] == 0) {
			cached_v1[v1p2] = jacobi(v1p2, testval);
		}
		if (cached_v1[v1p2] != -1) {
			continue;
		}

		/*
		 * found a x that satisfies Ref4 condition 1
		 */
		ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) +
			" v1= " + str(x) + " using tbl[ " +
			str(i) + " ]");
		return x;
	}

	/*
	 * We are in that rare case (less than 1 in 1 000 000) where none of the
	 * common X values satisfy Ref4 condition 1.  We start a linear search
	 * of odd values at next_x from here on.
	 */
	x = next_x;
	while (rodseth_xhn(x, h, n) != 1) {
		x += 2;
	}
	/* finally found a v(1) value */
	ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) +
	       " v1= " + str(x) + " beyond tbl");
	return x;
}

/*
 * 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;
}

/*
 ************************
 * Amdahl 6 legacy code *
 ************************
 *
 * NOTE: What follows is legacy code based on the method used by the
 *	 Amdahl 6 group:
 *
 *	    John Brown, Landon Curt Noll, Bodo Parady, Gene Smith,
 *	    Joel Smith and Sergio Zarantonello
 *
 * This method generated v(1) for nearly all values, except for a
 * few rare cases when h mod 3 == 0.  The code is NOT used by lucas.cal
 * above.  The gen_v1() function above is based on an improved method
 * outlined in Ref4.  That method generated v(1) for all h.
 *
 * The code below is kept for historical purposes only.  The functions
 * and global variables of the Amdahl 6 legacy code all begin with legacy_.
 */

/*
 * Trial tables used by legacy_gen_v1()
 *
 * When h mod 3 == 0, one needs particular values of D, a and b (see
 * legacy_gen_v1 documentation) in order to find a value of v(1).
 *
 * This table defines 'legacy_quickmax' possible tests to be taken in ascending
 * order.  The legacy_v1_qval[x] refers to a v(1) value from Ref1, Table 1.  A
 * related D value is found in legacy_d_qval[x].  All D values expect
 * legacy_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 legacy_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.
 */
legacy_quickmax = 8;
mat legacy_d_qval[legacy_quickmax];
mat legacy_v1_qval[legacy_quickmax];
legacy_d_qval[0] = 5;	legacy_v1_qval[0] = 3;	/* a=1	 b=1  r=4  */
legacy_d_qval[1] = 7;	legacy_v1_qval[1] = 5;	/* a=3	 b=1  r=12  D=21 */
legacy_d_qval[2] = 13;	legacy_v1_qval[2] = 11;	/* a=3	 b=1  r=4  */
legacy_d_qval[3] = 11;	legacy_v1_qval[3] = 20;	/* a=3	 b=1  r=2  */
legacy_d_qval[4] = 29;	legacy_v1_qval[4] = 27;	/* a=5	 b=1  r=4  */
legacy_d_qval[5] = 53;	legacy_v1_qval[5] = 51;	/* a=53	 b=1  r=4  */
legacy_d_qval[6] = 17;	legacy_v1_qval[6] = 66;	/* a=17	 b=1  r=1  */
legacy_d_qval[7] = 19;	legacy_v1_qval[7] = 74;	/* a=38	 b=1  r=2  */

/*
 * legacy_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(2) = (2+sqrt(3))^h + (2-sqrt(3))^h		(see Ref1, page 869)
 *	     = (2+sqrt(3))^h + (2+sqrt(3))^(-h)		(some call this u(0))
 *
 * and since Ref1, Theorem 5 states:
 *
 *	u(2) = alpha^h + alpha^(-h)
 *	r = abs(2^2 - 1^2*3) = 1
 *
 * where these values work for Case 1: 			(h mod 3 != 0)
 *
 *	a = 1
 *	b = 2
 *	D = 1
 *
 * Now at the bottom of Ref1, page 872 states:
 *
 *	v(x) = alpha^x + alpha^(-x)
 *
 * If we let:
 *
 *	alpha = (2+sqrt(3))
 *
 * then
 *
 *	u(2) = 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(2) = alpha^h + alpha^(-h)
 *
 * The bottom of Ref1, page 872 states:
 *
 *	v(x) = alpha^x + alpha^(-x)
 *
 * thus since:
 *
 *	u(2) = 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	(must be >= 1)
 *	n	n as in h*2^n-1	(must be >= 1)
 *
 * output:
 *	returns v(1), or -1 is there is no quick way
 */
define
legacy_gen_v1(h, n)
{
	local d;	/* the 'D' value to try */
	local val_mod;	/* h*2^n-1 mod 'D' */
	local i;

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

	/*
	 * 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 < legacy_quickmax; ++i) {

		/* grab our 'D' value */
		d = legacy_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 legacy_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 legacy_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;
}