Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/unix/volume27/calc-2.9.0/part16

Text File | 1993-12-07 | 60KB | 1,805 lines

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

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

#!/bin/sh
# this is part 16 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/lib/README continued
#
CurArch=16
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/README"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/lib/README
X
X    qpi(epsilon)
X
X    Calculate pi within the specified epsilon using the quartic convergence
X    iteration.
X
X
Xpollard.cal
X
X    factor(N, N, ai, af)
X
X    Factor using Pollard's p-1 method.
X
X
Xpoly.cal
X
X    Calculate with polynomials of one variable.  There are many functions.
X    Read the documentation in the library file.
X
X
Xpsqrt.cal
X
X    psqrt(u, p)
X
X    Calculate square roots modulo a prime
X
X
Xquat.cal
X
X    quat(a, b, c, d)
X    quat_print(a)
X    quat_norm(a)
X    quat_abs(a, e)
X    quat_conj(a)
X    quat_add(a, b)
X    quat_sub(a, b)
X    quat_inc(a)
X    quat_dec(a)
X    quat_neg(a)
X    quat_mul(a, b)
X    quat_div(a, b)
X    quat_inv(a)
X    quat_scale(a, b)
X    quat_shift(a, b)
X
X    Calculate using quaternions of the form: a + bi + cj + dk.  In these
X    functions, quaternians are manipulated in the form: s + v, where
X    s is a scalar and v is a vector of size 3.
X
X
Xregress.cal
X
X    Test the correct execution of the calculator by reading this library file.
X    Errors are reported with '****' mssages, or worse.  :-)
X
X
Xsolve.cal
X
X    solve(low, high, epsilon)
X
X    Solve the equation f(x) = 0 to within the desired error value for x.
X    The function 'f' must be defined outside of this routine, and the low
X    and high values are guesses which must produce values with opposite signs.
X
X
Xsumsq.cal
X
X    ss(p)
X
X    Determine the unique two positive integers whose squares sum to the
X    specified prime.  This is always possible for all primes of the form
X    4N+1, and always impossible for primes of the form 4N-1.
X
X
Xsurd.cal
X
X    surd(a, b)
X    surd_print(a)
X    surd_conj(a)
X    surd_norm(a)
X    surd_value(a, xepsilon)
X    surd_add(a, b)
X    surd_sub(a, b)
X    surd_inc(a)
X    surd_dec(a)
X    surd_neg(a)
X    surd_mul(a, b)
X    surd_square(a)
X    surd_scale(a, b)
X    surd_shift(a, b)
X    surd_div(a, b)
X    surd_inv(a)
X    surd_sgn(a)
X    surd_cmp(a, b)
X    surd_rel(a, b)
X
X    Calculate using quadratic surds of the form: a + b * sqrt(D).
X
X
Xunitfrac.cal
X
X    unitfrac(x)
X
X    Represent a fraction as sum of distinct unit fractions.
X
X
Xvarargs.cal
X
X    sc(a, b, ...)
X
X    Example program to use 'varargs'.  Program to sum the cubes of all
X    the specified numbers.
SHAR_EOF
echo "File calc2.9.0/lib/README is complete"
chmod 0644 calc2.9.0/lib/README || echo "restore of calc2.9.0/lib/README fails"
set `wc -c calc2.9.0/lib/README`;Sum=$1
if test "$Sum" != "5516"
then echo original size 5516, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/bernoulli.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/bernoulli.cal &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Calculate the Nth Bernoulli number B(n).
X * This uses the following symbolic formula to calculate B(n):
X *
X *    (b+1)^(n+1) - b^(n+1) = 0
X *
X * where b is a dummy value, and each power b^i gets replaced by B(i).
X * For example, for n = 3:
X *    (b+1)^4 - b^4 = 0
X *    b^4 + 4*b^3 + 6*b^2 + 4*b + 1 - b^4 = 0
X *    4*b^3 + 6*b^2 + 4*b + 1 = 0
X *    4*B(3) + 6*B(2) + 4*B(1) + 1 = 0
X *    B(3) = -(6*B(2) + 4*B(1) + 1) / 4
X *
X * The combinatorial factors in the expansion of the above formula are
X * calculated interatively, and we use the fact that B(2i+1) = 0 if i > 0.
X * Since all previous B(n)'s are needed to calculate a particular B(n), all
X * values obtained are saved in an array for ease in repeated calculations.
X */
Xstatic Bnmax;
Xstatic mat Bn[1001];
X
X
Xdefine B(n)
X{
X    local    nn, np1, i, sum, mulval, divval, combval;
X
X    if (!isint(n) || (n < 0))
X        quit "Non-negative integer required for Bernoulli";
X
X    if (n == 0)
X        return 1;
X    if (n == 1)
X        return -1/2;
X    if (isodd(n))
X        return 0;
X    if (n > 1000)
X        quit "Very large Bernoulli";
X
X    if (n <= Bnmax)
X        return Bn[n];
X
X    for (nn = Bnmax + 2; nn <= n; nn+=2) {
X        np1 = nn + 1;
X        mulval = np1;
X        divval = 1;
X        combval = 1;
X        sum = 1 - np1 / 2;
X        for (i = 2; i < np1; i+=2) {
X            combval = combval * mulval-- / divval++;
X            combval = combval * mulval-- / divval++;
X            sum += combval * Bn[i];
X        }
X        Bn[nn] = -sum / np1;
X    }
X    Bnmax = n;
X    return Bn[n];
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X    print "B(n) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/bernoulli.cal || echo "restore of calc2.9.0/lib/bernoulli.cal fails"
set `wc -c calc2.9.0/lib/bernoulli.cal`;Sum=$1
if test "$Sum" != "1492"
then echo original size 1492, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/bigprime.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/bigprime.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 * A prime test, base a, on p*2^x+1 for even x>m.
X */
X
Xdefine bigprime(a, m, p)
X{
X    local n1, n;
X
X    n1 = 2^m * p;
X    for (;;) {
X        m++;
X        n1 += n1;
X        n = n1 + 1;
X        if (isodd(m))
X            continue;
X        print m;
X        if (pmod(a, n1 / 2, n) != n1)
X            continue;
X        if (pmod(a, n1 / p, n) == 1)
X            continue;
X        print "    " : n;
X    }
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X    print "bigprime(a, m, p) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/bigprime.cal || echo "restore of calc2.9.0/lib/bigprime.cal fails"
set `wc -c calc2.9.0/lib/bigprime.cal`;Sum=$1
if test "$Sum" != "555"
then echo original size 555, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/bindings (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/bindings &&
X# Default key bindings for calc line editing functions
X
Xmap    base-map
Xdefault    insert-char
X^@    set-mark
X^A    start-of-line
X^B    backward-char
X^D    delete-char
X^E    end-of-line
X^F    forward-char
X^H    backward-kill-char
X^J    new-line
X^K    kill-line
X^L    refresh-line
X^M    new-line
X^N    forward-history
X^O    save-line
X^P    backward-history
X^R    reverse-search
X^T    swap-chars
X^U    flush-input
X^V    quote-char
X^W    kill-region
X^Y    yank
X^?    backward-kill-char
X^[    ignore-char    esc-map
X
Xmap    esc-map
Xdefault    ignore-char    base-map
XG    start-of-line
XH    backward-history
XP    forward-history
XK    backward-char
XM    forward-char
XO    end-of-line
XS    delete-char
Xg    goto-line
Xs    backward-word
Xt    forward-word
Xd    forward-kill-word
Xu    uppercase-word
Xl    lowercase-word
Xh    list-history
X^[    flush-input
X[    arrow-key
SHAR_EOF
chmod 0644 calc2.9.0/lib/bindings || echo "restore of calc2.9.0/lib/bindings fails"
set `wc -c calc2.9.0/lib/bindings`;Sum=$1
if test "$Sum" != "730"
then echo original size 730, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/chrem.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/chrem.cal &&
X/*
X * chrem - Chinese remainder theorem/problem solver
X *
X * When possible, chrem finds solutions for x of a set of congruences
X * of the form:
X *
X *            x = r1 (mod m1)
X *            x = r2 (mod m2)
X *               ...
X *
X * where the residues r1, r2, ... and the moduli m1, m2, ... are
X * given integers.  The Chinese remainder theorem states that if
X * m1, m2, ... are relatively prime in pairs, the above congruences
X * have a unique solution modulo  m1 * m2 * ...   If m1, m2, ...
X * are not relatively prime in pairs, it is possible that no solution
X * exists.  If solutions exist, the general solution is expressible as:
X *
X *                   x = r (mod m)
X *
X * where m = lcm(m1,m2,...), and if m > 0, 0 <= r < m.  This
X * solution may be interpreted as:
X *
X *                  x = r + k * m            [[NOTE 1]]
X *
X * where k is an arbitrary integer.
X *
X ***
X *
X * usage:
X *
X *    chrem(r1,m1 [,r2,m2, ...])
X *
X *         r1, r2, ...      remainder integers or null values
X *         m1, m2, ...      moduli integers
X *
X *    chrem(r_list, [m_list])
X *
X *        r_list    list (r1,r2, ...)
X *        m_list    list (m1,m2, ...)
X *
X *         If m_list is omitted, then 'defaultmlist' is used.
X *        This default list is a global value that may be changed
X *        by the user.  Initially it is the first 8 primes.
X *
X * If a remainder is null(), then the corresponding congruence is
X * ignored.  This is useful when working with a fixed list of moduli.
X *
X * If there are more remainders than moduli, then the later moduli are
X * ignored.
X *
X * The moduli may be any integers, not necessarily relatively prime in
X * pairs (as required for the Chinese remainder theorem).  Any moduli
X * may be zero;  x = r (mod 0) has the meaning of x = r.
X *
X * returns:
X *
X *    If args were integer pairs:
X *
X *          r        ('r' is defined above, see [[NOTE 1]])
X *
X *    If 1 or 2 list args were given:
X *
X *      (r, m)    ('r' and 'm' are defined above, see [[NOTE 1]])
X *
X * NOTE: In all cases, null() is returned if there is no solution.
X *
X ***
X *
X * This function may be used to solve the following historical problems:
X *
X *   Sun-Tsu, 1st century A.D.
X *
X *    To find a number for which the reminders after division by 3, 5, 7
X *    are 2, 3, 2, respectively:
X *
X *        chrem(2,3,3,5,2,7) ---> 23
X *
X *   Fibonacci, 13th century A.D.
X *
X *    To find a number divisible by 7 which leaves remainder 1 when
X *    divided by 2, 3, 4, 5, or 6:
X *
X *
X *        chrem(list(0,1,1,1,1,1),list(7,2,3,4,5,6)) ---> (301,420)
X *
X *    i.e., any value that is 301 mod 420.
X *
X * Written by: Ernest W Bowen <[e-mail]>
X * Interface by: Landon Curt Noll <[e-mail]>
X */
X
Xstatic defaultmlist = list(2,3,5,7,11,13,17,19); /* The first eight primes */
X
Xdefine chrem()
X{
X    local argc;        /* number of args given */
X    local rlist;    /* reminder list - ri */
X    local mlist;    /* modulus list - mi */
X    local list_args;      /* true => args given are lists, not r1,m1, ... */
X    local m,z,r,y,d,t,x,u,i;
X
X    /*
X     * parse args
X     */
X    argc = param(0);
X    if (argc == 0) {
X    quit "usage: chrem(r1,m1 [,r2,m2 ...]) or chrem(r_list, m_list)";
X    }
X    list_args = islist(param(1));
X    if (list_args) {
X    rlist = param(1);
X    mlist = (argc == 1) ? defaultmlist : param(2);
X    if (size(rlist) > size(mlist)) {
X        quit "too many residues";
X    }
X    } else {
X    if (argc % 2 == 1) {
X        quit "odd number integers given";
X    }
X    rlist = list();
X    mlist = list();
X    for (i=1; i <= argc; i+=2) {
X        push(rlist, param(i));
X            push(mlist, param(i+1));
X    }
X    }
X
X    /*
X     * solve the problem found in rlist & mlist
X     */
X    m = 1;
X    z = 0;
X    while (size(rlist)) {
X    r=pop(rlist);
X    y=abs(pop(mlist));
X    if (r==null())
X        continue;
X    if (m) {
X        if (y) {
X        d = t = z - r;
X        m = lcm(x=m, y);
X        while (d % y) {
X            u = x;
X            x %= y;
X            swap(x,y);
X            if (y==0)
X            return;
X            z += (t *= -u/y);
X        }
X        } else {
X        if ((r % m) != (z % m))
X            return;
X        else {
X            m = 0;
X            z = r;
X        }
X        }
X    } else if (((y) && (r % y != z % y)) || (r != z))
X        return;
X    }
X    if (m) {
X    z %= m;
X    if (z < 0)
X        z += m;
X    }
X
X    /*
X     * return information as required
X     */
X    if (list_args) {
X    return list(z,m);
X    } else {
X    return z;
X    }
X}
X
Xglobal lib_debug;
Xif (lib_debug >= 0) {
X    print "chrem(r1,m1 [,r2,m2 ...]) defined";
X    print "chrem(rlist [,mlist]) defined";
X}
SHAR_EOF
chmod 0644 calc2.9.0/lib/chrem.cal || echo "restore of calc2.9.0/lib/chrem.cal fails"
set `wc -c calc2.9.0/lib/chrem.cal`;Sum=$1
if test "$Sum" != "4346"
then echo original size 4346, current size $Sum;fi
echo "x - extracting calc2.9.0/lib/cryrand.cal (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/lib/cryrand.cal &&
X/*
X * cryrand - cryptographically strong pseudo-random number generator library
X */
X/*
X * Copyright (c) 1993 by Landon Curt Noll.  All Rights Reserved.
X *
X * Permission to use, copy, modify, and distribute this software and
X * its documentation for any purpose and without fee is hereby granted,
X * provided that the above copyright, this permission notice and text
X * this comment, and the disclaimer below appear in all of the following:
X *
X *    supporting documentation
X *    source copies
X *    source works derived from this source
X *    binaries derived from this source or from derived source
X *
X * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
X * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO
X * EVENT SHALL LANDON CURT NOLL BE LIABLE FOR ANY SPECIAL, INDIRECT OR
X * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
X * USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
X * OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
X * PERFORMANCE OF THIS SOFTWARE.
X *
X * chongo was here    /\../\
X */
X
X/*
X * These routines are written in the calc language.  At the time of this
X * notice, calc was at version 1.24.7 (We refer to calc, as in the C
X * program, not the Emacs subsystem).
X *
X * Calc is available by anonymous ftp from ftp.uu.net (137.39.1.9)
X * under the directory /pub/calc.
X *
X * If you can't get calc any other way, Email a request to my Email
X * address as shown below.
X *
X * Comments, suggestions, bug fixes and questions about these routines
X * are welcome.  Send Email to the address given below.
X *
X * Happy bit twiddling,
X *
X *            Landon Curt Noll
X *
X *            [e-mail]
X *            ...!{pyramid,sun,uunet}!hoptoad!chongo
X */
X
X/*
X * AN OVERVIEW OF THE FUNCTIONS:
X *
X * This calc library contains several pseudo-random number generators:
X *
X *   additive 55:
X *
X *    a55rand      - external interface to the additive 55 generator
X *    sa55rand  - seed the additive 55 generator
X *
X *    This is a generator based on the additive 55 generator as
X *    described in Knuth's "The Art of Computer Programming -
X *    Seminumerical Algorithms", vol 2, 2nd edition (1981),
X *    Section 3.2.2, page 27, Algorithm A.
X *
X *    The period and other properties of this generator make it very
X *    useful to 'seed' other generators.
X *
X *    This generator is used by other other generators to produce
X *    various internal values.  In fact, the lower 64 bits of seed
X *    given to other other generators is passed to sa55rand().
X *
X *   shuffle:
X *
X *    shufrand  - generate a pseudo-random value via a shuffle process
X *    sshufrand - seed the shufrand generator
X *    rand      - generate a pseudo-random value over a given range
X *    srand     - seed the random stream generator
X *
X *    This is a generator based on the shuffle generator as described in
X *    Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
X *    vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
X *
X *    The shuffle generator is fast and serves as a fairly good standard
X *    pseudo-random generator.
X *
X *    The shuffle generator is feed random values by the additive 55
X *    generator via a55rand().  Calling a55rand() or sa55rand() will
X *    affect the shuffle generator output.
X *
X *    The rand function is really another interface to the shuffle
X *    generator.  Unlike shufrand(), rand() can return a value of any
X *    given size.  The value returned by rand() may be confined to
X *    either a half open [0,a) (0 <= value < a) or closed interval
X *    [a,b] (a <= value <= b).  By default, a 64 bit value is returned.
X *
X *    Calling srand() simply calls sshufrand() with the same seed.
X *
X *   crypto:
X *
X *    cryrand   - produce a cryptographically strong pseudo-random number
X *    scryrand  - seed the crypto generator
X *    random    - produce a cryptographically strong pseudo-random number
X *            over a given range
X *     srandom   - seed random
X *
X *    This generator is described in the papers:
X *
X *        Blum, Blum, and Shub, "Comparison of Two Pseudorandom Number
X *        Generators", in Chaum, D. et. al., "Advances in Cryptology:
X *        Proceedings Crypto 82", pp. 61-79, Plenum Press, 1983.
X *
X *        "Probabilistic Encryption", Journal of Computer & System
X *        Sciences 28, pp. 270-299.
X *
X *    We also refer to this generator as the 'Blum' generator.
X *
X *    This generator is considered 'strong' in that it passes all
X *    polynomial-time statistical tests.  This means that there
X *    is no statistical test, which runs in polynomial time, that
X *    can distinguish the crypto generator from a truly random source.
X *
X *    The crypto generator is not as fast as most generators, though
X *    it is not painfully slow either.
X *
X *    One may fully seed this generator via scryrand().  Calling
X *    scryrand() with 1 or 3 arguments will result in the additive
X *    55 generator being seeded with the same seed.  Calling
X *     scryrand() with 4 arguments, where the first argument
X *    is >= 0 will also result in the additive 55 generator
X *    being seeded with the same seed.
X *
X *    The random() generator is really another interface to the
X *    crypto generator.  Unlike cryrand(), random() can return a
X *    value confined to either a half open (0 <= value < a) or closed
X *    interval (a <= value <= b).  By default, a 64 bit value is
X *    returned.
X *
X *    Calling srandom() simply calls scryrand(seed).  The additive
X *    55 generator will be seeded with the same seed.
X *
X * All generators come already seeded with precomputed initial constants.
X * Thus, it is not required to seed a generator before using it.
X *
X * Using a seed of '0' will reload generators with their initial states.
X * In the case of scryrand(), lengths of -1 must also be supplied.
X *
X *    sa55rand(0)
X *    sshufrand(0)
X *    srand(0)
X *    scryrand(0,-1,-1)
X *    srandom(0)
X *    scryrand(-1,ip,iq,ir)
X *
X * All of the above 'seed 0' calls are fairly fast.  In fact, passing
X * seeding with a non-zero seed, in the above cases, where seed is
X * not excessively large (many bits long), is also reasonably fast.
X * The 4 arg scryrand() call is fairly fast because no checking is
X * performed on the 'ip', 'iq', or 'ir' values.
X *
X * A call of scryrand(seed,len1,len2), with len1,len2 > 4, (3 args)
X * is a somewhat slow as the length args increase.  This type of
X * call selects 2 primes that are used internally by the crypto
X * generator.  A call of scryrand(seed,ip,iq,ir) where seed >= 0
X * is as slow as the 3 arg case.  See scryrand() for more information.
X *
X * A call of scryrand() (no args) may be used to quickly change the
X * internal state of the crypto and random generators.  Only one major
X * internal crypto generator value (a quadratic residue) is randomly
X * selected via the additive 55 generator.  The other 2 major internal
X * values (the 2 Blum primes) are preserved.  In this form, the additive
X * 55 is not seeded.
X *
X * Calling scryrand() with 1 or 3 args, or calling srandom() will leave
X * the additive 55 generator in a seeded state as if sa55rand(seed) has
X * been called.  Calling scryrand() with 4 args, where the first arg is
X * >0 will also leave the additive 55 generator in the same scryrand(seed)
X * state.  Calling scryrand() with no args will not seed the additive
X * 55 generator before or afterwards, however the additive 55 and shuffle
X * generators will be changed as a side effect of that call.
X *
X * The states of all generators (additive 55, shuffle and crypto) can be
X * saved and restored via the randstate() function.  Saving the state just
X * after * seeding a generator and restoring it later as a very fast way
X * to reseed a generator.
X *
X * As a bonus, the function 'nxtprime' is provided to produce a probable
X * prime number.
X *
X * TRUTH IN ADVERTISING:
X *
X * The word 'probable', in reference to the nxtprime() function, is used
X * because of an extremely extremely small chance that a composite (a
X * non-prime) may be returned.  In no cases will a prime be skipped.
X * For our purposes, this is sufficient as the chance of returning a
X * composite is much smaller than the chance that a hardware glitch
X * will cause nxtprime() to return a bogus result.
X *
X * Another "truth in advertising" issue is the use of the term
X * 'pseudo-random'.  All deterministic generators are pseudo random.
X * This is opposed to true random generators based on some special
X * physical device.
X *
X * The crypto generator is 'pseudo-random'.  There is no statistical
X * test, which runs in polynomial time, that can distinguish the crypto
X * generator from a truly random source.
X *
X * A final "truth in advertising" issue deals with how the magic numbers
X * found in this library were generated.  Detains can be found in the
X * various functions, while a overview can be found in the SOURCE FOR
X * MAGIC NUMBERS section below.
X *
X ****
X *
X * ON THE GENERATORS:
X *
X * The additive 55 generator has a good period, and is fast.  It is
X * reasonable as generators go, though there are better ones available.
X * We use it in seeding the crypto generator as its period and
X * other statistical properties are good enough for our purposes.
X *
X * This shuffle generator has a very good period, and is fast.  It is
X * fairly good as generators go, and is better than the additive 55
X * generator.  Casual direct use of the shuffle generator may be
X * acceptable.  Because of this, the interface to the shuffle generator,
X * but not the additive 55 generator, is advertised when this file is
X * loaded.
X *
X * The shuffle generator functions, shufrand() and rand() use the same
X * seed and tables.  The shuffle generator shuffles the values produced
X * by the additive 55 generator.  Calling or seeding the additive 55
X * generator will affect the output of the shuffle generator.
X *
X * The crypto generator is the best generator in this package.  It
X * produces a cryptographically strong pseudo-random bit sequence.
X * Internally, a fixed number of bits are generated after each
X * generator iteration.  Any unused bits are saved for the next call
X * to the generator.  The crypto generator is not too slow, though
X * seeding the generator from scratch is slow.  Shortcuts and
X * pre-computer seeds have been provided for this reason.  Use of
X * crypto should be more than acceptable for many applications.
X *
X * The crypto seed is in 4 parts, the additive 55 seed (lower 64
X * bits of seed), the shuffle seed (all but the lower 64 bits of
X * seed), and two lengths.  The two lengths specifies the minimum
X * bit size of two primes used internal to the crypto generator.
X * Not specifying the lengths, or using -1 will cause crypto to
X * use the default minimum lengths of 248 and 264 bits, respectively.
X *
X * The random() function just another interface to the crypto
X * generator.  Like rand(), random() provides an interval capability
X * (closed or open) as well as a 64 bit default return value.
X * The random() function as good as crypto, and produces numbers
X * that are equally cryptographically strong.  One may use the
X * seed functions srandom() or scryrand() for either the random()
X * or cryrand() functions.
X *
X * The seed for all of the generators may be of any size.  Only the
X * lower 64 bits of seed affect the additive 55 generator.  Bits
X * beyond the lower 64 bits affect the shuffle generators.  Excessively
X * large values of seed will result in increased memory usage as
X * well as a larger seed time for the shuffle and crypto generators.
X * See REGARDING SEEDS below, for more information.
X *
X * One may save and restore the state of all generators via the
X * randstate() function.
X *
X ****
X *
X * REGARDING SEEDS:
X *
X * Because the generators are interrelated, seeding one generator
X * will directly or indirectly affect the other generators.  Seeding
X * the shuffle generator seeds the additive 55 generator.  Seeding
X * the crypto generator seeds the shuffle generator.
X *
X * The seed of '0' implies that a generator should be seeded back
X * to its initial default state.
X *
X * For the moment, seeds < -1 are reserved for future use.  The
X * value of -1 is used as a special indicator to the fourth form
X * of scryrand(), and it not a real seed.
X *
X * A seed may be of any size.  The additive 55 generator uses only
X * the lower 64 bits, while the shuffle generator uses bytes beyond
X * the lower 64 bits.
X *
X * To help make the generator produced by seed S, significantly
X * different from S+1, seeds are scrambled prior to use.  The
X * function randreseed64() maps [0,2^64) into [0,2^64) in a 1-to-1
X * and onto fashion.
X *
X * The purpose of the randreseed64() is not to add security.  It
X * simply helps remove the human perception of the relationship
X * between the seed and the production of the generator.
X *
X * The randreseed64() process does not reduce the security of the
X * generators.  Every seed is converted into a different unique seed.
X * No seed is ignored or favored.
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, or 64 to
X * 780 bits long.
X *
X ****
X *
X * SOURCE OF MAGIC NUMBERS:
X *
X * Most of the magic constants used on this library ultimately are
X * based on the Rand book of random numbers.  The Rand book contains
X * 10^6 decimal digits, generated by a physical process.  This book,
X * produced by the Rand corporation in the 1950's is considered
X * a standard against which other generators may be measured.
X *
X * The Rand book of numbers was groups into groups of 20 digits.
X * The first 55 groups < 2^64 were used to initialize add55_init_tbl.
X * The size of 20 digits was used because 2^64 is 20 digits long.
X * The restriction of < 2^64 was used to prevent modulus biasing.
X * (see the note on  modulus biasing in rand()).
X *
X * The additive 55 generator during seeding is used 128 times to help
X * remove the initial seed state from the initial values produced.
X * The loop count of 128 was a power of 2 that permits each of the
X * 55 table entries to be processed at least twice.
X *
X * The function, randreseed64(), uses 4 primes to scramble 64 bits
X * into 64 bits.  These primes were also extracted from the Rand
X * book of numbers.  See sshufrand() for details.
X *
X * The default shuffle table size of 128 entries is the power of 2
X * that is longer than the 100 entries recommended by Knuth for
X * the shuffle algorithm (see the text cited in shufrand()).
X * We use a power of 2 shuffle table length so that the shuffle
X * process can select a table entry from a new additive 55 value
X * by extracting its top most bits.
X *
X * The quadratic residue search performed by cryres() starts at
X * a value that is in the interval [2^sqrpow,n-3), where 2^sqrpow
X * is the first power of 2 > n^(3/4) and n = p*q.  We look at 1009
X * random numbers in this interval for a quadratic residue.  If
X * none are found, sqrpow is decremented by 1 if sqrpow > 2.
X * In practice, decrementing sqrpow never happens.
X *
X * The use of n^(3/4) insures that the quadratic residue is
X * large, but not too large.  We want to avoid residues that are
X * near the square root of 'n', and are nearly 'n' itself (hence
X * the n-3 upper limit) as such residues are trivial or semi-trivial.
X *
X * The '1009' trials is simply an attempt to look for a while before
X * picking a new range.  The number '1009' comes from the first 4
X * digits of the rand book of numbers.
X *
X * Keeping sqrpow > 2 means that we do not look for quadratic
X * residues that are below 4.  This is in keeping with the
X * n-3 in the [2^sqrpow,n-3) interval.
X *
X * The final magic numbers: '248' and '264' are the exponents the
X * largest powers of 2 that are < the two default Blum primes 'p'
X * and 'q' used by the crypto generator.  The values of '248' and
X * '264' implies that the product n=p*q > 2^512.  Each iteration
X * of the crypto generator produces log2(log2(n=p*q)) random bits.
X * The crypto generator is the most efficient when n is slightly >
X * 2^(2^b).  The product n > 2^(2^9)) produces 9 bits as efficiently
X * as possible under the crypto generator process.
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 * The two len values differ slightly to avoid factorization attacks
X * that work on numbers that are a perfect square, or where the two
X * primes are nearly the same.  I elected to have the sizes differ
X * by 3% of the product size.  The difference between '248' and
X * '264', is '16', which is ~3.125% of '512'.  Now 3% of '512' is
X * '~15.36', and the next largest whole number is '16'.
X *
X * The product n=p*q > 2^512 implies a product if at least 155 digits.
X * A product of two primes that is at least 155 digits is slightly
X * beyond Number Theory and computer power of Nov 1992, though this
X * will likely change in the future.
X *
X * Again, the ability (or lack thereof) to factor 'n=p*q' does not
X * directly relate to the strength of the crypto generator.  We
X * selected n=p*q > 2^512 mainly because '512 was a power of 2 and
X * only slightly because it is up in the range where it is difficult
X * to factor.
X *
X ****
X *
X * FOR THE PARANOID:
X *
X * The truly paranoid might suggest that my claims in the MAGIC NUMBERS
X * section are a lie intended to entrap people.  Well they are not, but
X * you need not take my word for it.
X *
X * The random numbers from the Rand book of random numbers can be
X * verified by anyone who obtains the book.  As these numbers were
X * created before I (Landon Curt Noll) was born (you can look up
X * my birth record if you want), I claim to have no possible influence
X * on their generation.
X *
X * There is a very slight chance that the electronic copy of the
X * Rand book that I was given access to differs from the printed text.
X * I am willing to provide access to this electronic copy should
X * anyone wants to compare it to the printed text.
X *
X * One might object to the complexity of the seed scramble/mapping
X * via the randreseed64() function.  The randreseed64() function maps:
X *
X *    1 ==> 4967126403401436567
X *
X * calling srand(1) with the randreseed64() process would be the same
X * as calling srand(4967126403401436567) without it.  No extra security
X * is gained or reduced by using the randreseed64() process.  The meaning
X * of seeds are exchanged, but not lost or favored (used by more than
X * one input seed).
X *
X * One could take issue with my selection of the default sizes '248'
X * and '264'.   As far as I know, 155 digits (512 bits) is beyond the
X * state of the art of Number Theory and Computation as of 01 Sep 92.
X * It will likely be true that 155 digit products of two primes could
X * come within reach in the next few years, but so what?  If you are
X * truly paranoid, why would you want to use the default seed, which
X * is well known?
X *
X * The paranoid today might consider using the lengths of at least '345'
X * and '325' will produce a product of two primes that is 202 digits.
X * (the 2nd and 3rd args of scryrand > 345 & >325 respectively)  Factoring
X * 200+ digit product of two primes is well beyond the current hopes of
X * Number Theory and Computer power, though even this limit may be passed
X * someday.
X *
X * If all the above fails to pacify the truly paranoid, then one may
X * select by some independent means, 2 Blum primes (primes mod 4 == 3,
X * p < q), and a quadratic residue if p*q.  Then by calling:
X *
X *    scryrand(-1, p, q, r)
X *
X * and then calling cryrand() or random(), one may bypass all magic
X * numbers and use the pure generator.
X *
X * Note that randstate() may also be used by the truly paranoid.
X * Even though it holds state for the other generators, their states
X * are independent.
X *
X ****
X *
X * GOALS:
X *
X * The goals of this package are:
X *
X *    all magic numbers are explained
X *
X *        I distrust systems with constants (magic numbers) and tables
X *        that have no justification (e.g., DES).  I believe that I have
X *        done my best to justify all of the magic numbers used.
X *
X *     full documentation
X *
X *        You have this source file, plus background publications,
X *        what more could you ask?
X *
X *    large selection of seeds
X *
X *        Seeds are not limited to a small number of bits.  A seed
X *        may be of any size.
X *
X *    the strength of the generators may be tuned to meet the application need
X *
X *        By using the appropriate seed arguments, one may increase
X *        the strength of the generator to suit the need of the
X *        application.  One does not have just a few levels.
X *
X * This calc lib file is intended for demonstration purposes.  Writing
X * a C program (with possible assembly or libmp assist) would produce
X * a faster generator.
X *
X * Even though I have done my best to implement a good system, you still
X * must use these routines your own risk.
X *
X * Share and enjoy!  :-)
X */
X
X
X/*
X * These constants are used by all of the generators in various direct and
X * indirect forms.
X */
Xstatic cry_seed = 0;            /* master seed */
Xstatic cry_mask = 0xffffffffffffffff;    /* 64 bit work mask */
X
X
X/*
X * Initial magic values for the additive 55 generator - used by sa55rand()
X *
X * These values were generated from the Rand book of random numbers.
X * In groups of 20 digits, we took the first 55 groups < 2^64.
X */
Xstatic mat add55_init_tbl[55] = {
X    10097325337652013586,    8422689531964509303,
X    9376707153831131165,    12807999708015736147,
X    12171768336606574717,    2051656926866574818,
X    3529647783580834282,    13746700781847540610,
X    17468509505804776974,    14385537637435099817,
X    14225685144642756788,    11100020401286074697,
X    7207317906119690446,    15474452669527079953,
X    16868487670307112059,    4493524947524633824,
X    13021248927856520106,    15956600001874392423,
X    1758753794041921585,    1540354560501451176,
X    5335129695612719255,    9973334408846123356,
X    2295368703230757546,    15020099946907494138,
X    10518216150184876938,    9188200973282539527,
X    4220863048338987374,    682273982071453295,
X    7706178130835869910,    4618975533122308420,
X    397583911260717646,    5686731560708285046,
X    10123916228549657560,    1304775865627110086,
X    15501295782182641134,    3061180729620744156,
X    6958929830512809719,    10850627469959910507,
X    13499063195307571839,    6410193623982098952,
X    4111084083850807341,    17719042079595449953,
X    5462692006544395659,    18288274374963224041,
X    8337656769629990836,    7477446061798548911,
X    9815931464890815877,    6913451974267278601,
X    11883095286301198901,    14974403441045516019,
X    14210337129134237821,    12883973436502761184,
X    4285013921797415077,    16435915296724552670,
X    3742838738308012451
X};
X
X
X/*
X * additive 55 tables - used by a55rand() and sa55rand()
X */
Xstatic add55_j = 23;        /* the first walking table pointer */
Xstatic add55_k = 54;        /* the second walking table pointer */
Xstatic add55_seed64 = 0;    /* lower 64 bits of the reseeded seed */
Xstatic mat add55_tbl[55];    /* additive 55 state table */
X
X
X/*
X * cryobj - cryptographic pseudo-random state object
X */
Xobj cryobj {                            \
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
X
X/*
X * initial cryptographic pseudo-random values - used by scryrand()
X *
X * These values are what the crypto generator is initialized with
X * with this library first read.  These values may be reproduced the
X * hard way by calling scryrand(0,248,264) or scryrand(0,-1,-1).
X *
X * We will build up these values a piece at a time to avoid long lines
X * that are difficult to send via Email.
X */
X/* p, first Blum prime */
Xstatic cryrand_init_p = 0x1aa9e726a735044;
Xcryrand_init_p <<= 200;
Xcryrand_init_p |= 0x73b7457c5297122310880fcbfa8d4e38380a00396d533300bb;
X/* q, second Blum prime */
Xstatic cryrand_init_q = 0xa62ee0481aa12059c3;
Xcryrand_init_q <<= 200;
Xcryrand_init_q |= 0x79ef44e72ff58ea70cacabbe9d264a0b16db72117d96f77e17;
X/* quadratic residue of n=p*q */
Xstatic cryrand_init_r = 0x3d214853f9a5bb4b12f467c9052129a9;
Xcryrand_init_r <<= 200;
Xcryrand_init_r |= 0xd215cc6b3c2eae8c7090591b16d8044a886b3c6a58759b1a07;
Xcryrand_init_r <<= 200;
Xcryrand_init_r |= 0x2b50a914b42e1b6f9703be86742837c4bc637804c2dc668c5b;
X
X/*
X * cryptographic pseudo-random values - used by cryrand() and scryrand()
X */
X/* p, first Blum prime */
Xstatic cryrand_p = cryrand_init_p;
X/* q, second Blum prime */
Xstatic cryrand_q = cryrand_init_q;
X/* n = p*q */
Xstatic cryrand_n = cryrand_p*cryrand_q;
X/* quad residue of n */
Xstatic cryrand_r = cryrand_init_r;
X/* this cryrand() running exp used in computing crypto good bits */
Xstatic cryrand_exp = cryrand_r;
X/* good crypto bits unused from the last cryrand() call */
Xstatic cryrand_left = 0;
X/* the value cryrand_left contains cryrand_bitcnt crypto good bits */
Xstatic cryrand_bitcnt = 0;
X
X
X/*
X * shufrand - shuffle generator constants
X */
Xstatic shuf_size = 128;            /* entries in shuffle table */
Xstatic shuf_shift = (64-highbit(shuf_size));    /* shift a55 value to get tbl indx */
Xstatic shuf_y;                /* Y value (previous output) */
Xstatic mat shuf_tbl[shuf_size];        /* shuffle state table */
X
X
X/*
X * products of consecutive primes - used by nxtprime()
X *
X * We compute these products now to avoid recalculating them on each call.
X * These values help weed out numbers that have a prime factor < 1000.
X */
Xstatic nxtprime_pfact10 = pfact(10);
Xstatic nxtprime_pfact100 = pfact(100)/nxtprime_pfact10;
Xstatic nxtprime_pfact1000 = pfact(1000)/nxtprime_pfact100;
X
X
X/*
X * a55rand - additive 55 pseudo-random number generator
X *
X * returns:
X *    A number in the half open interval [0,2^64)
X *
X * This function implements the additive 55 pseudo-random number generator.
X *
X * This is a generator based on the additive 55 generator as described in
X * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
X * vol 2, 2nd edition (1981), Section 3.2.2, page 27, Algorithm A.
X *
X * This function is called from the shuffle generator function shufrand().
X *
X * NOTE: This is NOT Blum's method.  This is used by Blum's method to
X *       generate some internal values.
X *
X * NOTE: If you are looking for a fast generator and do not need a crypto
X *      strong generator, you should not use this generator.  Use of
X *     the shuffle generator functions shufrand() and rand() should
X *     be considered instead.
X */
Xdefine
Xa55rand()
X{
X    local ret;            /* the pseudo-random number to return */
X
X    /* generate the next value using the additive 55 generator method */
X    ret = add55_tbl[add55_k] + add55_tbl[add55_j];
X    ret &= cry_mask;
X    add55_tbl[add55_k] = ret;
X
X    /* post-process out data with the seed */
X    ret = xor(ret, add55_seed64);
X
X    /* step the pointers */
X    --add55_j;
X    if (add55_j < 0) {
X    add55_j = 54;
X    }
X    --add55_k;
X    if (add55_k < 0) {
X    add55_k = 54;
X    }
X
X    /* return the new pseudo-random number */
X    return(ret);
X}
X
X
X/*
X * sa55rand - initialize the additive 55 pseudo-random generator
X *
X * given:
X *    seed
X *
X * returns:
X *    old_seed
X *
X * This function seeds the additive 55 pseudo-random generator.
X *
X * This is a generator based on the additive 55 generator as described in
X * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
X * vol 2, 2nd edition (1981), Section 3.2.2, page 27, Algorithm A.
X *
X * Unlike Knuth's description, we will let a seed post process our data.
X *
X * We do not apply the seed processing to the additive 55 table
X * data as this would disturb its pseudo-random generator properties.
X * Instead, we xor the output with the low 64 bits of seed.
X *
X * If seed == 0:
X *
X *    This function produces values in exactly the same way
X *    described by Knuth.
X *
X * else:
X *
X *    Each value produced is xor-ed by a 64 bit value.  This value
X *    is the result of randreseed64(rand), which produces a 64
X *    bit value.
X *
X * Anyone comfortable with seed == 0 should also be comfortable with
X * non-zero seeds.  A non-zero seeded sequence will produce a values
X * that have the exact same pseudo-random properties as the algorithm
X * described by Knuth.  I.e., the sequence, while different, is as good
X * (or bad) as the sequence produced by a seed of 0.
X *
X * This function updates both the cry_seed and a55_seed64 global values.
X *
X * This function is called from the shuffle generator seed function sshufrand().
X *
X * NOTE: This is NOT Blum's method.  This is used by Blum's method to
X *       generate some internal values.
X *
X * NOTE: If you are looking for a fast generator and do not need a crypto
X *      strong generator, you should not use this generator.  Use of the
X *     shuffle generator seed functions sshufrand() and srand() should
X *     be considered instead.
X */
Xdefine
Xsa55rand(seed)
X{
X    local oldseed;        /* previous seed */
X    local junk;            /* discards the first few numbers */
X    local j;
X
X    /* firewall */
X    if (!isint(seed)) {
X    quit "bad arg: arg must be an integer";
X    }
X    if (seed < 0) {
X    quit "bad arg: seed < 0 is reserved for future use";
X    }
X
X    /* misc table setup */
X    oldseed = cry_seed;                /* remember the previous seed */
X    cry_seed = seed;                /* save the new seed */
X    add55_tbl = add55_init_tbl;    /* reload the table */
X    add55_j = 23;               /* reset first walking table pointer */
X    add55_k = 54;               /* reset second walking table pointer */
X
X    /* obtain our 64 bit xor seed */
X    add55_seed64 = randreseed64(seed);
X
X    /* spin the pseudo-random number generator a while */
X    for (j=0; j < 128; ++j) {
X    junk = a55rand();
X    }
X
X    /* return the old seed */
X    return(oldseed);
X}
X
X
X/*
X * shufrand - implement the shuffle pseudo-random number generator
X *
X * returns:
X *    A number in the half open interval [0,2^64)
X *
X * This function implements the shuffle number generator.
X *
X * This is a generator based on the shuffle generator as described in
X * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
X * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
X *
X * The function rand() calls this function.
X *
X * NOTE: This is NOT Blum's method.  This is used by Blum's method to
X *       generate some internal values.
X *
X * NOTE: If you do not need a crypto strong pseudo-random generator
X *     this function may very well serve your needs.
X */
Xdefine
Xshufrand()
X{
X    local j;        /* table index to replace */
X
X    /*
X     * obtain a new random value
X     * determine the table entry to shuffle
X     * shuffle out the value we will return
X     */
X    shuf_y = shuf_tbl[j = (shuf_y >> shuf_shift)];
X
X    /* shuffle in the new random value */
X    shuf_tbl[j] = a55rand();
X
X    /* return the shuffled out value */
X    return (shuf_y);
X}
X
X
X/*
X * sshufrand - seed the shuffle pseudo-random generator
X *
X * given:
X *    a seed
X *
X * returns:
X *    the previous seed
X *
X * This function implements the shuffle number generator.
X *
X * This is a generator based on the shuffle generator as described in
X * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
X * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
X *
X * The low 64 bits of seed are used by the additive 55 generator.
X * See the sa55rand() function for details.  The remaining bits of seed
X * are used to perform an initial shuffle on the shuffle state table.
X * The size of the seed also determines the size of the shuffle table.
X *
X * The shuffle table size is always a power of 2, and is at least 128
X * entries long.  Let the table size be:
X *
X *    shuf_size = 2^shuf_pow
X *
X * The number of ways one could shuffle that table is:
X *
X *    (2^shuf_pow)!
X *
X * We select the smallest 'shuf_pow' (and thus the size of the shuffle table)
X * such that the following are true:
X *
X *    (2^shuf_pow)! >= (seed / 2^64)    and    2^shuf_pow >= 128
X *
X * Given that we now have the table size of 'shuf_size', we must go about
X * loading the table and shuffling it.
X *
X * Loading is easy, we will generate random values via the additive 55
X * generator (a55rand()) and load them into successive entries.
X *
X * We enumerate the (2^shuf_pow)! values via:
X *
X *    shuf_seed = 2*(U[2] + 3*(U[3] + 4*(U[4] + ...
X *              + (U[shuf_pow-1]*(shuf_pow-1)) ... )))
X *    0 <= U[j] < j
X *
X * We swap the swap table entries shuf_tbl[U[j]] & shuf_tbl[j-1] for all
X * 1 < j < shuf_pow.
X *
X * We will convert 'seed / 2^64' into shuf_seed, by applying the 64 bit
X * scramble function on 64 bit chunks of 'seed / 2^64'.
X *
X * The function srand() calls this function.
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 is NOT Blum's method.  This is used by Blum's method to
X *       generate some internal values.
X *
X * NOTE: If you do not need a crypto strong pseudo-random generator
X *     this function may very well serve your needs.
X */
Xdefine
Xsshufrand(seed)
X{
X    local shuf_pow;        /* power of two - log2(shuf_size) */
X    local shuf_seed;        /* seed bits beyond low 64 bits */
X    local oldseed;        /* previous seed */
X    local shift;        /* shift factor to access 64 bit chunks */
X    local swap_indx;        /* what to swap shuf_tbl[0] with */
X    local rval;            /* random value form additive 55 generator */
X    local j;
X
X    /* firewall */
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
X    /*
X     * seed the additive 55 generator
X     */
X    oldseed = sa55rand(seed);
X
X    /*
X     * form the shuffle table size and constants
X     */
X    shuf_seed = seed >> 64;
X    for (shuf_pow = 7; shuf_seed > (j=fact(1<<(shuf_pow))) && shuf_pow < 64; \
X     ++shuf_pow) {
X    }
X    shuf_size = (1 << shuf_pow);
X    shuf_shift = 64 - shuf_pow;
X    /* reallocate the shuffle table */
X    mat shuf_tbl[shuf_size];
X
X    /*
X     * scramble the seed above the low 64 bits
X     */
X    if (shuf_seed > 0) {
X    j = 0;
X    for (shift=0; shift < highbit(shuf_seed)+1; shift += 64) {
X        j |= (randreseed64(shuf_seed >> shift) << shift);
X    }
X    shuf_seed = j;
X    }
X
X    /*
X     * load the shuffle table
X     */
X    for (j=0; j < shuf_size; ++j) {
X    shuf_tbl[j] = a55rand();
X    }
X    shuf_y = a55rand();        /* get the next Y value */
X
X    /*
X     * We will shuffle based the process outlined in:
X     *
X     * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
X     * vol 2, 2nd edition (1981), Section 3.4.2, page 139, Algorithm P.
X     *
X     * Here, we will let j run over the range [0,shuf_size) instead of
X     * [shuf_size,0) as outlined in algorithm P.  We will also generate
X     * U values from shuf_seed.
X     */
X    j = 0;
X    while (shuf_seed > 0 && ++j < shuf_size) {
X
X    /* determine what index we will swap with the '0' index */
X    quomod(shuf_seed, (j+1), shuf_seed, swap_indx);
X
X    /* swap table entries, if needed */
X    if (swap_indx != j) {
X        swap(shuf_tbl[j], shuf_tbl[swap_indx]);
X    }
X    }
X
X    /*
X     * run the generator for twice the table size
X     */
X    for (j=0; j < shuf_size*2; ++j) {
X    rval = shufrand();
X    }
X
X    /* return the old seed */
X    return (oldseed);
X}
X
X
X/*
X * rand - generate a pseudo-random value over a given range via additive 55
X *
X * usage:
X *    rand()         - generate a pseudo-random integer >=0 and < 2^64
X *    rand(a)     - generate a pseudo-random integer >=0 and < a
X *    rand(a,b)    - generate a pseudo-random integer >=a and <= b
X *
X * returns:
X *    a large pseudo-random integer over a give range (see usage)
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 * shufrand().
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 * The input values a and b, if given, must be integers.
X *
X * This generator is simply a different interface to the shuffle generator.
X * calling shufrand(), or seeding via sshufrand() will affect the output
X * of this function.
X *
X * NOTE: Unlike cryrand(), this function does not preserve unused bits for
X *     use by the next function call.
X *
X * NOTE: The Un*x rand() function returns only 16 bit or 31 bits, while we
X *     return a number of any given size (default is 64 bits).
X *
X * NOTE: This is NOT Blum's method.  This is used by Blum's method to
X *       generate some internal values.
X *
X * NOTE: If you do not need a crypto strong pseudo-random generator
X *     this function may very well serve your needs.
X */
Xdefine
Xrand(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 ret;            /* pseudo-random value */
X    local fullwords;        /* number of full 64 bit chunks in ret */
X    local finalmask;        /* mask of bits in final chunk of range */
X    local j;
X
X    /*
X     * setup and special cases
X     */
X    /* deal with the rand() case */
X    if (isnull(a) && isnull(b)) {
X    /* no args means same range as additive 55 generator */
X    return(a55rand());
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 rand(a,b), where a < b into rand(b,a) */
X    range = b-a+1;
X    offset = a;
X    }
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    fullwords = highbit(range-1)//64;
X    finalmask = (1 << (1+(highbit(range-1)%64)))-1;
X
X    /*
X     * loop until we get a value that is in range
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 frequently
X     * because the numbers [81,100) mod 80 wrap back into [0,20).
X     */
X    do {
X    /* load up all lower full 64 bit chunks with pseudo-random bits */
X    ret = 0;
X    for (j=0; j < fullwords; ++j) {
X        ret = (ret << 64) + shufrand();
X    }
X
X    /* load the highest chunk */
X    ret <<= (highbit(finalmask)+1);
X    ret |= (shufrand() >> (64-highbit(finalmask)-1));
X    } while (ret >= range);
X
X    /*
X     * return the adjusted range value
X     */
X    return(ret+offset);
X}
X
X
X/*
X * srand - seed the pseudo-random additive 55 generator
X *
X * input:
X *    seed
X *
X * returns:
X *    old_seed
X *
X * This function actually seeds the shuffle generator (and indirectly
X * the additive 55 generator used by rand() and a55rand().
X *
X * See sshufrand() and sa55rand() for information about a seed.
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 is NOT Blum's method.  This is used by Blum's method to
X *       generate some internal values.
X *
X * NOTE: If you do not need a crypto strong pseudo-random generator
X *     this function may very well serve your needs.
X */
Xdefine
Xsrand(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(sshufrand(seed));
X}
X
X
X/*
X * cryrand - cryptographically strong pseudo-random number generator
X *
X * usage:
X *    cryrand(len)
X *
X * given:
X *    len        number of pseudo-random bits to generate
X *
X * returns:
X *    a cryptographically strong pseudo-random number of len bits
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() or srandom() will clear out all
X * unused bits.  Thus:
X *
X *    scryrand(0);            <-- restore generator to initial state
X *    cryrand(16);            <-- 16 bits
X *
X * will produce the same value as:
X *
X *    scryrand(0);            <-- restore generator to initial state
X *    cryrand(4)<<12 | cryrand(12);    <-- 4+12 = 16 bits
X *
X * and will produce the same value as:
X *
X *    scryrand(0);            <-- restore generator to initial state
X *    cryrand(3)<<13 | cryrand(7)<<6 | cryrand(6);    <-- 3+7+6 = 16 bits
X *
X * The crypto generator is not as fast as most generators, though it is not
X * painfully slow either.
X *
X * NOTE: This function is the Blum cryptographically strong
X *     pseudo-random number generator!
X */
Xdefine
Xcryrand(len)
X{
X    local goodbits;    /* the number of good bits generated each pass */
X    local goodmask;    /* mask for the low order good bits */
X    local randval;    /* pseudo-random value being generated */
X
X    /*
X     * firewall
X     */
X    if (!isint(len) || len < 1) {
X    quit "bad arg: len must be an integer > 0";
X    }
X
X    /*
X     * Determine how many bits may be generated each pass.
X     *
X     * The result by Alexi et. al., says that the log2(log2(n=p*q))
X     * least significant bits are secure, where log2(x) is log base 2.
X     */
X    goodbits = highbit(highbit(cryrand_n));
X    goodmask = (1 << goodbits)-1;
X
X    /*
X     * If we have bits left over from the previous call, collect
X     * them now.
X     */
X    if (cryrand_bitcnt > 0) {
X
X    /* case where the left over bits are enough for this call */
X    if (len <= cryrand_bitcnt) {
X
X        /* we need only len bits */
X        randval = (cryrand_left >> (cryrand_bitcnt-len));
X
X        /* save the unused bits for later use */
X        cryrand_left &= ((1 << (cryrand_bitcnt-len))-1);
X
X        /* save away the number of bits that we will not use */
X        cryrand_bitcnt -= len;
X
X        /* return our complete result */
X        return(randval);
X
X    /* case where we need more than just the left over bits */
X    } else {
X
X        /* clear out the number of left over bits */
X        len -= cryrand_bitcnt;
X        cryrand_bitcnt = 0;
X
X        /* collect all of the left over bits for now */
X        randval = cryrand_left;
X    }
X
X    /* case where we have no previously left over bits */
X    } else {
X    randval = 0;
X    }
X
X    /*
X     * Pump out len cryptographically strong pseudo-random bits,
X     * 'goodbits' at a time using Blum's process.
X     */
X    while (len >= goodbits) {
X
X    /* generate the bits */
X    cryrand_exp = (cryrand_exp^2) % cryrand_n;
X    randval <<= goodbits;
X    randval |= (cryrand_exp & goodmask);
X
X    /* reduce the need count */
X    len -= goodbits;
X    }
X
X    /* if needed, save the unused bits for later use */
X    if (len > 0) {
X
X    /* generate the bits */
X    cryrand_exp = (cryrand_exp^2) % cryrand_n;
X    randval <<= len;
X    randval |= ((cryrand_exp&goodmask) >> (goodbits-len));
X
X    /* save away the number of bits that we will not use */
X    cryrand_left = cryrand_exp & ((1 << (goodbits-len))-1);
X    cryrand_bitcnt = goodbits-len;
X    }
X
X    /*
X     * return our pseudo-random bits
X     */
X     return(randval);
X}
X
X
X/*
X * scryrand - seed the cryptographically strong pseudo-random number generator
X *
X * usage:
X *    scryrand(seed)
X *    scryrand()
X *    scryrand(seed,len1,len2)
X *    scryrand(seed,ip,iq,ir)
X *
X * input:
X *    [seed        pseudo-random seed
X *    [len1 len2]    minimum bit length of the Blum primes 'p' and 'q'
X *            -1 => default lengths
X *    [ip iq ir]    Initial search values for Blum primes 'p', 'q' and
X *            a quadratic residue 'r'
X *
X * returns:
X *    the previous seed
X *
X *
X * This function will seed and setup the generator needed to produce
X * cryptographically strong pseudo-random numbers.  See the function
X * a55rand() and sshufrand() for information about how 'seed' works.
X *
X * The first form of this function are fairly fast if the seed is not
X * excessively large.  The second form is also fairly fast if the internal
X * primes are not too large.  The third form, can take a long time to call.
X * (see below)   The fourth form, if the 'seed' arg is not -1, can take
X * as long as the third form to call.  If the fourth form is called with
X * a 'seed' arg of -1, then it is fairly fast.
X *
X * Calling scryrand() with 1 or 3 args (first and third forms), or
X * calling srandom(), or calling scryrand() with 4 args with the first
X * arg >0, will leave the shuffle generator in a seeded state as if
X * sshufrand(seed) has been called.
X *
X * Calling scryrand() with no args will not seed the shuffle generator,
X * before or afterwards, however the shuffle generator will have been
X * changed as a side effect of that call.
X *
X * Calling scryrand() with 4 args where the first arg is 0 or '-1'
X * will not change the other generators.
X *
X *
X * First form of call:  scryrand(seed)
X *
X * The first form of this function will seed the shuffle generator
X * (via srand).  The default precomputed constants will be used.
X *
X *
X * Second form of call:  scryrand()
X *
X * Only a new quadratic residue of n=p*q is recomputed.  The previous prime
SHAR_EOF
echo "End of part 16"
echo "File calc2.9.0/lib/cryrand.cal is continued in part 17"
echo "17" > s2_seq_.tmp
exit 0