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

Text File | 1992-05-09 | 34.6 KB | 1,386 lines

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

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

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 16 (of 21)."
# Contents:  zmod.c
# Wrapped by dbell@elm on Tue Feb 25 15:21:14 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'zmod.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'zmod.c'\"
else
echo shar: Extracting \"'zmod.c'\" \(32468 characters\)
sed "s/^X//" >'zmod.c' <<'END_OF_FILE'
X/*
X * Copyright (c) 1992 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Routines to do modulo arithmetic both normally and also using the REDC
X * algorithm given by Peter L. Montgomery in Mathematics of Computation,
X * volume 44, number 170 (April, 1985).  For multiple multiplies using
X * the same large modulus, the REDC algorithm avoids the usual division
X * by the modulus, instead replacing it with two multiplies or else a
X * special algorithm.  When these two multiplies or the special algorithm
X * are faster then the division, then the REDC algorithm is better.
X */
X
X#include "math.h"
X
X
X#define    POWBITS    4        /* bits for power chunks (must divide BASEB) */
X#define    POWNUMS    (1<<POWBITS)    /* number of powers needed in table */
X
X
XLEN _pow2_ = POW_ALG2;        /* modulo size to use REDC for powers */
XLEN _redc2_ = REDC_ALG2;    /* modulo size to use second REDC algorithm */
X
Xstatic REDC *powermodredc = NULL;    /* REDC info for raising to power */
X
X#if 0
Xextern void zaddmod proto((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
Xextern void znegmod proto((ZVALUE z1, ZVALUE z2, ZVALUE *res));
X
X/*
X * Multiply two numbers together and then mod the result with a third number.
X * The two numbers to be multiplied can be negative or out of modulo range.
X * The result will be in the range 0 to the modulus - 1.
X */
Xvoid
Xzmulmod(z1, z2, z3, res)
X    ZVALUE z1;        /* first number to be multiplied */
X    ZVALUE z2;        /* second number to be multiplied */
X    ZVALUE z3;        /* number to take mod with */
X    ZVALUE *res;        /* result */
X{
X    ZVALUE tmp;
X    FULL prod;
X    FULL digit;
X    BOOL neg;
X
X    if (iszero(z3) || isneg(z3))
X        error("Mod of non-positive integer");
X    if (iszero(z1) || iszero(z2) || isunit(z3)) {
X        *res = _zero_;
X        return;
X    }
X
X    /*
X     * If the modulus is a single digit number, then do the result
X     * cheaply.  Check especially for a small power of two.
X     */
X    if (istiny(z3)) {
X        neg = (z1.sign != z2.sign);
X        digit = z3.v[0];
X        if ((digit & -digit) == digit) {    /* NEEDS 2'S COMP */
X            prod = ((FULL) z1.v[0]) * ((FULL) z2.v[0]);
X            prod &= (digit - 1);
X        } else {
X            z1.sign = 0;
X            z2.sign = 0;
X            prod = (FULL) zmodi(z1, (long) digit);
X            prod *= (FULL) zmodi(z2, (long) digit);
X            prod %= digit;
X        }
X        if (neg && prod)
X            prod = digit - prod;
X        itoz((long) prod, res);
X        return;
X    }
X
X    /*
X     * The modulus is more than one digit.
X     * Actually do the multiply and divide if necessary.
X     */
X    zmul(z1, z2, &tmp);
X    if (ispos(tmp) && ((tmp.len < z3.len) || ((tmp.len == z3.len) &&
X        (tmp.v[tmp.len-1] < z2.v[z3.len-1]))))
X    {
X        *res = tmp;
X        return;
X    }
X    zmod(tmp, z3, res);
X    freeh(tmp.v);
X}
X
X
X/*
X * Square a number and then mod the result with a second number.
X * The number to be squared can be negative or out of modulo range.
X * The result will be in the range 0 to the modulus - 1.
X */
Xvoid
Xzsquaremod(z1, z2, res)
X    ZVALUE z1;        /* number to be squared */
X    ZVALUE z2;        /* number to take mod with */
X    ZVALUE *res;        /* result */
X{
X    ZVALUE tmp;
X    FULL prod;
X    FULL digit;
X
X    if (iszero(z2) || isneg(z2))
X        error("Mod of non-positive integer");
X    if (iszero(z1) || isunit(z2)) {
X        *res = _zero_;
X        return;
X    }
X
X    /*
X     * If the modulus is a single digit number, then do the result
X     * cheaply.  Check especially for a small power of two.
X     */
X    if (istiny(z2)) {
X        digit = z2.v[0];
X        if ((digit & -digit) == digit) {    /* NEEDS 2'S COMP */
X            prod = (FULL) z1.v[0];
X            prod = (prod * prod) & (digit - 1);
X        } else {
X            z1.sign = 0;
X            prod = (FULL) zmodi(z1, (long) digit);
X            prod = (prod * prod) % digit;
X        }
X        itoz((long) prod, res);
X        return;
X    }
X
X    /*
X     * The modulus is more than one digit.
X     * Actually do the square and divide if necessary.
X     */
X    zsquare(z1, &tmp);
X    if ((tmp.len < z2.len) ||
X        ((tmp.len == z2.len) && (tmp.v[tmp.len-1] < z2.v[z2.len-1]))) {
X            *res = tmp;
X            return;
X    }
X    zmod(tmp, z2, res);
X    freeh(tmp.v);
X}
X
X
X/*
X * Add two numbers together and then mod the result with a third number.
X * The two numbers to be added can be negative or out of modulo range.
X * The result will be in the range 0 to the modulus - 1.
X */
Xstatic void
Xzaddmod(z1, z2, z3, res)
X    ZVALUE z1;        /* first number to be added */
X    ZVALUE z2;        /* second number to be added */
X    ZVALUE z3;        /* number to take mod with */
X    ZVALUE *res;        /* result */
X{
X    ZVALUE tmp;
X    FULL sumdigit;
X    FULL moddigit;
X
X    if (iszero(z3) || isneg(z3))
X        error("Mod of non-positive integer");
X    if ((iszero(z1) && iszero(z2)) || isunit(z3)) {
X        *res = _zero_;
X        return;
X    }
X    if (istwo(z2)) {
X        if ((z1.v[0] + z2.v[0]) & 0x1)
X            *res = _one_;
X        else
X            *res = _zero_;
X        return;
X    }
X    zadd(z1, z2, &tmp);
X    if (isneg(tmp) || (tmp.len > z3.len)) {
X        zmod(tmp, z3, res);
X        freeh(tmp.v);
X        return;
X    }
X    sumdigit = tmp.v[tmp.len - 1];
X    moddigit = z3.v[z3.len - 1];
X    if ((tmp.len < z3.len) || (sumdigit < moddigit)) {
X        *res = tmp;
X        return;
X    }
X    if (sumdigit < 2 * moddigit) {
X        zsub(tmp, z3, res);
X        freeh(tmp.v);
X        return;
X    }
X    zmod(tmp, z2, res);
X    freeh(tmp.v);
X}
X
X
X/*
X * Subtract two numbers together and then mod the result with a third number.
X * The two numbers to be subtract can be negative or out of modulo range.
X * The result will be in the range 0 to the modulus - 1.
X */
Xvoid
Xzsubmod(z1, z2, z3, res)
X    ZVALUE z1;        /* number to be subtracted from */
X    ZVALUE z2;        /* number to be subtracted */
X    ZVALUE z3;        /* number to take mod with */
X    ZVALUE *res;        /* result */
X{
X    if (iszero(z3) || isneg(z3))
X        error("Mod of non-positive integer");
X    if (iszero(z2)) {
X        zmod(z1, z3, res);
X        return;
X    }
X    if (iszero(z1)) {
X        znegmod(z2, z3, res);
X        return;
X    }
X    if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
X        (z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0)) {
X            *res = _zero_;
X            return;
X    }
X    z2.sign = !z2.sign;
X    zaddmod(z1, z2, z3, res);
X}
X
X
X/*
X * Calculate the negative of a number modulo another number.
X * The number to be negated can be negative or out of modulo range.
X * The result will be in the range 0 to the modulus - 1.
X */
Xstatic void
Xznegmod(z1, z2, res)
X    ZVALUE z1;        /* number to take negative of */
X    ZVALUE z2;        /* number to take mod with */
X    ZVALUE *res;        /* result */
X{
X    int sign;
X    int cv;
X
X    if (iszero(z2) || isneg(z2))
X        error("Mod of non-positive integer");
X    if (iszero(z1) || isunit(z2)) {
X        *res = _zero_;
X        return;
X    }
X    if (istwo(z2)) {
X        if (z1.v[0] & 0x1)
X            *res = _one_;
X        else
X            *res = _zero_;
X        return;
X    }
X
X    /*
X     * If the absolute value of the number is within the modulo range,
X     * then the result is just a copy or a subtraction.  Otherwise go
X     * ahead and negate and reduce the result.
X     */
X    sign = z1.sign;
X    z1.sign = 0;
X    cv = zrel(z1, z2);
X    if (cv == 0) {
X        *res = _zero_;
X        return;
X    }
X    if (cv < 0) {
X        if (sign)
X            zcopy(z1, res);
X        else
X            zsub(z2, z1, res);
X        return;
X    }
X    z1.sign = !sign;
X    zmod(z1, z2, res);
X}
X#endif
X
X
X/*
X * Calculate the number congruent to the given number whose absolute
X * value is minimal.  The number to be reduced can be negative or out of
X * modulo range.  The result will be within the range -int((modulus-1)/2)
X * to int(modulus/2) inclusive.  For example, for modulus 7, numbers are
X * reduced to the range [-3, 3], and for modulus 8, numbers are reduced to
X * the range [-3, 4].
X */
Xvoid
Xzminmod(z1, z2, res)
X    ZVALUE z1;        /* number to find minimum congruence of */
X    ZVALUE z2;        /* number to take mod with */
X    ZVALUE *res;        /* result */
X{
X    ZVALUE tmp1, tmp2;
X    int sign;
X    int cv;
X
X    if (iszero(z2) || isneg(z2))
X        error("Mod of non-positive integer");
X    if (iszero(z1) || isunit(z2)) {
X        *res = _zero_;
X        return;
X    }
X    if (istwo(z2)) {
X        if (isodd(z1))
X            *res = _one_;
X        else
X            *res = _zero_;
X        return;
X    }
X
X    /*
X     * Do a quick check to see if the number is very small compared
X     * to the modulus.  If so, then the result is obvious.
X     */
X    if (z1.len < z2.len - 1) {
X        zcopy(z1, res);
X        return;
X    }
X
X    /*
X     * Now make sure the input number is within the modulo range.
X     * If not, then reduce it to be within range and make the
X     * quick check again.
X     */
X    sign = z1.sign;
X    z1.sign = 0;
X    cv = zrel(z1, z2);
X    if (cv == 0) {
X        *res = _zero_;
X        return;
X    }
X    tmp1 = z1;
X    if (cv > 0) {
X        z1.sign = (BOOL)sign;
X        zmod(z1, z2, &tmp1);
X        if (tmp1.len < z2.len - 1) {
X            *res = tmp1;
X            return;
X        }
X        sign = 0;
X    }
X
X    /*
X     * Now calculate the difference of the modulus and the absolute
X     * value of the original number.  Compare the original number with
X     * the difference, and return the one with the smallest absolute
X     * value, with the correct sign.  If the two values are equal, then
X     * return the positive result.
X     */
X    zsub(z2, tmp1, &tmp2);
X    cv = zrel(tmp1, tmp2);
X    if (cv < 0) {
X        freeh(tmp2.v);
X        tmp1.sign = (BOOL)sign;
X        if (tmp1.v == z1.v)
X            zcopy(tmp1, res);
X        else
X            *res = tmp1;
X    } else {
X        if (cv)
X            tmp2.sign = !sign;
X        if (tmp1.v != z1.v)
X            freeh(tmp1.v);
X        *res = tmp2;
X    }
X}
X
X
X/*
X * Compare two numbers for equality modulo a third number.
X * The two numbers to be compared can be negative or out of modulo range.
X * Returns TRUE if the numbers are not congruent, and FALSE if they are
X * congruent.
X */
XBOOL
Xzcmpmod(z1, z2, z3)
X    ZVALUE z1;        /* first number to be compared */
X    ZVALUE z2;        /* second number to be compared */
X    ZVALUE z3;        /* modulus */
X{
X    ZVALUE tmp1, tmp2, tmp3;
X    FULL digit;
X    LEN len;
X    int cv;
X
X    if (isneg(z3) || iszero(z3))
X        error("Non-positive modulus in zcmpmod");
X    if (istwo(z3))
X        return (((z1.v[0] + z2.v[0]) & 0x1) != 0);
X
X    /*
X     * If the two numbers are equal, then their mods are equal.
X     */
X    if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
X        (z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0))
X            return FALSE;
X
X    /*
X     * If both numbers are negative, then we can make them positive.
X     */
X    if (isneg(z1) && isneg(z2)) {
X        z1.sign = 0;
X        z2.sign = 0;
X    }
X
X    /*
X     * For small negative numbers, make them positive before comparing.
X     * In any case, the resulting numbers are in tmp1 and tmp2.
X     */
X    tmp1 = z1;
X    tmp2 = z2;
X    len = z3.len;
X    digit = z3.v[len - 1];
X
X    if (isneg(z1) && ((z1.len < len) ||
X        ((z1.len == len) && (z1.v[z1.len - 1] < digit))))
X            zadd(z1, z3, &tmp1);
X
X    if (isneg(z2) && ((z2.len < len) ||
X        ((z2.len == len) && (z2.v[z2.len - 1] < digit))))
X            zadd(z2, z3, &tmp2);
X
X    /*
X     * Now compare the two numbers for equality.
X     * If they are equal we are all done.
X     */
X    if (zcmp(tmp1, tmp2) == 0) {
X        if (tmp1.v != z1.v)
X            freeh(tmp1.v);
X        if (tmp2.v != z2.v)
X            freeh(tmp2.v);
X        return FALSE;
X    }
X
X    /*
X     * They are not identical.  Now if both numbers are positive
X     * and less than the modulus, then they are definitely not equal.
X     */
X    if ((tmp1.sign == tmp2.sign) &&
X        ((tmp1.len < len) || (zrel(tmp1, z3) < 0)) &&
X        ((tmp2.len < len) || (zrel(tmp2, z3) < 0)))
X    {
X        if (tmp1.v != z1.v)
X            freeh(tmp1.v);
X        if (tmp2.v != z2.v)
X            freeh(tmp2.v);
X        return TRUE;
X    }
X
X    /*
X     * Either one of the numbers is negative or is large.
X     * So do the standard thing and subtract the two numbers.
X     * Then they are equal if the result is 0 (mod z3).
X     */
X    zsub(tmp1, tmp2, &tmp3);
X    if (tmp1.v != z1.v)
X        freeh(tmp1.v);
X    if (tmp2.v != z2.v)
X        freeh(tmp2.v);
X
X    /*
X     * Compare the result with the modulus to see if it is equal to
X     * or less than the modulus.  If so, we know the mod result.
X     */
X    tmp3.sign = 0;
X    cv = zrel(tmp3, z3);
X    if (cv == 0) {
X        freeh(tmp3.v);
X        return FALSE;
X    }
X    if (cv < 0) {
X        freeh(tmp3.v);
X        return TRUE;
X    }
X
X    /*
X     * We are forced to actually do the division.
X     * The numbers are congruent if the result is zero.
X     */
X    zmod(tmp3, z3, &tmp1);
X    freeh(tmp3.v);
X    if (iszero(tmp1)) {
X        freeh(tmp1.v);
X        return FALSE;
X    } else {
X        freeh(tmp1.v);
X        return TRUE;
X    }
X}
X
X
X/*
X * Compute the result of raising one number to a power modulo another number.
X * That is, this computes:  a^b (modulo c).
X * This calculates the result by examining the power POWBITS bits at a time,
X * using a small table of POWNUMS low powers to calculate powers for those bits,
X * and repeated squaring and multiplying by the partial powers to generate
X * the complete power.  If the power being raised to is high enough, then
X * this uses the REDC algorithm to avoid doing many divisions.  When using
X * REDC, multiple calls to this routine using the same modulus will be
X * slightly faster.
X */
Xvoid
Xzpowermod(z1, z2, z3, res)
X    ZVALUE z1, z2, z3, *res;
X{
X    HALF *hp;        /* pointer to current word of the power */
X    REDC *rp;        /* REDC information to be used */
X    ZVALUE *pp;        /* pointer to low power table */
X    ZVALUE ans, temp;    /* calculation values */
X    ZVALUE modpow;        /* current small power */
X    ZVALUE lowpowers[POWNUMS];    /* low powers */
X    int sign;        /* original sign of number */
X    int curshift;        /* shift value for word of power */
X    HALF curhalf;        /* current word of power */
X    unsigned int curpow;    /* current low power */
X    unsigned int curbit;    /* current bit of low power */
X    int i;
X
X    if (isneg(z3) || iszero(z3))
X        error("Non-positive modulus in zpowermod");
X    if (isneg(z2))
X        error("Negative power in zpowermod");
X
X    sign = z1.sign;
X    z1.sign = 0;
X
X    /*
X     * Check easy cases first.
X     */
X    if (iszero(z1) || isunit(z3)) {        /* 0^x or mod 1 */
X        *res = _zero_;
X        return;
X    }
X    if (istwo(z3)) {            /* mod 2 */
X        if (isodd(z1))
X            *res = _one_;
X        else
X            *res = _zero_;
X        return;
X    }
X    if (isunit(z1) && (!sign || iseven(z2))) {
X        /* 1^x or (-1)^(2x) */
X        *res = _one_;
X        return;
X    }
X
X    /*
X     * Normalize the number being raised to be non-negative and to lie
X     * within the modulo range.  Then check for zero or one specially.
X     */
X    zmod(z1, z3, &temp);
X    if (iszero(temp)) {
X        freeh(temp.v);
X        *res = _zero_;
X        return;
X    }
X    z1 = temp;
X    if (sign) {
X        zsub(z3, z1, &temp);
X        freeh(z1.v);
X        z1 = temp;
X    }
X    if (isunit(z1)) {
X        freeh(z1.v);
X        *res = _one_;
X        return;
X    }
X
X    /*
X     * If the modulus is odd, large enough, is not one less than an
X     * exact power of two, and if the power is large enough, then use
X     * the REDC algorithm.  The size where this is done is configurable.
X     */
X    if ((z2.len > 1) && (z3.len >= _pow2_) && isodd(z3)
X        && !zisallbits(z3))
X    {
X        if (powermodredc && zcmp(powermodredc->mod, z3)) {
X            zredcfree(powermodredc);
X            powermodredc = NULL;
X        }
X        if (powermodredc == NULL)
X            powermodredc = zredcalloc(z3);
X        rp = powermodredc;
X        zredcencode(rp, z1, &temp);
X        zredcpower(rp, temp, z2, &z1);
X        freeh(temp.v);
X        zredcdecode(rp, z1, res);
X        freeh(z1.v);
X        return;
X    }
X
X    /*
X     * Modulus or power is small enough to perform the power raising
X     * directly.  Initialize the table of powers.
X     */
X    for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
X        pp->len = 0;
X    lowpowers[0] = _one_;
X    lowpowers[1] = z1;
X    ans = _one_;
X
X    hp = &z2.v[z2.len - 1];
X    curhalf = *hp;
X    curshift = BASEB - POWBITS;
X    while (curshift && ((curhalf >> curshift) == 0))
X        curshift -= POWBITS;
X
X    /*
X     * Calculate the result by examining the power POWBITS bits at a time,
X     * and use the table of low powers at each iteration.
X     */
X    for (;;) {
X        curpow = (curhalf >> curshift) & (POWNUMS - 1);
X        pp = &lowpowers[curpow];
X
X        /*
X         * If the small power is not yet saved in the table, then
X         * calculate it and remember it in the table for future use.
X         */
X        if (pp->len == 0) {
X            if (curpow & 0x1)
X                zcopy(z1, &modpow);
X            else
X                modpow = _one_;
X
X            for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
X                pp = &lowpowers[curbit];
X                if (pp->len == 0) {
X                    zsquare(lowpowers[curbit/2], &temp);
X                    zmod(temp, z3, pp);
X                    freeh(temp.v);
X                }
X                if (curbit & curpow) {
X                    zmul(*pp, modpow, &temp);
X                    freeh(modpow.v);
X                    zmod(temp, z3, &modpow);
X                    freeh(temp.v);
X                }
X            }
X            pp = &lowpowers[curpow];
X            *pp = modpow;
X        }
X
X        /*
X         * If the power is nonzero, then accumulate the small power
X         * into the result.
X         */
X        if (curpow) {
X            zmul(ans, *pp, &temp);
X            freeh(ans.v);
X            zmod(temp, z3, &ans);
X            freeh(temp.v);
X        }
X
X        /*
X         * Select the next POWBITS bits of the power, if there is
X         * any more to generate.
X         */
X        curshift -= POWBITS;
X        if (curshift < 0) {
X            if (hp-- == z2.v)
X                break;
X            curhalf = *hp;
X            curshift = BASEB - POWBITS;
X        }
X
X        /*
X         * Square the result POWBITS times to make room for the next
X         * chunk of bits.
X         */
X        for (i = 0; i < POWBITS; i++) {
X            zsquare(ans, &temp);
X            freeh(ans.v);
X            zmod(temp, z3, &ans);
X            freeh(temp.v);
X        }
X    }
X
X    for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) {
X        if (pp->len)
X            freeh(pp->v);
X    }
X    *res = ans;
X}
X
X
X/*
X * Initialize the REDC algorithm for a particular modulus,
X * returning a pointer to a structure that is used for other
X * REDC calls.  An error is generated if the structure cannot
X * be allocated.  The modulus must be odd and positive.
X */
XREDC *
Xzredcalloc(z1)
X    ZVALUE z1;        /* modulus to initialize for */
X{
X    REDC *rp;        /* REDC information */
X    ZVALUE tmp;
X    long bit;
X
X    if (iseven(z1) || isneg(z1))
X        error("REDC requires positive odd modulus");
X
X    rp = (REDC *) malloc(sizeof(REDC));
X    if (rp == NULL)
X        error("Cannot allocate REDC structure");
X
X    /*
X     * Round up the binary modulus to the next power of two
X     * which is at a word boundary.  Then the shift and modulo
X     * operations mod the binary modulus can be done very cheaply.
X     * Calculate the REDC format for the number 1 for future use.
X     */
X    bit = zhighbit(z1) + 1;
X    if (bit % BASEB)
X        bit += (BASEB - (bit % BASEB));
X    zcopy(z1, &rp->mod);
X    zbitvalue(bit, &tmp);
X    z1.sign = 1;
X    (void) zmodinv(z1, tmp, &rp->inv);
X    zmod(tmp, rp->mod, &rp->one);
X    freeh(tmp.v);
X    rp->len = bit / BASEB;
X    return rp;
X}
X
X
X/*
X * Free any numbers associated with the specified REDC structure,
X * and then the REDC structure itself.
X */
Xvoid
Xzredcfree(rp)
X    REDC *rp;        /* REDC information to be cleared */
X{
X    freeh(rp->mod.v);
X    freeh(rp->inv.v);
X    freeh(rp->one.v);
X    free(rp);
X}
X
X
X/*
X * Convert a normal number into the specified REDC format.
X * The number to be converted can be negative or out of modulo range.
X * The resulting number can be used for multiplying, adding, subtracting,
X * or comparing with any other such converted numbers, as if the numbers
X * were being calculated modulo the number which initialized the REDC
X * information.  When the final value is unconverted, the result is the
X * same as if the usual operations were done with the original numbers.
X */
Xvoid
Xzredcencode(rp, z1, res)
X    REDC *rp;        /* REDC information */
X    ZVALUE z1;        /* number to be converted */
X    ZVALUE *res;        /* returned converted number */
X{
X    ZVALUE tmp1, tmp2;
X
X    /*
X     * Handle the cases 0, 1, -1, and 2 specially since these are
X     * easy to calculate.  Zero transforms to zero, and the others
X     * can be obtained from the precomputed REDC format for 1 since
X     * addition and subtraction act normally for REDC format numbers.
X     */
X    if (iszero(z1)) {
X        *res = _zero_;
X        return;
X    }
X    if (isone(z1)) {
X        zcopy(rp->one, res);
X        return;
X    }
X    if (isunit(z1)) {
X        zsub(rp->mod, rp->one, res);
X        return;
X    }
X    if (istwo(z1)) {
X        zadd(rp->one, rp->one, &tmp1);
X        if (zrel(tmp1, rp->mod) < 0) {
X            *res = tmp1;
X            return;
X        }
X        zsub(tmp1, rp->mod, res);
X        freeh(tmp1.v);
X        return;
X    }
X
X    /*
X     * Not a trivial number to convert, so do the full transformation.
X     * Convert negative numbers to positive numbers before converting.
X     */
X    tmp1.len = 0;
X    if (isneg(z1)) {
X        zmod(z1, rp->mod, &tmp1);
X        z1 = tmp1;
X    }
X    zshift(z1, rp->len * BASEB, &tmp2);
X    if (tmp1.len)
X        freeh(tmp1.v);
X    zmod(tmp2, rp->mod, res);
X    freeh(tmp2.v);
X}
X
X
X/*
X * The REDC algorithm used to convert numbers out of REDC format and also
X * used after multiplication of two REDC numbers.  Using this routine
X * avoids any divides, replacing the divide by two multiplications.
X * If the numbers are very large, then these two multiplies will be
X * quicker than the divide, since dividing is harder than multiplying.
X */
Xvoid
Xzredcdecode(rp, z1, res)
X    REDC *rp;        /* REDC information */
X    ZVALUE z1;        /* number to be transformed */
X    ZVALUE *res;        /* returned transformed number */
X{
X    ZVALUE tmp1, tmp2;    /* temporaries */
X    HALF *hp;        /* saved pointer to tmp2 value */
X
X    if (isneg(z1))
X        error("Negative number for zredc");
X
X    /*
X     * Check first for the special values for 0 and 1 that are easy.
X     */
X    if (iszero(z1)) {
X        *res = _zero_;
X        return;
X    }
X    if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
X        (zcmp(z1, rp->one) == 0)) {
X            *res = _one_;
X            return;
X    }
X
X    /*
X     * First calculate the following:
X     *     tmp2 = ((z1 % 2^bitnum) * inv) % 2^bitnum.
X     * The mod operations can be done with no work since the bit
X     * number was selected as a multiple of the word size.  Just
X     * reduce the sizes of the numbers as required.
X     */
X    tmp1 = z1;
X    if (tmp1.len > rp->len)
X        tmp1.len = rp->len;
X    zmul(tmp1, rp->inv, &tmp2);
X    if (tmp2.len > rp->len)
X        tmp2.len = rp->len;
X
X    /*
X     * Next calculate the following:
X     *    res = (z1 + tmp2 * modulus) / 2^bitnum
X     * The division by a power of 2 is always exact, and requires no
X     * work.  Just adjust the address and length of the number to do
X     * the divide, but save the original pointer for freeing later.
X     */
X    zmul(tmp2, rp->mod, &tmp1);
X    freeh(tmp2.v);
X    zadd(z1, tmp1, &tmp2);
X    freeh(tmp1.v);
X    hp = tmp2.v;
X    if (tmp2.len <= rp->len) {
X        freeh(hp);
X        *res = _zero_;
X        return;
X    }
X    tmp2.v += rp->len;
X    tmp2.len -= rp->len;
X
X    /*
X     * Finally do a final modulo by a simple subtraction if necessary.
X     * This is all that is needed because the previous calculation is
X     * guaranteed to always be less than twice the modulus.
X     */
X    if (zrel(tmp2, rp->mod) < 0)
X        zcopy(tmp2, res);
X    else
X        zsub(tmp2, rp->mod, res);
X    freeh(hp);
X}
X
X
X/*
X * Multiply two numbers in REDC format together producing a result also
X * in REDC format.  If the result is converted back to a normal number,
X * then the result is the same as the modulo'd multiplication of the
X * original numbers before they were converted to REDC format.  This
X * calculation is done in one of two ways, depending on the size of the
X * modulus.  For large numbers, the REDC definition is used directly
X * which involves three multiplies overall.  For small numbers, a
X * complicated routine is used which does the indicated multiplication
X * and the REDC algorithm at the same time to produce the result.
X */
Xvoid
Xzredcmul(rp, z1, z2, res)
X    REDC *rp;        /* REDC information */
X    ZVALUE z1;        /* first REDC number to be multiplied */
X    ZVALUE z2;        /* second REDC number to be multiplied */
X    ZVALUE *res;        /* resulting REDC number */
X{
X    FULL mulb;
X    FULL muln;
X    HALF *h1;
X    HALF *h2;
X    HALF *h3;
X    HALF *hd;
X    HALF Ninv;
X    HALF topdigit;
X    LEN modlen;
X    LEN len;
X    LEN len2;
X    SIUNION sival1;
X    SIUNION sival2;
X    SIUNION sival3;
X    SIUNION carry;
X    ZVALUE tmp;
X
X    if (isneg(z1) || (z1.len > rp->mod.len) ||
X        isneg(z2) || (z2.len > rp->mod.len))
X            error("Negative or too large number in zredcmul");
X
X    /*
X     * Check for special values which we easily know the answer.
X     */
X    if (iszero(z1) || iszero(z2)) {
X        *res = _zero_;
X        return;
X    }
X
X    if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
X        (zcmp(z1, rp->one) == 0)) {
X            zcopy(z2, res);
X            return;
X    }
X
X    if ((z2.len == rp->one.len) && (z2.v[0] == rp->one.v[0]) &&
X        (zcmp(z2, rp->one) == 0)) {
X            zcopy(z1, res);
X            return;
X    }
X
X    /*
X     * If the size of the modulus is large, then just do the multiply,
X     * followed by the two multiplies contained in the REDC routine.
X     * This will be quicker than directly doing the REDC calculation
X     * because of the O(N^1.585) speed of the multiplies.  The size
X     * of the number which this is done is configurable.
X     */
X    if (rp->mod.len >= _redc2_) {
X        zmul(z1, z2, &tmp);
X        zredcdecode(rp, tmp, res);
X        freeh(tmp.v);
X        return;
X    }
X
X    /*
X     * The number is small enough to calculate by doing the O(N^2) REDC
X     * algorithm directly.  This algorithm performs the multiplication and
X     * the reduction at the same time.  Notice the obscure facts that
X     * only the lowest word of the inverse value is used, and that
X     * there is no shifting of the partial products as there is in a
X     * normal multiply.
X     */
X    modlen = rp->mod.len;
X    Ninv = rp->inv.v[0];
X
X    /*
X     * Allocate the result and clear it.
X     * The size of the result will be equal to or smaller than
X     * the modulus size.
X     */
X    res->sign = 0;
X    res->len = modlen;
X    res->v = alloc(modlen);
X
X    hd = res->v;
X    len = modlen;
X    while (len--)
X        *hd++ = 0;
X
X    /*
X     * Do this outermost loop over all the digits of z1.
X     */
X    h1 = z1.v;
X    len = z1.len;
X    while (len--) {
X        /*
X         * Start off with the next digit of z1, the first
X         * digit of z2, and the first digit of the modulus.
X         */
X        mulb = (FULL) *h1++;
X        h2 = z2.v;
X        h3 = rp->mod.v;
X        hd = res->v;
X        sival1.ivalue = mulb * ((FULL) *h2++) + ((FULL) *hd++);
X        muln = ((HALF) (sival1.silow * Ninv));
X        sival2.ivalue = muln * ((FULL) *h3++);
X        sival3.ivalue = ((FULL) sival1.silow) + ((FULL) sival2.silow);
X        carry.ivalue = ((FULL) sival1.sihigh) + ((FULL) sival2.sihigh)
X            + ((FULL) sival3.sihigh);
X
X        /*
X         * Do this innermost loop for each digit of z2, except
X         * for the first digit which was just done above.
X         */
X        len2 = z2.len;
X        while (--len2 > 0) {
X            sival1.ivalue = mulb * ((FULL) *h2++);
X            sival2.ivalue = muln * ((FULL) *h3++);
X            sival3.ivalue = ((FULL) sival1.silow)
X                + ((FULL) sival2.silow)
X                + ((FULL) *hd)
X                + ((FULL) carry.silow);
X            carry.ivalue = ((FULL) sival1.sihigh)
X                + ((FULL) sival2.sihigh)
X                + ((FULL) sival3.sihigh)
X                + ((FULL) carry.sihigh);
X
X            hd[-1] = sival3.silow;
X            hd++;
X        }
X
X        /*
X         * Now continue the loop as necessary so the total number
X         * of interations is equal to the size of the modulus.
X         * This acts as if the innermost loop was repeated for
X         * high digits of z2 that are zero.
X         */
X        len2 = modlen - z2.len;
X        while (len2--) {
X            sival2.ivalue = muln * ((FULL) *h3++);
X            sival3.ivalue = ((FULL) sival2.silow)
X                + ((FULL) *hd)
X                + ((FULL) carry.silow);
X            carry.ivalue = ((FULL) sival2.sihigh)
X                + ((FULL) sival3.sihigh)
X                + ((FULL) carry.sihigh);
X
X            hd[-1] = sival3.silow;
X            hd++;
X        }
X
X        res->v[modlen - 1] = carry.silow;
X        topdigit = carry.sihigh;
X    }
X
X    /*
X     * Now continue the loop as necessary so the total number
X     * of interations is equal to the size of the modulus.
X     * This acts as if the outermost loop was repeated for high
X     * digits of z1 that are zero.
X     */
X    len = modlen - z1.len;
X    while (len--) {
X        /*
X         * Start off with the first digit of the modulus.
X         */
X        h3 = rp->mod.v;
X        hd = res->v;
X        muln = ((HALF) (*hd * Ninv));
X        sival2.ivalue = muln * ((FULL) *h3++);
X        sival3.ivalue = ((FULL) *hd++) + ((FULL) sival2.silow);
X        carry.ivalue = ((FULL) sival2.sihigh) + ((FULL) sival3.sihigh);
X
X        /*
X         * Do this innermost loop for each digit of the modulus,
X         * except for the first digit which was just done above.
X         */
X        len2 = modlen;
X        while (--len2 > 0) {
X            sival2.ivalue = muln * ((FULL) *h3++);
X            sival3.ivalue = ((FULL) sival2.silow)
X                + ((FULL) *hd)
X                + ((FULL) carry.silow);
X            carry.ivalue = ((FULL) sival2.sihigh)
X                + ((FULL) sival3.sihigh)
X                + ((FULL) carry.sihigh);
X
X            hd[-1] = sival3.silow;
X            hd++;
X        }
X        res->v[modlen - 1] = carry.silow;
X        topdigit = carry.sihigh;
X    }
X
X    /*
X     * Determine the true size of the result, taking the top digit of
X     * the current result into account.  The top digit is not stored in
X     * the number because it is temporary and would become zero anyway
X     * after the final subtraction is done.
X     */
X    if (topdigit == 0) {
X        len = modlen;
X        hd = &res->v[len - 1];
X        while ((*hd == 0) && (len > 1)) {
X            hd--;
X            len--;
X        }
X        res->len = len;
X    }
X
X    /*
X     * Compare the result with the modulus.
X     * If it is less than the modulus, then the calculation is complete.
X     */
X    if ((topdigit == 0) && ((len < modlen)
X        || (res->v[len - 1] < rp->mod.v[len - 1])
X        || (zrel(*res, rp->mod) < 0)))
X            return;
X
X    /*
X     * Do a subtraction to reduce the result to a value less than
X     * the modulus.  The REDC algorithm guarantees that a single subtract
X     * is all that is needed.  Ignore any borrowing from the possible
X     * highest word of the current result because that would affect
X     * only the top digit value that was not stored and would become
X     * zero anyway.
X     */
X    carry.ivalue = 0;
X    h1 = rp->mod.v;
X    hd = res->v;
X    len = modlen;
X    while (len--) {
X        carry.ivalue = BASE1 - ((FULL) *hd) + ((FULL) *h1++)
X            + ((FULL) carry.silow);
X        *hd++ = BASE1 - carry.silow;
X        carry.silow = carry.sihigh;
X    }
X
X    /*
X     * Now finally recompute the size of the result.
X     */
X    len = modlen;
X    hd = &res->v[len - 1];
X    while ((*hd == 0) && (len > 1)) {
X        hd--;
X        len--;
X    }
X    res->len = len;
X}
X
X
X/*
X * Square a number in REDC format producing a result also in REDC format.
X */
Xvoid
Xzredcsquare(rp, z1, res)
X    REDC *rp;        /* REDC information */
X    ZVALUE z1;        /* REDC number to be squared */
X    ZVALUE *res;        /* resulting REDC number */
X{
X    ZVALUE tmp;
X
X    if (isneg(z1))
X        error("Negative number in zredcsquare");
X    if (iszero(z1)) {
X        *res = _zero_;
X        return;
X    }
X    if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
X        (zcmp(z1, rp->one) == 0)) {
X            zcopy(z1, res);
X            return;
X    }
X
X    /*
X     * If the modulus is small enough, then call the multiply
X     * routine to produce the result.  Otherwise call the O(N^1.585)
X     * routines to get the answer.
X     */
X    if (rp->mod.len < _redc2_) {
X        zredcmul(rp, z1, z1, res);
X        return;
X    }
X    zsquare(z1, &tmp);
X    zredcdecode(rp, tmp, res);
X    freeh(tmp.v);
X}
X
X
X/*
X * Compute the result of raising a REDC format number to a power.
X * The result is within the range 0 to the modulus - 1.
X * This calculates the result by examining the power POWBITS bits at a time,
X * using a small table of POWNUMS low powers to calculate powers for those bits,
X * and repeated squaring and multiplying by the partial powers to generate
X * the complete power.
X */
Xvoid
Xzredcpower(rp, z1, z2, res)
X    REDC *rp;        /* REDC information */
X    ZVALUE z1;        /* REDC number to be raised */
X    ZVALUE z2;        /* normal number to raise number to */
X    ZVALUE *res;        /* result */
X{
X    HALF *hp;        /* pointer to current word of the power */
X    ZVALUE *pp;        /* pointer to low power table */
X    ZVALUE ans, temp;    /* calculation values */
X    ZVALUE modpow;        /* current small power */
X    ZVALUE lowpowers[POWNUMS];    /* low powers */
X    int curshift;        /* shift value for word of power */
X    HALF curhalf;        /* current word of power */
X    unsigned int curpow;    /* current low power */
X    unsigned int curbit;    /* current bit of low power */
X    int i;
X
X    if (isneg(z1))
X        error("Negative number in zredcpower");
X    if (isneg(z2))
X        error("Negative power in zredcpower");
X
X    /*
X     * Check for zero or the REDC format for one.
X     */
X    if (iszero(z1) || isunit(rp->mod)) {
X        *res = _zero_;
X        return;
X    }
X    if (zcmp(z1, rp->one) == 0) {
X        zcopy(rp->one, res);
X        return;
X    }
X
X    /*
X     * See if the number being raised is the REDC format for -1.
X     * If so, then the answer is the REDC format for one or minus one.
X     * To do this check, calculate the REDC format for -1.
X     */
X    if (((HALF)(z1.v[0] + rp->one.v[0])) == rp->mod.v[0]) {
X        zsub(rp->mod, rp->one, &temp);
X        if (zcmp(z1, temp) == 0) {
X            if (isodd(z2)) {
X                *res = temp;
X                return;
X            }
X            freeh(temp.v);
X            zcopy(rp->one, res);
X            return;
X        }
X        freeh(temp.v);
X    }
X
X    for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
X        pp->len = 0;
X    zcopy(rp->one, &lowpowers[0]);
X    zcopy(z1, &lowpowers[1]);
X    zcopy(rp->one, &ans);
X
X    hp = &z2.v[z2.len - 1];
X    curhalf = *hp;
X    curshift = BASEB - POWBITS;
X    while (curshift && ((curhalf >> curshift) == 0))
X        curshift -= POWBITS;
X
X    /*
X     * Calculate the result by examining the power POWBITS bits at a time,
X     * and use the table of low powers at each iteration.
X     */
X    for (;;) {
X        curpow = (curhalf >> curshift) & (POWNUMS - 1);
X        pp = &lowpowers[curpow];
X
X        /*
X         * If the small power is not yet saved in the table, then
X         * calculate it and remember it in the table for future use.
X         */
X        if (pp->len == 0) {
X            if (curpow & 0x1)
X                zcopy(z1, &modpow);
X            else
X                zcopy(rp->one, &modpow);
X
X            for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
X                pp = &lowpowers[curbit];
X                if (pp->len == 0)
X                    zredcsquare(rp, lowpowers[curbit/2],
X                        pp);
X                if (curbit & curpow) {
X                    zredcmul(rp, *pp, modpow, &temp);
X                    freeh(modpow.v);
X                    modpow = temp;
X                }
X            }
X            pp = &lowpowers[curpow];
X            *pp = modpow;
X        }
X
X        /*
X         * If the power is nonzero, then accumulate the small power
X         * into the result.
X         */
X        if (curpow) {
X            zredcmul(rp, ans, *pp, &temp);
X            freeh(ans.v);
X            ans = temp;
X        }
X
X        /*
X         * Select the next POWBITS bits of the power, if there is
X         * any more to generate.
X         */
X        curshift -= POWBITS;
X        if (curshift < 0) {
X            if (hp-- == z2.v)
X                break;
X            curhalf = *hp;
X            curshift = BASEB - POWBITS;
X        }
X
X        /*
X         * Square the result POWBITS times to make room for the next
X         * chunk of bits.
X         */
X        for (i = 0; i < POWBITS; i++) {
X            zredcsquare(rp, ans, &temp);
X            freeh(ans.v);
X            ans = temp;
X        }
X    }
X
X    for (pp = lowpowers; pp < &lowpowers[POWNUMS]; pp++) {
X        if (pp->len)
X            freeh(pp->v);
X    }
X    *res = ans;
X}
X
X/* END CODE */
END_OF_FILE
if test 32468 -ne `wc -c <'zmod.c'`; then
    echo shar: \"'zmod.c'\" unpacked with wrong size!
fi
# end of 'zmod.c'
fi
echo shar: End of archive 16 \(of 21\).
cp /dev/null ark16isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 21 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0