lucas.cal

Calc Software Package for Number Calc

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