lucas.cal

Calc Software Package for Number Calc

/*
 * 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.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.6 $
 * @(#) $Id: lucas.cal,v 29.6 2002/07/10 09:43:46 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/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	(h mod 2 != 0)
 *	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 mod rules
	 */
	if (h%2 == 0) {
		quit "h must not be even";
	}

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