Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/unix/volume26/calc/part13

Text File | 1992-05-09 | 29.9 KB | 1,063 lines

Newsgroups: comp.sources.unix
From: [e-mail] (David I. Bell)
Subject: v26i039: CALC - An arbitrary precision C-like calculator, Part13/21
Sender: [e-mail]
Approved: [e-mail]

Submitted-By: [e-mail] (David I. Bell)
Posting-Number: Volume 26, Issue 39
Archive-Name: calc/part13

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