Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/unix/volume27/calc-2.9.0/part17

Text File | 1993-12-07 | 59.2 KB | 2,027 lines

Newsgroups: comp.sources.unix
From: [e-mail] (David I. Bell)
Subject: v27i144: calc-2.9.0 - arbitrary precision C-like programmable calculator, Part17/19
References: <[e-mail]>
Sender: [e-mail]
Approved: [e-mail]

Submitted-By: [e-mail] (David I. Bell)
Posting-Number: Volume 27, Issue 144
Archive-Name: calc-2.9.0/part17

#!/bin/sh
# this is part 17 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/lib/cryrand.cal continued
#
CurArch=17
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
     exit 1; fi
( read Scheck
  if test "$Scheck" != $CurArch
  then echo "Please unpack part $Scheck next!"
       exit 1;
  else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc2.9.0/lib/cryrand.cal"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/lib/cryrand.cal
X * values are kept.
X *
X * Unlike the first and second forms of this function, the shuffle
X * generator function is not seeded before or after the call.  The
X * current state is used to generate a new quadratic residue of n=p*q.
X *
X *
X * Third form of call:  scryrand(seed,len1,len2)
X *
X * In the third form, 'len1' and 'len2' guide this function in selecting
X * internally used prime numbers.  The larger the lengths, the longer
X * the time this function will take.  The impact on execution time of
X * cryrand() and random() may also be noticed, but not as much.
X *
X * If a length is '-1', then the default lengths (248 for len1, and 264
X * for len2) are used.  The call scryrand(0,-1,-1) recreates the initial
X * crypto state the slow and hard way.  (use scryrand(0) or srandom(0))
X *
X * This function can take a long time to call given reasonable values
X * of len1 and len2.  On a Pyramid R3000, the time to seed was:
X *
X *    Approx value    digits   seed time
X *      of len1+len2   in n=p*q       in sec
X *    ------------   --------       ------
X *          8           3         0.53
X *         16           5         0.54
X *         32          10         0.79
X *         64          20         1.17
X *        128          39         2.89
X *        200          61         4.68
X *        256          78         7.49
X *        322         100        12.47
X *        464         140        35.56
X *        512         155        53.57
X *        664         200        83.97
X *        830         250       122.93
X *        996         300       242.49
X *       1024         309       295.66
X *       1328         400       663.44
X *       1586         478      2002.10
X *       1660         500      1643.45  (Faster mult/square methods kick in
X *       1992         600      2885.81   in certain cases. Type  help config
X *       2048         617      1578.06   in calc for more details.)
X *
X *     NOTE: The small lengths above are given for comparison
X *           purposes and are NOT recommended for actual use.
X *
X *     NOTE: Generating crypto pseudo-random numbers is MUCH
X *           faster than seeding a crypto generator.
X *
X *     NOTE: This calc lib file is intended for demonstration
X *           purposes.  Writing a C program (with possible assembly
X *           or libmp assist) would produce a faster generator.
X *
X *
X * Fourth form of call:  scryrand(seed,ip,iq,ir)
X *
X * In the fourth form, 'ip', 'iq' and 'ir' serve as initial search
X * values for the two Blum primes 'p' and 'q' and an associated
X * quadratic residue 'r' respectively.  Unlike the 3rd form, where
X * lengths are given, the fourth form allows one to specify minimum
X * search values.
X *
X * The 'seed' value is interpreted as follows:
X *
X *   If seed > 0:
X *
X *    Seed and use the shuffle generator to generate 3 jump values
X *    that are in the range '[0,ip)', '[0,iq)' and '[0,ir)' respectively.
X *    Start searching for legal 'p', 'q' and 'r' values by adding
X *    the jump values to their respective argument values.
X *
X *   If seed == 0:
X *
X *    Start searching for legal 'p', 'q' and 'r' values from
X *    'ip', 'iq' and 'ir' respectively.
X *
X *    This form does not change/seed the other generators.
X *
X *   If seed == -1:
X *
X *    Let 'p' == 'ip', 'q' == 'iq' and 'r' == 'ir'.  Do not check
X *    if the value given are legal Blum primes or an associated
X *    quadratic residue respectively.
X *
X *    This form does not change/seed the other generators.
X *
X *    WARNING: No checks are performed on the args passed.
X *         Passing improper values will likely produce
X *         poor results, or worse!
X *
X *
X * It should be noted that calling scryrand() while using the default
X * primes took only 0.04 seconds.  Calling scryrand(0,-1,-1) took
X * 47.19 seconds.
X *
X * The paranoid, when giving explicit lengths, should keep in mind that
X * len1 and len2 are the largest powers of 2 that are less than the two
X * probable primes ('p' and 'q').  These two primes  will be used
X * internally to cryrand().  For simplicity, we refer to len1 and len2
X * as bit lengths, even though they are actually 1 less then the
X * minimum possible prime length.
X *
X * The actual lengths may exceed the lengths by slightly more than 3%.
X * Furthermore, part of the strength of this generator rests on the
X * difficultly to factor 'p*q'.  Thus one should select 'len1' and 'len2'
X * (from which 'p' and 'q' are selected) such that factoring a 'len1+len2'
X * bit number is difficult.
X *
X * Not being able to factor 'n=p*q' into 'p' and 'q' does not directly
X * improve the crypto generator.  On the other hand, it can't hurt.
X *
X * There is no limit on the size of a seed.  On the other hand,
X * extremely large seeds require large tables and long seed times.
X * Using a seed in the range of [2^64, 2^64 * 128!) should be
X * sufficient for most purposes.  An easy way to stay within this
X * range to to use seeds that are between 21 and 215 digits long, or
X * 64 to 780 bits long.
X *
X * NOTE: This function will clear any internally buffer bits.  See
X *     cryrand() for details.
X *
X * NOTE: This function seeds the Blum cryptographically strong
X *     pseudo-random number generator.
X */
Xdefine
Xscryrand(seed,len1,len2,arg4)
X{
X    local rval;        /* a temporary pseudo-random value */
X    local oldseed;    /* the previous seed */
X    local newres;    /* the new quad res */
X    local ip;        /* initial Blum prime 'p' search value */
X    local iq;        /* initial Blum prime 'q' search value */
X    local ir;        /* initial quadratic residue search value */
X    local rlim;        /* high quadratic res search value */
X
X    /*
X     * firewall - avoid bogus args and very trivial lengths
X     */
X    /* catch the case of no args - compute a different quadratic residue */
X    if (isnull(seed) && isnull(len1) && isnull(len2)) {
X
X    /* generate the next quadratic residue */
X    do {
X        newres = cryres(cryrand_p, cryrand_q);
X    } while (newres == cryrand_r);
X    cryrand_r = newres;
X    cryrand_exp = cryrand_r;
X
X    /* clear the internal bits */
X    cryrand_left = 0;
X    cryrand_bitcnt = 0;
X
X    /* return the current seed early */
X    return (cry_seed);
X    }
X    if (!isint(seed)) {
X    quit "bad arg: seed arg (1st) must be an integer";
X    }
X    if (param(0) == 4) {
X    if (seed < -1) {
X        quit "bad arg: with 4 args: a seed < -1 is reserved for future use";
X    }
X    } else {
X    if (4 && seed < 0) {
X        quit "bad arg: a seed arg (1st) < 0 is reserved for future use";
X    }
X    }
X
X    /*
X     * 4 arg case: select or search for 'p', 'q' and 'r' from given values
X     */
X    if (param(0) == 4) {
X
X    /* set initial values */
X    ip = len1;
X    iq = len2;
X    ir = arg4;
X
X    /*
X     * Unless prohibited by a seed if -1, force minimum values on
X     * 'ip', 'iq' and 'ir'.
X     */
X    if (seed >= 0) {
X        if (!isint(ip) || ip < 3) {
X           /* smallest Blum prime */
X           ip = 3;
X        }
X        if (!isint(iq) || iq < 3) {
X           /* smallest Blum prime */
X           iq = 3;
X        }
X        if (!isint(ir) || ir < 4) {
X           /* cryres() never selects a value < 4 */
X           ir = 4;
X        }
X    }
X
X    /*
X     * Determine our prime search points
X     */
X    if (seed > 0) {
X        /* add in a random value */
X        oldseed = srand(seed);
X        ip += rand(ip);
X        iq += rand(iq);
X    }
X
X    /*
X     * force ip <= iq
X     */
X    if (ip > iq) {
X        swap(ip, iq);
X    }
X
X    /*
X     * find the first Blum prime 'p'
X     */
X    if (seed >= 0) {
X        cryrand_p = nxtprime(ip,3,4);
X    } else {
X        /* seed == -1: assume 'ip' is a Blum prime */
X        cryrand_p = ip;
X    }
X
X    /*
X     * find the 2nd Blum prime 'q' > 'p', if needed
X     */
X    if (seed >= 0) {
X        if (iq <= cryrand_p) {
X        iq = cryrand_p + 2;
X        }
X        cryrand_q = nxtprime(iq,3,4);
X    } else {
X        /* seed == -1: assume 'iq' is a Blum prime */
X        cryrand_q = iq;
X    }
X
X    /* remember our p*q Blum prime product as well */
X    cryrand_n = cryrand_p*cryrand_q;
X
X    /*
X     * search for a quadratic residue
X     */
X    if (seed >= 0) {
X
X        /* add in a random value to 'ir' if seeded */
X        if (seed > 0) {
X        ir += rand(ir);
X        }
X
X        /*
X         * increment until we find a quadratic value
X         *
X         * p is prime and J(r,p) == 1  ==>  r is a quadratic residue of p
X         * q is prime and J(q,p) == 1  ==>  r is a quadratic residue of q
X         *
X         * J(r,p)*J(r,q) == J(r,p*q)
X         * J(r,p) and J(q,p) == 1       ==>  J(r,p*q) == 1
X         * J(r,p*q) == 1           ==>  r is a quadratic residue of p*q
X         *
X         * Try to avoid trivial quadratic residues by forcing the search
X         * over the interval [4, n-3).
X         */
X        rlim = cryrand_n-3;
X        /* if ir is too big, drop down to the bottom of the range */
X        if (ir >= rlim) {
X        ir = 4;
X        }
X        while (jacobi(ir,cryrand_p) != 1 || jacobi(ir,cryrand_q) != 1) {
X        /* if ir is too big, drop down to the bottom of the range */
X        if (++ir >= rlim) {
X            ir = 4;
X        }
X        }
X    }
X    cryrand_r = ir;
X
X    /*
X     * clear the previously unused cryrand bits & other misc setup
X     */
X    cryrand_left = 0;
X    cryrand_bitcnt = 0;
X    cryrand_exp = cryrand_r;
X
X    /*
X     * reseed the generator, if needed
X     *
X     * The crypto generator no longer needs the additive 55 and shuffle
X     * generators.  We however, restore the additive 55 and shuffle
X     * generators back to its seeded state in order to be sure that it
X     * will be left in the same state.
X     *
X     * This will make more reproducible, calls to the additive 55 and
X     * shuffle generator; or more reproducible, calls to this function
X     * without args.
X     */
X    if (seed > 0) {
X        ir = srand(seed);   /* ignore this return value */
X        return(oldseed);
X    } else {
X        /* no seed change, return the current seed */
X        return (cry_seed);
X    }
X    }
X
X    /*
X     * If not the 4 arg case:
X     *
X     * convert explicit -1 args into default values
X     * convert missing args into -1 values (use precomputed tables)
X     */
X    if ((isint(len1) && len1 != -1 && len1 < 4) ||
X    (isint(len2) && len2 != -1 && len2 < 4) ||
X    (!isint(len1) && isint(len2)) ||
X    (isint(len1) && !isint(len2))) {
X    quit "bad args: 2 & 3: if 2nd is given, must be -1 or ints > 3";
X    }
X    if (isint(len1) && len1 == -1) {
X    len1 = 248;    /* default len1 value */
X    }
X    if (isint(len2) && len2 == -1) {
X    len2 = 264;    /* default len2 value */
X    }
X    if (!isint(len1) && !isint(len2)) {
X    /* from here down, -1 means use precomputed values */
X    len1 = -1;
X    len2 = -1;
X    }
X
X    /*
X     * force len1 <= len2
X     */
X    if (len1 > len2) {
X    swap(len1, len2);
X    }
X
X    /*
X     * seed the generator
X     */
X    oldseed = srand(seed);
X
X    /*
X     * generate p and q Blum primes
X     *
X     * The Blum process requires the primes to be of the form 3 mod 4.
X     * We also generate n=p*q for future reference.
X     *
X     * We make sure that the lengths are the minimum lengths possible.
X     * We want some range to select a random prime from, so we
X     * go at least 3 bits higher, and as much as 3% plus 3 bits
X     * higher.  Since the section is a random, how high really
X     * does not matter that much, but we want to avoid going to
X     * an extreme to keep the execution time from getting too long.
X     *
X     * Finally, we generate a quadratic residue of n=p*q.
X     */
X    if (len1 > 0) {
X    /* generate a pseudo-random prime ~len1 bits long */
X    rval = rand(2^(len1-1), 2^((int(len1*1.03))+3));
X    cryrand_p = nxtprime(rval,3,4);
X    } else {
X    /* use precomputed 'p' value */
X    cryrand_p = cryrand_init_p;
X    }
X    if (len2 > 0) {
X    /* generate a pseudo-random prime ~len1 bits long */
X    rval = rand(2^(len2-1), 2^((int(len2*1.03))+3));
X    cryrand_q = nxtprime(rval,3,4);
X    } else {
X    /* use precomputed 'q' value */
X    cryrand_q = cryrand_init_q;
X    }
X
X    /*
X     * find the quadratic residue
X     */
X    cryrand_n = cryrand_p*cryrand_q;
X    cryrand_r = cryres(cryrand_p, cryrand_q);
X
X    /*
X     * clear the previously unused cryrand bits & other misc setup
X     */
X    cryrand_left = 0;
X    cryrand_bitcnt = 0;
X    cryrand_exp = cryrand_r;
X
X    /*
X     * reseed the generator
X     *
X     * The crypto generator no longer needs the additive 55 and shuffle
X     * generators.  We however, restore the additive 55 and shuffle
X     * generators back to its seeded state in order to be sure that it
X     * will be left in the same state.
X     *
X     * This will make more reproducible, calls to the additive 55 and
X     * shuffle generator; or more reproducible, calls to this function
X     * without args.
X     */
X    /* we do not care about this old seed */
X    rval = srand(seed);
X
X    /*
X     * return the old seed
X     */
X    return(oldseed);
X}
X
X
X/*
X * random - a cryptographically strong pseudo-random number generator
X *
X * usage:
X *    random()     - generate a pseudo-random integer >=0 and < 2^64
X *    random(a)     - generate a pseudo-random integer >=0 and < a
X *    random(a,b)    - generate a pseudo-random integer >=a and <= b
X *
X * returns:
X *    a large cryptographically strong pseudo-random number  (see usage)
X *
X * This function is just another interface to the crypto generator.
X * (see the cryrand() function).
X *
X * When no arguments are given, a pseudo-random number in the half open
X * interval [0,2^64) is produced.  This form is identical to calling
X * cryrand(64).
X *
X * When 1 argument is given, a pseudo-random number in the half open interval
X * [0,a) is produced.
X *
X * When 2 arguments are given, a pseudo-random number in the closed interval
X * [a,b] is produced.
X *
X * This generator uses the crypto to return a large pseudo-random number.
X *
X * The input values a and b, if given, must be integers.
X *
X * Internally, bits are produced log2(log2(n=p*q)) at a time.  If a
X * call to this function does not exhaust all of the collected bits,
X * the unused bits will be saved away and used at a later call.
X * Setting the seed via scryrand(), srandom() or cryrand(len,1)
X * will clear out all unused bits.
X *
X * NOTE: The BSD random() function returns only 31 bits, while we return 64.
X *
X * NOTE: This function is the Blum cryptographically strong
X *     pseudo-random number generator!
X */
Xdefine
Xrandom(a,b)
X{
X    local range;        /* we must generate [0,range) first */
X    local offset;        /* what to add to get a adjusted range */
X    local rangebits;        /* the number of bits in range */
X    local ret;            /* pseudo-random bit value */
X
X    /*
X     * setup and special cases
X     */
X    /* deal with the rand() case */
X    if (isnull(a) && isnull(b)) {
X    /* no args means return 64 bits */
X    return(cryrand(64));
X    }
X    /* firewall - args, if given must be in integers */
X    if (!isint(a) || (!isnull(b) && !isint(b))) {
X    quit "bad args: args, if given, must be integers";
X    }
X    /* convert rand(x) into rand(0,x-1) */
X    if (isnull(b)) {
X    /* convert call into a closed interval */
X    b = a-1;
X    a = 0;
X    /* firewall - rand(0) should act like rand(0,0) */
X    if (b == -1) {
X        return(0);
X    }
X    }
X    /* determine the range and offset */
X    if (a >= b) {
X    /* deal with the case of rand(a,a) */
X    if (a == b) {
X        /* not very random, but it is true! */
X        return(a);
X    }
X    range = a-b+1;
X    offset = b;
X    } else {
X    /* convert random(a,b), where a<b, into random(b,a) */
X    range = b-a+1;
X    offset = a;
X    }
X    rangebits = highbit(range-1)+1;
X    /*
X     * At this point, we seek a pseudo-random number [0,range) to which
X     * we will add offset to produce a number [offset,range+offset).
X     */
X
X    /*
X     * loop until we get a value that is in range
X     *
X     * We will obtain pseudo-random values over the range [0,2^rangebits)
X     * where 2^rangebits >= range and 2^(rangebits-1) < range.  We
X     * will ignore any results that are > the range that we want.
X     *
X     * A note in modulus biasing:
X     *
X     * We will not fall into the trap of thinking that we can simply take
X     * a value mod 'range'.  Consider the case where 'range' is '80'
X     * and we are given pseudo-random numbers [0,100).  If we took them
X     * mod 80, then the numbers [0,20) would be produced more often
X     * because the numbers [81,100) mod 80 wrap back into [0,20).
X     */
X    do {
X    /* obtain a pseudo-random value */
X    ret = cryrand(rangebits);
X    } while (ret >= range);
X
X    /*
X     * return the adjusted range value
X     */
X    return(ret+offset);
X}
X
X
X/*
X * srandom - seed the cryptographically strong pseudo-random number generator
X *
X * given:
X *    seed    a random number seed
X *
X * returns:
X *      the previous seed
X *
X * This function is just another interface to the crypto generator.
X * (see the scryrand() function).
X *
X * This function makes indirect use of the additive 55 and shuffle
X * generator.
X *
X * There is no limit on the size of a seed.  On the other hand,
X * extremely large seeds require large tables and long seed times.
X * Using a seed in the range of [2^64, 2^64 * 128!) should be
X * sufficient for most purposes.  An easy way to stay within this
X * range to to use seeds that are between 21 and 215 digits long, or
X * 64 to 780 bits long.
X *
X * NOTE: Calling this function will clear any internally buffer bits.
X *     See cryrand() for details.
X *
X * NOTE: This function seeds the Blum cryptographically strong
X *     pseudo-random number generator!
X */
Xdefine
Xsrandom(seed)
X{
X    if (!isint(seed)) {
X    quit "bad arg: seed must be an integer";
X    }
X    if (seed < 0) {
X    quit "bad arg: seed < 0 is reserved for future use";
X    }
X    return(scryrand(seed));
X}
X
X
X/*
X * randstate - set/get the state of all of the generators
X *
X * usage:
X *    randstate()    return the current state
X *    randstate(0)    return the previous state, set the default state
X *    randstate(cobj)    return the previous state, set a new state
X *
X * In the first form: randstate()
X *
X *    This function returns an cryobj object containing information
X *    about the current state of all of the generators.
X *
X * In the second form: randstate(0)
X *
X *    This function sets all of the generators to the default initial
X *    state (i.e., the state when this library was loaded).
X *
X *    This function returns an cryobj object containing information
X *    about the previous state of all of the generators.
X *
X * In the third form: randstate(cobj)
X *
X *    This function sets all of the generators to the state as found
X *    in the cryobj object.
X *
X *    This function returns an cryobj object containing information
X *    about the previous state of all of the generators.
X *
X * This function may be used to save and restore cryrand() & random()
X * generator states.  For example:
X *
X *    state = randstate()        <-- save the current state
X *    random()            <-- print the next 64 bits
X *    randstate(state)        <-- restore previous state
X *    random()            <-- print the same 64 bits
X *
X * One may quickly reseed a generator.  For example:
X *
X *    srandom(1,330,350)        <-- seed the generator
X *    seed1state = randstate()    <-- remember this 1st seeded state
X *    random()            <-- print 1st 64 bits seed 1 generator
X *    srandom(2,331,351)        <-- seed the generator again
X *    seed2state = randstate()    <-- remember this 2nd seeded state
X *    random()            <-- print 1st 64 bits seed 2 generator
X *
X *    randstate(seed1state)        <-- reseed to the 1st seeded state
X *    random()            <-- reprint 1st 64 bits seed 1 generator
X *    randstate(seed2state)        <-- reseed to the 2nd seeded state
X *    random()            <-- reprint 1st 64 bits seed 2 generator
X *
X *    oldstate = randstate(0)        <-- seed to the default generator
X *    random()            <-- print 1st 64 bits from default
X *    a55rand()            <-- print 1st 64 bits a55 generator
X *    prevstate = randstate(oldstate)    <-- restore seed 2 generator
X *    random()            <-- print 2nd 64 bits seed 2 generator
X *    randstate(prevstate)        <-- restore def generator in progress
X *    random()            <-- print 2nd 64 bits default generator
X *    a55rand()            <-- print 2nd 64 bits a55 generator
X *
X * given:
X *    cobj    if a cryobj object, use that object to set the current state
X *        if 0, set to the default state
X *
X * return:
X *    return the state of the crypto pseudo-random number generator in
X *           the form of an cryobj object, as it was prior to this call
X *
X * NOTE: No checking is performed on the data the 3rd form (cryobj object
X *     arg) is used.  The user must ensure that the arg represents a valid
X *     generator state.
X *
X * NOTE: When using the second form (passing an integer arg), only 0 is
X *     defined.  All other integer values are reserved for future use.
X */
Xdefine
Xrandstate(arg)
X{
X    /* declare our objects */
X    local obj cryobj x;        /* firewall comparator */
X    local obj cryobj prev;    /* previous states of the generators */
X    local junk;            /* dummy holder of random values */
X
X    /* firewall */
X    if (!isint(arg) && !istype(arg,x) && !isnull(arg)) {
X    quit "bad arg: argument must be integer, an cryobj object or missing";
X    }
X    if (isint(arg) && arg != 0) {
X        quit "bad arg:  non-zero integer arguments are reserved for future use";
X    }
X
X    /*
X     * save the current state
X     */
X    prev.p = cryrand_p;
X    prev.q = cryrand_q;
X    prev.r = cryrand_r;
X    prev.exp = cryrand_exp;
X    prev.left = cryrand_left;
X    prev.bitcnt = cryrand_bitcnt;
X    prev.a55j = add55_j;
X    prev.a55k = add55_k;
X    prev.seed = cry_seed;
X    prev.shufy = shuf_y;
X    prev.shufsz = shuf_size;
X    prev.shuftbl = shuf_tbl;
X    prev.a55tbl = add55_tbl;
X    if (isnull(x)) {
X    /* if no args, just return current state */
X    return (prev);
X    }
X
X    /*
X     * deal with the cryobj arg - set the state
X     */
X    if (istype(arg, x)) {
X    /* set the state from this object */
X    cryrand_p = arg.p;
X    cryrand_q = arg.q;
X    cryrand_n = cryrand_p*cryrand_q;
X    cryrand_r = arg.r;
X    cryrand_exp = arg.exp;
X    cryrand_left = arg.left;
X    cryrand_bitcnt = arg.bitcnt;
X    add55_j = arg.a55j;
X    add55_k = arg.a55k;
X    cry_seed = arg.seed;
X    add55_seed64 = randreseed64(cry_seed);
X    shuf_y = arg.shufy;
X    shuf_size = arg.shufsz;
X    shuf_shift = (64-highbit(shuf_size));
X    mat shuf_tbl[shuf_size];
X    shuf_tbl = arg.shuftbl;
X    add55_tbl = arg.a55tbl;
X
X    /*
X     * deal with the 0 integer arg - set the default initial state
X     */
X    } else if (isint(arg) && arg == 0) {
X    cryrand_p = cryrand_init_p;
X    cryrand_q = cryrand_init_q;
X    cryrand_n = cryrand_p * cryrand_q;
X    cryrand_r = cryrand_init_r;
X    cryrand_exp = cryrand_r;
X    cryrand_left = 0;
X    cryrand_bitcnt = 0;
X    cry_seed = 0;
X    cry_seed = sshufrand(0);
X    }
X
X    /*
X     * return the previous state
X     */
X    return (prev);
X}
X
X
X/*
X * nxtprime - find a probable prime >= n_arg
X *
X * usage:
X *    nxtprime(n_arg)
X *    nxtprime(n_arg, modval, modulus)
X *
X * given:
X *    n_arg            lower bound of the search
X *    [modval modulus]    if given, look for numbers mod modulus == modval
X *
X * returns:
X *    A number is that is very likely prime.
X *
X * In the first case 'nxtprime(n_arg)', this function returns a probable
X * prime >= n_arg.  In the second case 'nxtprime(n_arg, v, u)', this
X * function returns a probable prime >= n_arg that is also == v mod u.
X *
X * This function will not skip over a prime, through there is a
X * extremely unlikely chance that it will return a composite.
X * The odds that a number returned by this function is not prime
X * are 1 in 4^50.  The failure rate of this function is many orders
X * or magnitude lower than the failure rate due to a hardware error.
X *
X * NOTE: This function can take a long time, given a large value of n_arg.
X */
Xdefine
Xnxtprime(n_arg, modval, modulus)
X{
X    local modgcd;        /* gcd(modulus,modval) */
X    local n;            /* value >= n_arg that is being tested */
X    local j;
X
X    /*
X     * firewall
X     */
X    if (!isint(n_arg)) {
X    quit "bad args: 1st arg must be an integer";
X    }
X    if (!isnull(modulus) && !isint(modval)) {
X    quit "bad args: 3rd arg, if 2nd arg is given, must be an integer";
X    }
X    if (!isnull(modulus) && (!isint(modulus) || modulus <= 0)) {
X    quit "bad args: 3nd arg, if given, must be an integer > 0";
X    }
X
X    /*
X     * get values < 3 out of the way
X     */
X    n = n_arg;
X    if (n < 3) {
X    /* get the even prime out of the way, if possible */
X    if (isnull(modulus) ||
X        modulus == 1 ||
X        (n%modulus == modval%modulus)) {
X        /*
X         * 2 is the greatest odd prime, because
X         * 2 is the least even prime  :-)
X         */
X        return(2);
X    }
X    /* we have eliminated everything < 3 */
X    n = 3;
X    }
X
X    /*
X     * convert nxtprime(n) to nxtprime(n,1,2)
X     * convert nxtprime(n,x,1) to nxtprime(n,1,2)
X     * convert nxtprime(n,a,b) to nxtprime(n,a mod b,b)
X     */
X    if (isnull(modulus) || modulus < 2) {
X    modulus = 2;
X    modval = 1;
X    }
X    modval %= modulus;
X
X    /*
X     * catch cases where no more primes == 'modval' mod 'modulus' exist
X     */
X    modgcd = gcd(modval,modulus);
X    if (modgcd > 1) {
X
X    /* if beyond the modgcd, then no primes can exist */
X        if (n > modgcd) {
X        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
X        quit "no such prime of that form exists";
X    }
X
X    /* else n <= modgcd, then our only chance is if modgcd is prime */
X    /* reject if modgcd has an obvious prime factor */
X    if (modgcd > 10 && gcd(modgcd,nxtprime_pfact10) != 1) {
X        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
X        quit "no such prime of that form exists";
X    }
X    if (modgcd > 100 && gcd(modgcd,nxtprime_pfact100) != 1) {
X        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
X        quit "no such prime of that form exists";
X    }
X    if (modgcd > 1000 && gcd(modgcd,nxtprime_pfact1000) != 1) {
X        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
X        quit "no such prime of that form exists";
X    }
X
X    /* do 50 probable prime tests, for a 1 in 4^50 false prime rate */
X    if (!ptest(modgcd,50)) {
X        print "n_arg:",n_arg,"  modval:",modval,"  modulus:",modulus;
X        quit "no such prime of that form exists";
X    }
X
X    /* modgcd is the only prime >= n */
X    return(modgcd);
X    }
X
X    /*
X     * bump n up to the next possible case
X     *
X     * n will be an odd number == 'modval' mod 'modulus'
X     */
X    if (n%modulus != modval) {
X    j = n - (n%modulus) + modval;
X    if (j < n) {
X        n = j+modulus;
X    } else {
X        n = j;
X    }
X    }
X    if (n%2 == 0) {
X    n += modulus;
X    }
X
X    /* look for a prime */
X    n = n-modulus;
X    do {
X    /* try the next integer */
X    n = n+modulus;
X
X    /* reject if it has an obvious prime factor */
X    if (n > 10 && gcd(n,nxtprime_pfact10) != 1) {
X        continue;
X    }
X    if (n > 100 && gcd(n,nxtprime_pfact100) != 1) {
X        continue;
X    }
X    if (n > 1000 && gcd(n,nxtprime_pfact1000) != 1) {
X        continue;
X    }
X
X    /* do 50 probable prime tests */
X    if (!ptest(n,50)) {
X        continue;
X    }
X
X    /* n is very likely a prime number */
X    return(n);
X
X    } while (1);
X}
X
X
X/*
X * cryobj - how to initialize a cryobj object
X *
X * given:
X *    p        first Blum prime (prime 3 mod 4)
X *    q        second Blum prime (prime 3 mod 4)
X *    r        quadratic residue of n=p*q
X *    exp        used in computing crypto good bits
X *    left        bits unused from the last cryrand() call
X *    bitcnt        left contains bitcnt crypto good bits
X *    a55j        1st additive 55 table pointer
X *    a55k        2nd additive 55 table pointer
X *    seed        last seed set by sa55rand() or 0
X *    shufy        Y (previous a55rand() output for shuffle)
X *    shufsz        size of the shuffle table
X *    shuftbl        a matrix of shufsz entries
X *    a55tbl        additive 55 generator state table
X *
X * return:
X *    an cryobj object
X *
X * NOTE: This function, by convention, returns an cryobj object.
X */
Xdefine
Xcryobj(p,q,r,exp,left,bitcnt,a55j,a55k,seed,shufy,shufsz,shuftbl,a55tbl)
X{
X    /* declare our objects */
X    local obj cryobj x;
X
X    /* firewall */
X    if (!isint(p) || !isint(q) || !isint(r) || !isint(exp) || \
X    !isint(left) || !isint(bitcnt) || !isint(a55j) || \
X    !isint(a55k) || !isint(seed) || !isint(shufy) || !isint(shufsz)) {
X    quit "bad args: first 11 args must be integers";
X    }
X    if (!ismat(shuftbl) || matdim(shuftbl) != 1 || \
X    matmin(shuftbl,1) != 0 || matmax(shuftbl,1) != shuf_size-1) {
X    quit "bad arg: 12th is not a mat[0:shuf_size-1]";
X    }
X    if (!ismat(a55tbl) || matdim(a55tbl) != 1 || \
X    matmin(a55tbl,1) != 0 || matmax(a55tbl,1) != 54) {
X    quit "bad arg: 13th is not a mat[0:54]";
X    }
X
X    /* initialize object with default startup values */
X    x.p = p;
X    x.q = q;
X    x.r = r;
X    x.exp = exp;
X    x.left = left;
X    x.bitcnt = bitcnt;
X    x.a55j = a55j;
X    x.a55k = a55k;
X    x.seed = seed;
X    x.shufy = shuf_y;
X    x.shufsz = shuf_size;
X    x.shuftbl = shuf_tbl;
X    x.a55tbl = a55tbl;
X
X    /* return the initialized object */
X    return (x);
X}
X
X
X/*
X * cryobj_print - print the value of a cryobj object
X *
X * usage:
X *    a    an cryobj object
X *
X * NOTE: This function is called automatically when an cryobj object
X *       is displayed.
X */
Xdefine
Xcryobj_print(a)
X{
X    /* declare our objects */
X    local obj cryobj x;        /* firewall comparator */
X
X    /* firewall */
X    if (!istype(a, x)) {
X    quit "bad arg: arg is not an cryobj object";
X    }
X    if (!ismat(a.shuftbl) || matdim(a.shuftbl) != 1 || \
X    matmin(a.shuftbl,1) != 0 || matmax(a.shuftbl,1) != a.shufsz-1) {
X    quit "bad arg: 12th is not a mat[0:shuf_size-1]";
X    }
X    if (!ismat(a.a55tbl) || matdim(a.a55tbl) != 1 || \
X    matmin(a.a55tbl,1) != 0 || matmax(a.a55tbl,1) != 54) {
X    quit "bad arg: 13th is not a mat[0:54]";
X    }
X
X    /* print the value */
X    print "cryobj(" : a.p : "," : a.q : "," : a.r : "," : a.exp : "," : \
X      a.left : "," : a.bitcnt : "," : a.a55j : "," : a.a55k : "," : \
X      a.seed : "," : a.shufy : "," : a.shufsz : \
X      ",[" : a.shuftbl[0] : "," : a.shuftbl[1] : "," : \
X      a.shuftbl[2] : ",...," : a.shuftbl[52] : "," : \
X      a.shuftbl[53] : "," : a.shuftbl[54] : "]" : \
X      ",[" : a.a55tbl[0] : "," : a.a55tbl[1] : "," : \
X      a.a55tbl[2] : ",...," : a.a55tbl[52] : "," : \
X      a.a55tbl[53] : "," : a.a55tbl[54] : "])" : ;
X}
X
X
X/*
X * cryres - find a pseudo-random quadratic residue for scryrand() and cryrand()
X *
X * given:
X *    p    first prime
X *    q    second prime
X *
X * returns:
X *    a number that is a quadratic residue of n=p*q
X *
X * This function is returns the pseudo-random quadratic residue of
X * the product of two primes.  Normally this function is called
X * only by the crypto generator.
X *
X * NOTE: No check is made to ensure that the values passed are primes.
X *
X * NOTE: This is NOT Blum's method.  This is used by Blum's method to
X *       generate some internal values.
X */
Xdefine
Xcryres(p,q)
X{
X    local rval;        /* a temporary pseudo-random value */
X    local sqrpow;    /* a power of 2 starting > n^(3/4) */
X    local quadres;    /* 0 or a quadratic residue of n = p*q */
X    local n;        /* n=p*q */
X    local j;
X
X    /*
X     * firewall
X     */
X    if (!isint(p) || !isint(q) || p < 3 || q < 3) {
X    quit "bad args: p and q must be integers > 2";
X    }
X
X    /*
X     * find a pseudo-random quadratic residue of n = p*q
X     *
X     * To find a pseudo-random quadratic residue, we will generate
X     * pseudo-random numbers to try in the interval [2^sqrpow,n-3),
X     * where 2^sqrpow is the first power of 2 > n^(3/4).  This range
X     * helps us avoid selecting trivial residues abs(quadres mod n=p*q)
X     * is tiny.  We will never select a residue < 4.
X     *
X     * When we fail to find a quadratic residue after 1009 tries, we will
X     * drop sqrpow by 1 and start at another pseudo-random location.
X     *
X     * It is very unlikely that we will need to search more than a
X     * few numbers to find a quadratic residue of n = p*q.
X     *
X     * We will reject any quadratic residue that is a square of
X     * itself, mod n=p*q.
X     */
X    n = p*q;
X    quadres = 0;
X    sqrpow = highbit(n)//4*3;
X    do {
X    /* generate a pseudo-random starting point */
X    rval = rand(2^sqrpow, n-4);
X
X    /* look for 1009 times or until >= cryrand_n */
X    for (j=0; j < 1009; ++j, ++rval) {
X
X        /*
X         * check if we have a quadratic residue of n = p*q
X         *
X         * p is prime and J(r,p) == 1  ==>  r is a quadratic residue of p
X         * q is prime and J(q,p) == 1  ==>  r is a quadratic residue of q
X         *
X         * J(r,p)*J(r,q) == J(r,p*q)
X         * J(r,p) and J(q,p) == 1       ==>  J(r,p*q) == 1
X         * J(r,p*q) == 1           ==>  r is a quadratic residue of p*q
X         */
X        if (jacobi(rval,p) == 1 && jacobi(rval,q) == 1) {
X        quadres = rval;
X        break;
X        }
X    }
X
X    /* setup for a new pseudo-random starting point if we found nothing */
X    if (quadres == 0) {
X        /* reduce sqrpow if reasonable */
X        if (sqrpow > 2) {
X        --sqrpow;
X        }
X    }
X    } while (quadres == 0 || ((quadres^2)%n) == quadres);
X
X    /*
X     * return the quadratic residue of n = p*q
X     */
X    return (quadres);
X}
X
X
X/*
X * randreseed64 - scramble a 64 bit seed
X *
X * given:
X *    a 64 bit seed
X *
X * returns:
X *    a 64 scrambled seed, or 0 if seed was 0
X *
X * It is 'nice' when a seed of "n" produces a 'significantly different'
X * sequence than a seed of "n+1".  Generators, by convention, assign
X * special significance to the seed of '0'.  It is an unfortunate that
X * people often pick small seed values, particularly when large seed
X * are of significance to the generators found in this file.
X *
X * If 'seed' is 0, then 0 is returned.  If 'seed' is non-zero, we will
X * produce a different and unique new scrambled 'seed'.  This scrambling
X * will effectively eliminate the human factors and perceptions that
X * are noted above.
X *
X * It should be noted that the purpose of this process to scramble a seed
X * ONLY.  We do not care if these generators produce good random numbers.
X * We only want to help eliminate the human factors and perceptions
X * noted above.
X *
X * This function scrambles the low 64 bits of a seed, by mapping [0,2^64)
X * into [0,2^64).  This map is one-to-one and onto.  Mapping is performed
X * using  a linear congruence generator of the form:
X *
X *        X1 <-- (a*X0 + c) mod m
X *
X * The generator are based on the linear congruential generators found in
X * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
X * vol 2, 2nd edition (1981), Section 3.6, pages 170-171.
X *
X * Because we process 64 bits we will take:
X *
X *        m = 2^64            (based on note ii)
X *
X * We will scan the Rand book of numbers to look for an 'a' such that:
X *
X *        a mod 8 == 5            (based on note iii)
X *        0.01*m < a < 0.99*m        (based on note iv)
X *        0.01*2^64 < a < 0.99*2^64
X *
X * To help keep the generators independent, we want:
X *
X *        a is prime
X *
X * The choice of an adder 'c' is considered immaterial according (based
X * in note v).  Knuth suggests 'c==1' or 'c==a'.  We elect to select 'c'
X * using the same process as we used to select 'a'.  The choice is
X * 'immaterial' after all, and as long as:
X *
X *        gcd(c, m) == 1        (based on note v)
X *        gcd(c, 2^64) == 1
X *
X * the concerns are met.   It can be shown that if we have:
X *
X *        gcd(a, c) == 1
X *
X * then the adders and multipliers will be more independent.
X *
X * We will obtain the values 'a' and 'c for our generator from the
X * Rand book of numbers.  Because m=2^64 is 20 decimal digits long, we
X * will search the Rand book of numbers 20 at a time.  We will skip any
X * of the 55 values that were used to initialize the additive 55 generators.
X * The values obtained from the Rand book are:
X *
X *        a = 6316878969928993981
X *        c = 1363042948800878693
X *
X * As we stated before, we must map 0 ==> 0.  The linear congruence
X * generator would normally map as follows:
X *
X *    0 ==> 1363042948800878693    (0 ==> c)
X *
X * To determine which value maps back into 0, we compute:
X *
X *    (-c*minv(a,m)) % m
X *
X * and thus we find that the congruence generator would also normally map:
X *
X *    10239951819489363767 ==> 0
X *
X * To overcome this, and preserve the 1-to-1 and onto map, we force:
X *
X *    0 ==> 0
X *    10239951819489363767 ==> 1363042948800878693
X *
X * To repeat, this function converts a values into a seed value.  With the
X * except of 'seed == 0', every value is mapped into a unique seed value.
X * This mapping need not be complex, random or secure.  All we attempt
X * to do here is to allow humans who pick small or successive seed values
X * to obtain reasonably different sequences from the generators below.
X *
X * NOTE: This is NOT a random number generator.  This function is intended
X *     to be used internally by sa55rand() and sshufrand().
X */
Xdefine
Xrandreseed64(seed)
X{
X    local ret;            /* return value */
X    local a;            /* generator 0 multiplier */
X    local c;            /* generator 0 adder */
X
X    /* firewall */
X    if (!isint(seed)) {
X    quit "bad args: seed must be an integer";
X    }
X    if (seed < 0) {
X    quit "bad arg: seed < 0 is reserved for future use";
X    }
X
X    /* if seed is 0, we will return 0 */
X    if (seed == 0) {
X    return (0);
X    }
X
X    /*
X     * process the low 64 bits of the seed
X     */
X    a = 6316878969928993981;
X    c = 1363042948800878693;
X    ret = (((seed & cry_mask)*a + c) & cry_mask);
X
X    /*
X     * Normally, the above equation would map:
X     *
X     *        f(0) == 1363042948800878693
X     *        f(10239951819489363767) == 0
X     *
X     * However, we have already forced f(0) == 0.  To preserve the
X     * 1-to-1 and onto map property, we force:
X     *
X     *        f(10239951819489363767) ==> 1363042948800878693
X     */
X    if (ret == 0) {
X    ret = c;
X    }
X
X    /* return the scrambled value */
X    return (ret);
X}
X
X
X/*
X * Initial read execution code
X */
Xcry_seed = srand(0);        /* pre-initialize the tables */
X
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X    print "ver: 10.5 06:19:34 5/9/93";
X    print "shufrand() defined";
X    print "sshufrand(seed) defined";
X    print "rand([a, [b]]) defined";
X    print "srand(seed) defined";
X    print "cryrand([a, [b]]) defined";
X    print "scryrand([seed, [len1, len2]]) defined";
X    print "scryrand(seed, ip, iq, ir) defined";
X    print "random([a, [b]]) defined";
X    print "srandom(seed) defined";
X    print "obj cryobj defined";
X    print "randstate([cryobj | 0]) defined";
X    print "nxtprime(n, [val, modulus]) defined";
X}
SHAR_EOF
echo "File calc2.9.0/lib/cryrand.cal is complete"
chmod 0644 calc2.9.0/lib/cryrand.cal || echo "restore of calc2.9.0/lib/cryrand.cal fails"
set `wc -c calc2.9.0/lib/cryrand.cal`;Sum=$1
if test "$Sum" != "82965"
then echo original size 82965, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/deg.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/deg.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Calculate in degrees, minutes, and seconds.
X */
X
Xobj dms {deg, min, sec};
X
Xdefine dms(deg, min, sec)
X{
X    local ans;
X
X    if (isnull(sec))
X        sec = 0;
X    if (isnull(min))
X        min = 0;
X    obj dms ans;
X    ans.deg = deg;
X    ans.min = min;
X    ans.sec = sec;
X    fixdms(&ans);
X    return ans;
X}
X
X
Xdefine dms_add(a, b)
X{
X    local obj dms    ans;
X
X    ans.deg = 0;
X    ans.min = 0;
X    ans.sec = 0;
X    if (istype(a, ans)) {
X        ans.deg += a.deg;
X        ans.min += a.min;
X        ans.sec += a.sec;
X    } else
X        ans.deg += a;
X    if (istype(b, ans)) {
X        ans.deg += b.deg;
X        ans.min += b.min;
X        ans.sec += b.sec;
X    } else
X        ans.deg += b;
X    fixdms(&ans);
X    return ans;
X}
X
X
Xdefine dms_neg(a)
X{
X    local obj dms    ans;
X
X    ans.deg = -ans.deg;
X    ans.min = -ans.min;
X    ans.sec = -ans.sec;
X    return ans;
X}
X
X
Xdefine dms_sub(a, b)
X{
X    return a - b;
X}
X
X
Xdefine dms_mul(a, b)
X{
X    local obj dms    ans;
X
X    if (istype(a, ans) && istype(b, ans))
X        quit "Cannot multiply degrees together";
X    if (istype(a, ans)) {
X        ans.deg = a.deg * b;
X        ans.min = a.min * b;
X        ans.sec = a.sec * b;
X    } else {
X        ans.deg = b.deg * a;
X        ans.min = b.min * a;
X        ans.sec = b.sec * a;
X    }
X    fixdms(&ans);
X    return ans;
X}
X
X
Xdefine dms_print(a)
X{
X    print a.deg : 'd' : a.min : 'm' : a.sec : 's' :;
X}
X
X
Xdefine dms_abs(a)
X{
X    return a.deg + a.min / 60 + a.sec / 3600;
X}
X
X
Xdefine fixdms(a)
X{
X    a.min += frac(a.deg) * 60;
X    a.deg = int(a.deg);
X    a.sec += frac(a.min) * 60;
X    a.min = int(a.min);
X    a.min += a.sec // 60;
X    a.sec %= 60;
X    a.deg += a.min // 60;
X    a.min %= 60;
X    a.deg %= 360;
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X    print "obj dms {deg, min, sec} defined";
X    print "dms(deg, min, sec) defined";
X    print "dms_add(a, b) defined";
X    print "dms_neg(a) defined";
X    print "dms_sub(a, b) defined";
X    print "dms_mul(a, b) defined";
X    print "dms_print(a) defined";
X    print "dms_abs(a) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/deg.cal || echo "restore of calc2.9.0/lib/deg.cal fails"
set `wc -c calc2.9.0/lib/deg.cal`;Sum=$1
if test "$Sum" != "1946"
then echo original size 1946, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/ellip.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/ellip.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Attempt to factor numbers using elliptic functions.
X *     y^2 = x^3 + a*x + b   (mod N).
X *
X * Many points (x,y) (mod N) are found that solve the above equation,
X * starting from a trivial solution and 'multiplying' that point together
X * to generate high powers of the point, looking for such a point whose
X * order contains a common factor with N.  The order of the group of points
X * varies almost randomly within a certain interval for each choice of a
X * and b, and thus each choice provides an independent opportunity to
X * factor N.  To generate a trivial solution, a is chosen and then b is
X * selected so that (1,1) is a solution.  The multiplication is done using
X * the basic fact that the equation is a cubic, and so if a line hits the
X * curve in two rational points, then the third intersection point must
X * also be rational.  Thus by drawing lines between known rational points
X * the number of rational solutions can be made very large.  When modular
X * arithmetic is used, solving for the third point requires the taking of a
X * modular inverse (instead of division), and if this fails, then the GCD
X * of the failing value and N provides a factor of N.  This description is
X * only an approximation, read "A Course in Number Theory and Cryptography"
X * by Neal Koblitz for a good explanation.
X *
X * factor(iN, ia, B, force)
X *    iN is the number to be factored.
X *    ia is the initial value of a in the equation, and each successive
X *    value of a is an independent attempt at factoring (default 1).
X *    B is the limit of the primes that make up the high power that the
X *    point is raised to for each factoring attempt (default 100).
X *    force is a flag to attempt to factor numbers even if they are
X *    thought to already be prime (default FALSE).
X *
X * Making B larger makes the power the point being raised to contain more
X * prime factors, thus increasing the chance that the order of the point
X * will be made up of those factors.  The higher B is then, the greater
X * the chance that any individual attempt will find a factor.  However,
X * a higher B also slows down the number of independent functions being
X * examined.  The order of the point for any particular function might
X * contain a large prime and so won't succeed even for a really large B,
X * whereas the next function might have an order which is quickly found.
X * So you want to trade off the depth of a particular search with the
X * number of searches made.  For example, for factoring 30 digits, I make
X * B be about 1000 (probably still too small).
X *
X * If you have lots of machines available, then you can run parallel
X * factoring attempts for the same number by giving different starting
X * values of ia for each machine (e.g. 1000, 2000, 3000).
X *
X * The output as the function is running is (occasionally) the value of a
X * when a new function is started, the prime that is being included in the
X * high power being calculated, and the current point which is the result
X * of the powers so far.
X *
X * If a factor is found, it is returned and is also saved in the global
X * variable f.  The number being factored is also saved in the global
X * variable N.
X */
X
Xobj point {x, y};
Xglobal    N;        /* number to factor */
Xglobal    a;        /* first coefficient */
Xglobal    b;        /* second coefficient */
Xglobal    f;        /* found factor */
X
X
Xdefine factor(iN, ia, B, force)
X{
X    local    C, x, p;
X
X    if (!force && ptest(iN, 50))
X        return 1;
X    if (isnull(B))
X        B = 100;
X    if (isnull(ia))
X        ia = 1;
X    obj point x;
X    a = ia;
X    b = -ia;
X    N = iN;
X    C = isqrt(N);
X    C = 2 * C + 2 * isqrt(C) + 1;
X    f = 0;
X    while (f == 0) {
X        print "A =", a;
X        x.x = 1;
X        x.y = 1;
X        print 2, x;
X        x = x ^ (2 ^ (highbit(C) + 1));
X        for (p = 3; ((p < B) && (f == 0)); p += 2) {
X            if (!ptest(p, 1))
X                continue;
X            print p, x;
X            x = x ^ (p ^ ((highbit(C) // highbit(p)) + 1));
X        }
X        a++;
X        b--;
X    }
X    return f;
X}
X
X
Xdefine point_print(p)
X{
X    print "(" : p.x : "," : p.y : ")" :;
X}
X
X
Xdefine point_mul(p1, p2)
X{
X    local    r, m;
X
X    if (p2 == 1)
X        return p1;
X    if (p1 == p2)
X        return point_square(&p1);
X    obj point r;
X    m = (minv(p2.x - p1.x, N) * (p2.y - p1.y)) % N;
X    if (m == 0) {
X        if (f == 0)
X            f = gcd(p2.x - p1.x, N);
X        r.x = 1;
X        r.y = 1;
X        return r;
X    }
X    r.x = (m^2 - p1.x - p2.x) % N;
X    r.y = ((m * (p1.x - r.x)) - p1.y) % N;
X    return r;
X}
X
X
Xdefine point_square(p)
X{
X    local    r, m;
X
X    obj point r;
X    m = ((3 * p.x^2 + a) * minv(p.y << 1, N)) % N;
X    if (m == 0) {
X        if (f == 0)
X            f = gcd(p.y << 1, N);
X        r.x = 1;
X        r.y = 1;
X        return r;
X    }
X    r.x = (m^2 - p.x - p.x) % N;
X    r.y = ((m * (p.x - r.x)) - p.y) % N;
X    return r;
X}
X
X
Xdefine point_pow(p, pow)
X{
X    local bit, r, t;
X
X    r = 1;
X    if (isodd(pow))
X        r = p;
X    t = p;
X    for (bit = 2; ((bit <= pow) && (f == 0)); bit <<= 1) {
X        t = point_square(&t);
X        if (bit & pow)
X            r = point_mul(&t, &r);
X    }
X    return r;
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X    print "factor(N, I, B, force) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/ellip.cal || echo "restore of calc2.9.0/lib/ellip.cal fails"
set `wc -c calc2.9.0/lib/ellip.cal`;Sum=$1
if test "$Sum" != "5017"
then echo original size 5017, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/lucas.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/lucas.cal &&
X/*
X * Copyright (c) 1993 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}!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 * At the time of discovery, this number was the largest known prime.
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 28 Aug 1993, this prime was the 2nd
X * largest known prime and the largest known non-Mersenne 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 28 Aug 1993, 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 Amdahl 6 team used a 'ultra-high speed
X * multi-precision' package that a machine dependent collection of routines
X * tuned for a long trace vector processor to work with very large numbers.
X * The heart of the package was a multiplication and square routine that
X * was based on the PFA Fast Fourier Transform and on Winograd's radix FFTs.
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 divisibility
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 divisible by a small primes.  The ideas behind this
X * point is beyond the scope of this program.
X */
X
Xglobal pprod256;    /* product of  "primes up to 256" / "primes up to 46" */
Xglobal lib_debug;    /* 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 divisible 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)) {
SHAR_EOF
echo "End of part 17"
echo "File calc2.9.0/lib/lucas.cal is continued in part 18"
echo "18" > s2_seq_.tmp
exit 0