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

Text File | 1992-05-09 | 36.8 KB | 1,724 lines

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

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

#! /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 18 (of 21)."
# Contents:  zfunc.c
# Wrapped by dbell@elm on Tue Feb 25 15:21:16 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'zfunc.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'zfunc.c'\"
else
echo shar: Extracting \"'zfunc.c'\" \(34345 characters\)
sed "s/^X//" >'zfunc.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 * Extended precision integral arithmetic non-primitive routines
X */
X
X#include "math.h"
X
Xstatic ZVALUE primeprod;        /* product of primes under 100 */
XZVALUE _tenpowers_[32];            /* table of 10^2^n */
X
X#if 0
Xstatic char *abortmsg = "Calculation aborted";
Xstatic char *memmsg = "Not enough memory";
X#endif
X
X
X/*
X * Compute the factorial of a number.
X */
Xvoid
Xzfact(z, dest)
X    ZVALUE z, *dest;
X{
X    long ptwo;        /* count of powers of two */
X    long n;            /* current multiplication value */
X    long m;            /* reduced multiplication value */
X    long mul;        /* collected value to multiply by */
X    ZVALUE res, temp;
X
X    if (isneg(z))
X        error("Negative argument for factorial");
X    if (isbig(z))
X        error("Very large factorial");
X    n = (istiny(z) ? z1tol(z) : z2tol(z));
X    ptwo = 0;
X    mul = 1;
X    res = _one_;
X    /*
X     * Multiply numbers together, but squeeze out all powers of two.
X     * We will put them back in at the end.  Also collect multiple
X     * numbers together until there is a risk of overflow.
X     */
X    for (; n > 1; n--) {
X        for (m = n; ((m & 0x1) == 0); m >>= 1)
X            ptwo++;
X        mul *= m;
X        if (mul < BASE1/2)
X            continue;
X        zmuli(res, mul, &temp);
X        freeh(res.v);
X        res = temp;
X        mul = 1;
X    }
X    /*
X     * Multiply by the remaining value, then scale result by
X     * the proper power of two.
X     */
X    if (mul > 1) {
X        zmuli(res, mul, &temp);
X        freeh(res.v);
X        res = temp;
X    }
X    zshift(res, ptwo, &temp);
X    freeh(res.v);
X    *dest = temp;
X}
X
X
X/*
X * Compute the product of the primes up to the specified number.
X */
Xvoid
Xzpfact(z, dest)
X    ZVALUE z, *dest;
X{
X    long n;            /* limiting number to multiply by */
X    long p;            /* current prime */
X    long i;            /* test value */
X    long mul;        /* collected value to multiply by */
X    ZVALUE res, temp;
X
X    if (isneg(z))
X        error("Negative argument for factorial");
X    if (isbig(z))
X        error("Very large factorial");
X    n = (istiny(z) ? z1tol(z) : z2tol(z));
X    /*
X     * Multiply by the primes in order, collecting multiple numbers
X     * together until there is a change of overflow.
X     */
X    mul = 1 + (n > 1);
X    res = _one_;
X    for (p = 3; p <= n; p += 2) {
X        for (i = 3; (i * i) <= p; i += 2) {
X            if ((p % i) == 0)
X                goto next;
X        }
X        mul *= p;
X        if (mul < BASE1/2)
X            continue;
X        zmuli(res, mul, &temp);
X        freeh(res.v);
X        res = temp;
X        mul = 1;
Xnext: ;
X    }
X    /*
X     * Multiply by the final value if any.
X     */
X    if (mul > 1) {
X        zmuli(res, mul, &temp);
X        freeh(res.v);
X        res = temp;
X    }
X    *dest = res;
X}
X
X
X/*
X * Compute the least common multiple of all the numbers up to the
X * specified number.
X */
Xvoid
Xzlcmfact(z, dest)
X    ZVALUE z, *dest;
X{
X    long n;            /* limiting number to multiply by */
X    long p;            /* current prime */
X    long pp;        /* power of prime */
X    long i;            /* test value */
X    ZVALUE res, temp;
X
X    if (isneg(z) || iszero(z))
X        error("Non-positive argument for lcmfact");
X    if (isbig(z))
X        error("Very large lcmfact");
X    n = (istiny(z) ? z1tol(z) : z2tol(z));
X    /*
X     * Multiply by powers of the necessary odd primes in order.
X     * The power for each prime is the highest one which is not
X     * more than the specified number.
X     */
X    res = _one_;
X    for (p = 3; p <= n; p += 2) {
X        for (i = 3; (i * i) <= p; i += 2) {
X            if ((p % i) == 0)
X                goto next;
X        }
X        i = p;
X        while (i <= n) {
X            pp = i;
X            i *= p;
X        }
X        zmuli(res, pp, &temp);
X        freeh(res.v);
X        res = temp;
Xnext: ;
X    }
X    /*
X     * Finish by scaling by the necessary power of two.
X     */
X    zshift(res, zhighbit(z), dest);
X    freeh(res.v);
X}
X
X
X/*
X * Compute the permuation function  M! / (M - N)!.
X */
Xvoid
Xzperm(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    long count;
X    ZVALUE cur, tmp, ans;
X
X    if (isneg(z1) || isneg(z2))
X        error("Negative argument for permutation");
X    if (zrel(z1, z2) < 0)
X        error("Second arg larger than first in permutation");
X    if (isbig(z2))
X        error("Very large permutation");
X    count = (istiny(z2) ? z1tol(z2) : z2tol(z2));
X    zcopy(z1, &ans);
X    zsub(z1, _one_, &cur);
X    while (--count > 0) {
X        zmul(ans, cur, &tmp);
X        freeh(ans.v);
X        ans = tmp;
X        zsub(cur, _one_, &tmp);
X        freeh(cur.v);
X        cur = tmp;
X    }
X    freeh(cur.v);
X    *res = ans;
X}
X
X
X/*
X * Compute the combinatorial function  M! / ( N! * (M - N)! ).
X */
Xvoid
Xzcomb(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE ans;
X    ZVALUE mul, div, temp;
X    FULL count, i;
X    HALF dh[2];
X
X    if (isneg(z1) || isneg(z2))
X        error("Negative argument for combinatorial");
X    zsub(z1, z2, &temp);
X    if (isneg(temp)) {
X        freeh(temp.v);
X        error("Second arg larger than first for combinatorial");
X    }
X    if (isbig(z2) && isbig(temp)) {
X        freeh(temp.v);
X        error("Very large combinatorial");
X    }
X    count = (istiny(z2) ? z1tol(z2) : z2tol(z2));
X    i = (istiny(temp) ? z1tol(temp) : z2tol(temp));
X    if (isbig(z2) || (!isbig(temp) && (i < count)))
X        count = i;
X    freeh(temp.v);
X    mul = z1;
X    div.sign = 0;
X    div.v = dh;
X    ans = _one_;
X    for (i = 1; i <= count; i++) {
X        dh[0] = i & BASE1;
X        dh[1] = i / BASE;
X        div.len = 1 + (dh[1] != 0);
X        zmul(ans, mul, &temp);
X        freeh(ans.v);
X        zquo(temp, div, &ans);
X        freeh(temp.v);
X        zsub(mul, _one_, &temp);
X        if (mul.v != z1.v)
X            freeh(mul.v);
X        mul = temp;
X    }
X    if (mul.v != z1.v)
X        freeh(mul.v);
X    *res = ans;
X}
X
X
X/*
X * Perform a probabilistic primality test (algorithm P in Knuth).
X * Returns FALSE if definitely not prime, or TRUE if probably prime.
X * Count determines how many times to check for primality.
X * The chance of a non-prime passing this test is less than (1/4)^count.
X * For example, a count of 100 fails for only 1 in 10^60 numbers.
X */
XBOOL
Xzprimetest(z, count)
X    ZVALUE z;        /* number to test for primeness */
X    long count;
X{
X    long ij, ik, ix;
X    ZVALUE zm1, z1, z2, z3, ztmp;
X    HALF val[2];
X
X    z.sign = 0;
X    if (iseven(z))        /* if even, not prime if not 2 */
X        return (istwo(z) != 0);
X    /*
X     * See if the number is small, and is either a small prime,
X     * or is divisible by a small prime.
X     */
X    if (istiny(z) && (*z.v <= (HALF)(101*101-1))) {
X        ix = *z.v;
X        for (ik = 3; (ik <= 97) && ((ik * ik) <= ix); ik += 2)
X            if ((ix % ik) == 0)
X                return FALSE;
X        return TRUE;
X    }
X    /*
X     * See if the number is divisible by one of the primes 3, 5,
X     * 7, 11, or 13.  This is a very easy check.
X     */
X    ij = zmodi(z, 15015L);
X    if (!(ij % 3) || !(ij % 5) || !(ij % 7) || !(ij % 11) || !(ij % 13))
X        return FALSE;
X    /*
X     * Check the gcd of the number and the product of more of the first
X     * few odd primes.  We must build the prime product on the first call.
X     */
X    ztmp.sign = 0;
X    ztmp.len = 1;
X    ztmp.v = val;
X    if (primeprod.len == 0) {
X        val[0] = 101;
X        zpfact(ztmp, &primeprod);
X    }
X    zgcd(z, primeprod, &z1);
X    if (!isunit(z1)) {
X        freeh(z1.v);
X        return FALSE;
X    }
X    freeh(z1.v);
X    /*
X     * Not divisible by a small prime, so onward with the real test.
X     * Make sure the count is limited by the number of odd numbers between
X     * three and the number being tested.
X     */
X    ix = ((istiny(z) ? z1tol(z) : z2tol(z) - 3) / 2);
X    if (count > ix) count = ix;
X    zsub(z, _one_, &zm1);
X    ik = zlowbit(zm1);
X    zshift(zm1, -ik, &z1);
X    /*
X     * Loop over various "random" numbers, testing each one.
X     * These numbers are the odd numbers starting from three.
X     */
X    for (ix = 0; ix < count; ix++) {
X        val[0] = (ix * 2) + 3;
X        ij = 0;
X        zpowermod(ztmp, z1, z, &z3);
X        for (;;) {
X            if (isone(z3)) {
X                if (ij)    /* number is definitely not prime */
X                    goto notprime;
X                break;
X            }
X            if (zcmp(z3, zm1) == 0)
X                break;
X            if (++ij >= ik)
X                goto notprime;    /* number is definitely not prime */
X            zsquare(z3, &z2);
X            freeh(z3.v);
X            zmod(z2, z, &z3);
X            freeh(z2.v);
X        }
X        freeh(z3.v);
X    }
X    freeh(zm1.v);
X    freeh(z1.v);
X    return TRUE;    /* number might be prime */
X
Xnotprime:
X    freeh(z3.v);
X    freeh(zm1.v);
X    freeh(z1.v);
X    return FALSE;
X}
X
X
X/*
X * Compute the Jacobi function (p / q) for odd q.
X * If q is prime then the result is:
X *    1 if p == x^2 (mod q) for some x.
X *    -1 otherwise.
X * If q is not prime, then the result is not meaningful if it is 1.
X * This function returns 0 if q is even or q < 0.
X */
XFLAG
Xzjacobi(z1, z2)
X    ZVALUE z1, z2;
X{
X    ZVALUE p, q, tmp;
X    long lowbit;
X    int val;
X
X    if (iseven(z2) || isneg(z2))
X        return 0;
X    val = 1;
X    if (iszero(z1) || isone(z1))
X        return val;
X    if (isunit(z1)) {
X        if ((*z2.v - 1) & 0x2)
X            val = -val;
X        return val;
X    }
X    zcopy(z1, &p);
X    zcopy(z2, &q);
X    for (;;) {
X        zmod(p, q, &tmp);
X        freeh(p.v);
X        p = tmp;
X        if (iszero(p)) {
X            freeh(p.v);
X            p = _one_;
X        }
X        if (iseven(p)) {
X            lowbit = zlowbit(p);
X            zshift(p, -lowbit, &tmp);
X            freeh(p.v);
X            p = tmp;
X            if ((lowbit & 1) && (((*q.v & 0x7) == 3) || ((*q.v & 0x7) == 5)))
X                val = -val;
X        }
X        if (isunit(p)) {
X            freeh(p.v);
X            freeh(q.v);
X            return val;
X        }
X        if ((*p.v & *q.v & 0x3) == 3)
X            val = -val;
X        tmp = q;
X        q = p;
X        p = tmp;
X    }
X}
X
X
X/*
X * Return the Fibonacci number F(n).
X * This is evaluated by recursively using the formulas:
X *    F(2N+1) = F(N+1)^2 + F(N)^2
X * and
X *    F(2N) = F(N+1)^2 - F(N-1)^2
X */
Xvoid
Xzfib(z, res)
X    ZVALUE z, *res;
X{
X    unsigned long i;
X    long n;
X    int sign;
X    ZVALUE fnm1, fn, fnp1;        /* consecutive fibonacci values */
X    ZVALUE t1, t2, t3;
X
X    if (isbig(z))
X        error("Very large Fibonacci number");
X    n = (istiny(z) ? z1tol(z) : z2tol(z));
X    if (n == 0) {
X        *res = _zero_;
X        return;
X    }
X    sign = z.sign && ((n & 0x1) == 0);
X    if (n <= 2) {
X        *res = _one_;
X        res->sign = (BOOL)sign;
X        return;
X    }
X    i = TOPFULL;
X    while ((i & n) == 0)
X        i >>= 1L;
X    i >>= 1L;
X    fnm1 = _zero_;
X    fn = _one_;
X    fnp1 = _one_;
X    while (i) {
X        zsquare(fnm1, &t1);
X        zsquare(fn, &t2);
X        zsquare(fnp1, &t3);
X        freeh(fnm1.v);
X        freeh(fn.v);
X        freeh(fnp1.v);
X        zadd(t2, t3, &fnp1);
X        zsub(t3, t1, &fn);
X        freeh(t1.v);
X        freeh(t2.v);
X        freeh(t3.v);
X        if (i & n) {
X            fnm1 = fn;
X            fn = fnp1;
X            zadd(fnm1, fn, &fnp1);
X        } else
X            zsub(fnp1, fn, &fnm1);
X        i >>= 1L;
X    }
X    freeh(fnm1.v);
X    freeh(fnp1.v);
X    *res = fn;
X    res->sign = (BOOL)sign;
X}
X
X
X/*
X * Compute the result of raising one number to the power of another
X * The second number is assumed to be non-negative.
X * It cannot be too large except for trivial cases.
X */
Xvoid
Xzpowi(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    int sign;        /* final sign of number */
X    unsigned long power;    /* power to raise to */
X    unsigned long bit;    /* current bit value */
X    long twos;        /* count of times 2 is in result */
X    ZVALUE ans, temp;
X
X    sign = (z1.sign && isodd(z2));
X    z1.sign = 0;
X    z2.sign = 0;
X    if (iszero(z2)) {    /* number raised to power 0 */
X        if (iszero(z1))
X            error("Zero raised to zero power");
X        *res = _one_;
X        return;
X    }
X    if (isleone(z1)) {    /* 0, 1, or -1 raised to a power */
X        ans = _one_;
X        ans.sign = (BOOL)sign;
X        if (*z1.v == 0)
X            ans = _zero_;
X        *res = ans;
X        return;
X    }
X    if (isbig(z2))
X        error("Raising to very large power");
X    power = (istiny(z2) ? z1tol(z2) : z2tol(z2));
X    if (istwo(z1)) {    /* two raised to a power */
X        zbitvalue((long) power, res);
X        return;
X    }
X    /*
X     * See if this is a power of ten
X     */
X    if (istiny(z1) && (*z1.v == 10)) {
X        ztenpow((long) power, res);
X        res->sign = (BOOL)sign;
X        return;
X    }
X    /*
X     * Handle low powers specially
X     */
X    if (power <= 4) {
X        switch ((int) power) {
X            case 1:
X                ans.len = z1.len;
X                ans.v = alloc(ans.len);
X                copyval(z1, ans);
X                ans.sign = (BOOL)sign;
X                *res = ans;
X                return;
X            case 2:
X                zsquare(z1, res);
X                return;
X            case 3:
X                zsquare(z1, &temp);
X                zmul(z1, temp, res);
X                freeh(temp.v);
X                res->sign = (BOOL)sign;
X                return;
X            case 4:
X                zsquare(z1, &temp);
X                zsquare(temp, res);
X                freeh(temp.v);
X                return;
X        }
X    }
X    /*
X     * Shift out all powers of twos so the multiplies are smaller.
X     * We will shift back the right amount when done.
X     */
X    twos = 0;
X    if (iseven(z1)) {
X        twos = zlowbit(z1);
X        ans.v = alloc(z1.len);
X        ans.len = z1.len;
X        copyval(z1, ans);
X        shiftr(ans, twos);
X        trim(&ans);
X        z1 = ans;
X        twos *= power;
X    }
X    /*
X     * Compute the power by squaring and multiplying.
X     * This uses the left to right method of power raising.
X     */
X    bit = TOPFULL;
X    while ((bit & power) == 0)
X        bit >>= 1L;
X    bit >>= 1L;
X    zsquare(z1, &ans);
X    if (bit & power) {
X        zmul(ans, z1, &temp);
X        freeh(ans.v);
X        ans = temp;
X    }
X    bit >>= 1L;
X    while (bit) {
X        zsquare(ans, &temp);
X        freeh(ans.v);
X        ans = temp;
X        if (bit & power) {
X            zmul(ans, z1, &temp);
X            freeh(ans.v);
X            ans = temp;
X        }
X        bit >>= 1L;
X    }
X    /*
X     * Scale back up by proper power of two
X     */
X    if (twos) {
X        zshift(ans, twos, &temp);
X        freeh(ans.v);
X        ans = temp;
X        freeh(z1.v);
X    }
X    ans.sign = (BOOL)sign;
X    *res = ans;
X}
X
X
X/*
X * Compute ten to the specified power
X * This saves some work since the squares of ten are saved.
X */
Xvoid
Xztenpow(power, res)
X    long power;
X    ZVALUE *res;
X{
X    long i;
X    ZVALUE ans;
X    ZVALUE temp;
X
X    if (power <= 0) {
X        *res = _one_;
X        return;
X    }
X    ans = _one_;
X    _tenpowers_[0] = _ten_;
X    for (i = 0; power; i++) {
X        if (_tenpowers_[i].len == 0)
X            zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
X        if (power & 0x1) {
X            zmul(ans, _tenpowers_[i], &temp);
X            freeh(ans.v);
X            ans = temp;
X        }
X        power /= 2;
X    }
X    *res = ans;
X}
X
X
X/*
X * Calculate modular inverse suppressing unnecessary divisions.
X * This is based on the Euclidian algorithm for large numbers.
X * (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
X * Returns TRUE if there is no solution because the numbers
X * are not relatively prime.
X */
XBOOL
Xzmodinv(u, v, res)
X    ZVALUE u, v;
X    ZVALUE *res;
X{
X    FULL    q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
X    ZVALUE    u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
X
X    if (isneg(u) || isneg(v) || (zrel(u, v) >= 0))
X        zmod(u, v, &v3);
X    else
X        zcopy(u, &v3);
X    zcopy(v, &u3);
X    u2 = _zero_;
X    v2 = _one_;
X
X    /*
X     * Loop here while the size of the numbers remain above
X     * the size of a FULL.  Throughout this loop u3 >= v3.
X     */
X    while ((u3.len > 1) && !iszero(v3)) {
X        uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
X        vh = 0;
X        if ((v3.len + 1) >= u3.len)
X            vh = v3.v[v3.len - 1];
X        if (v3.len == u3.len)
X            vh = (vh << BASEB) + v3.v[v3.len - 2];
X        A = 1;
X        B = 0;
X        C = 0;
X        D = 1;
X
X        /*
X         * Calculate successive quotients of the continued fraction
X         * expansion using only single precision arithmetic until
X         * greater precision is required.
X         */
X        while ((vh + C) && (vh + D)) {
X            q1 = (uh + A) / (vh + C);
X            q2 = (uh + B) / (vh + D);
X            if (q1 != q2)
X                break;
X            T = A - q1 * C;
X            A = C;
X            C = T;
X            T = B - q1 * D;
X            B = D;
X            D = T;
X            T = uh - q1 * vh;
X            uh = vh;
X            vh = T;
X        }
X
X        /*
X         * If B is zero, then we made no progress because
X         * the calculation requires a very large quotient.
X         * So we must do this step of the calculation in
X         * full precision
X         */
X        if (B == 0) {
X            zquo(u3, v3, &qz);
X            zmul(qz, v2, &tmp1);
X            zsub(u2, tmp1, &tmp2);
X            freeh(tmp1.v);
X            freeh(u2.v);
X            u2 = v2;
X            v2 = tmp2;
X            zmul(qz, v3, &tmp1);
X            zsub(u3, tmp1, &tmp2);
X            freeh(tmp1.v);
X            freeh(u3.v);
X            u3 = v3;
X            v3 = tmp2;
X            freeh(qz.v);
X            continue;
X        }
X        /*
X         * Apply the calculated A,B,C,D numbers to the current
X         * values to update them as if the full precision
X         * calculations had been carried out.
X         */
X        zmuli(u2, (long) A, &tmp1);
X        zmuli(v2, (long) B, &tmp2);
X        zadd(tmp1, tmp2, &tmp3);
X        freeh(tmp1.v);
X        freeh(tmp2.v);
X        zmuli(u2, (long) C, &tmp1);
X        zmuli(v2, (long) D, &tmp2);
X        freeh(u2.v);
X        freeh(v2.v);
X        u2 = tmp3;
X        zadd(tmp1, tmp2, &v2);
X        freeh(tmp1.v);
X        freeh(tmp2.v);
X        zmuli(u3, (long) A, &tmp1);
X        zmuli(v3, (long) B, &tmp2);
X        zadd(tmp1, tmp2, &tmp3);
X        freeh(tmp1.v);
X        freeh(tmp2.v);
X        zmuli(u3, (long) C, &tmp1);
X        zmuli(v3, (long) D, &tmp2);
X        freeh(u3.v);
X        freeh(v3.v);
X        u3 = tmp3;
X        zadd(tmp1, tmp2, &v3);
X        freeh(tmp1.v);
X        freeh(tmp2.v);
X    }
X
X    /*
X     * Here when the remaining numbers become single precision in size.
X     * Finish the procedure using single precision calculations.
X     */
X    if (iszero(v3) && !isone(u3)) {
X        freeh(u3.v);
X        freeh(v3.v);
X        freeh(u2.v);
X        freeh(v2.v);
X        return TRUE;
X    }
X    ui3 = (istiny(u3) ? z1tol(u3) : z2tol(u3));
X    vi3 = (istiny(v3) ? z1tol(v3) : z2tol(v3));
X    freeh(u3.v);
X    freeh(v3.v);
X    while (vi3) {
X        q1 = ui3 / vi3;
X        zmuli(v2, (long) q1, &tmp1);
X        zsub(u2, tmp1, &tmp2);
X        freeh(tmp1.v);
X        freeh(u2.v);
X        u2 = v2;
X        v2 = tmp2;
X        q2 = ui3 - q1 * vi3;
X        ui3 = vi3;
X        vi3 = q2;
X    }
X    freeh(v2.v);
X    if (ui3 != 1) {
X        freeh(u2.v);
X        return TRUE;
X    }
X    if (isneg(u2)) {
X        zadd(v, u2, res);
X        freeh(u2.v);
X        return FALSE;
X    }
X    *res = u2;
X    return FALSE;
X}
X
X
X#if 0
X/*
X * Approximate the quotient of two integers by another set of smaller
X * integers.  This uses continued fractions to determine the smaller set.
X */
Xvoid
Xzapprox(z1, z2, res1, res2)
X    ZVALUE z1, z2, *res1, *res2;
X{
X    int sign;
X    ZVALUE u1, v1, u3, v3, q, t1, t2, t3;
X
X    sign = ((z1.sign != 0) ^ (z2.sign != 0));
X    z1.sign = 0;
X    z2.sign = 0;
X    v3 = z2;
X    u3 = z1;
X    u1 = _one_;
X    v1 = _zero_;
X    while (!iszero(v3)) {
X        zdiv(u3, v3, &q, &t1);
X        zmul(v1, q, &t2);
X        zsub(u1, t2, &t3);
X        freeh(q.v);
X        freeh(t2.v);
X        freeh(u1.v);
X        if ((u3.v != z1.v) && (u3.v != z2.v))
X            freeh(u3.v);
X        u1 = v1;
X        u3 = v3;
X        v1 = t3;
X        v3 = t1;
X    }
X    if (!isunit(u3))
X        error("Non-relativly prime numbers for approx");
X    if ((u3.v != z1.v) && (u3.v != z2.v))
X        freeh(u3.v);
X    if ((v3.v != z1.v) && (v3.v != z2.v))
X        freeh(v3.v);
X    freeh(v1.v);
X    zmul(u1, z1, &t1);
X    zsub(t1, _one_, &t2);
X    freeh(t1.v);
X    zquo(t2, z2, &t1);
X    freeh(t2.v);
X    u1.sign = (BOOL)sign;
X    t1.sign = 0;
X    *res1 = t1;
X    *res2 = u1;
X}
X#endif
X
X
X/*
X * Binary gcd algorithm
X * This algorithm taken from Knuth
X */
Xvoid
Xzgcd(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE u, v, t;
X    register long j, k, olen, mask;
X    register HALF h;
X    HALF *oldv1, *oldv2;
X
X    /*
X     * First see if one number is very much larger than the other.
X     * If so, then divide as necessary to get the numbers near each other.
X     */
X    z1.sign = 0;
X    z2.sign = 0;
X    oldv1 = z1.v;
X    oldv2 = z2.v;
X    if (z1.len < z2.len) {
X        t = z1;
X        z1 = z2;
X        z2 = t;
X    }
X    while ((z1.len > (z2.len + 5)) && !iszero(z2)) {
X        zmod(z1, z2, &t);
X        if ((z1.v != oldv1) && (z1.v != oldv2))
X            freeh(z1.v);
X        z1 = z2;
X        z2 = t;
X    }
X    /*
X     * Ok, now do the binary method proper
X     */
X    u.len = z1.len;
X    v.len = z2.len;
X    u.sign = 0;
X    v.sign = 0;
X    if (!ztest(z1)) {
X        v.v = alloc(v.len);
X        copyval(z2, v);
X        *res = v;
X        goto done;
X    }
X    if (!ztest(z2)) {
X        u.v = alloc(u.len);
X        copyval(z1, u);
X        *res = u;
X        goto done;
X    }
X    u.v = alloc(u.len);
X    v.v = alloc(v.len);
X    copyval(z1, u);
X    copyval(z2, v);
X    k = 0;
X    while (u.v[k] == 0 && v.v[k] == 0)
X        ++k;
X    mask = 01;
X    h = u.v[k] | v.v[k];
X    k *= BASEB;
X    while (!(h & mask)) {
X        mask <<= 1;
X        k++;
X    }
X    shiftr(u, k);
X    shiftr(v, k);
X    trim(&u);
X    trim(&v);
X    if (isodd(u)) {
X        t.v = alloc(v.len);
X        t.len = v.len;
X        copyval(v, t);
X        t.sign = !v.sign;
X    } else {
X        t.v = alloc(u.len);
X        t.len = u.len;
X        copyval(u, t);
X        t.sign = u.sign;
X    }
X    while (ztest(t)) {
X        j = 0;
X        while (!t.v[j])
X            ++j;
X        mask = 01;
X        h = t.v[j];
X        j *= BASEB;
X        while (!(h & mask)) {
X            mask <<= 1;
X            j++;
X        }
X        shiftr(t, j);
X        trim(&t);
X        if (ztest(t) > 0) {
X            freeh(u.v);
X            u = t;
X        } else {
X            freeh(v.v);
X            v = t;
X            v.sign = !t.sign;
X        }
X        zsub(u, v, &t);
X    }
X    freeh(t.v);
X    freeh(v.v);
X    if (k) {
X        olen = u.len;
X        u.len += k / BASEB + 1;
X        u.v = (HALF *) realloc(u.v, u.len * sizeof(HALF));
X        while (olen != u.len)
X            u.v[olen++] = 0;
X        shiftl(u, k);
X    }
X    trim(&u);
X    *res = u;
X
Xdone:
X    if ((z1.v != oldv1) && (z1.v != oldv2))
X        freeh(z1.v);
X    if ((z2.v != oldv1) && (z2.v != oldv2))
X        freeh(z2.v);
X}
X
X
X/*
X * Compute the lcm of two integers (least common multiple).
X * This is done using the formula:  gcd(a,b) * lcm(a,b) = a * b.
X */
Xvoid
Xzlcm(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE temp1, temp2;
X
X    zgcd(z1, z2, &temp1);
X    zquo(z1, temp1, &temp2);
X    freeh(temp1.v);
X    zmul(temp2, z2, res);
X    freeh(temp2.v);
X}
X
X
X/*
X * Return whether or not two numbers are relatively prime to each other.
X */
XBOOL
Xzrelprime(z1, z2)
X    ZVALUE z1, z2;            /* numbers to be tested */
X{
X    FULL rem1, rem2;        /* remainders */
X    ZVALUE rem;
X    BOOL result;
X
X    z1.sign = 0;
X    z2.sign = 0;
X    if (iseven(z1) && iseven(z2))    /* false if both even */
X        return FALSE;
X    if (isunit(z1) || isunit(z2))    /* true if either is a unit */
X        return TRUE;
X    if (iszero(z1) || iszero(z2))    /* false if either is zero */
X        return FALSE;
X    if (istwo(z1) || istwo(z2))    /* true if either is two */
X        return TRUE;
X    /*
X     * Try reducing each number by the product of the first few odd primes
X     * to see if any of them are a common factor.
X     */
X    rem1 = zmodi(z1, 3L * 5 * 7 * 11 * 13);
X    rem2 = zmodi(z2, 3L * 5 * 7 * 11 * 13);
X    if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
X        return FALSE;
X    if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
X        return FALSE;
X    if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
X        return FALSE;
X    if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
X        return FALSE;
X    if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
X        return FALSE;
X    /*
X     * Try a new batch of primes now
X     */
X    rem1 = zmodi(z1, 17L * 19 * 23);
X    rem2 = zmodi(z2, 17L * 19 * 23);
X    if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
X        return FALSE;
X    if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
X        return FALSE;
X    if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
X        return FALSE;
X    /*
X     * Yuk, we must actually compute the gcd to know the answer
X     */
X    zgcd(z1, z2, &rem);
X    result = isunit(rem);
X    freeh(rem.v);
X    return result;
X}
X
X
X/*
X * Compute the log of one number base another, to the closest integer.
X * This is the largest integer which when the second number is raised to it,
X * the resulting value is less than or equal to the first number.
X * Example:  zlog(123456, 10) = 5.
X */
Xlong
Xzlog(z1, z2)
X    ZVALUE z1, z2;
X{
X    register ZVALUE *zp;        /* current square */
X    long power;            /* current power */
X    long worth;            /* worth of current square */
X    ZVALUE val;            /* current value of power */
X    ZVALUE temp;            /* temporary */
X    ZVALUE squares[32];        /* table of squares of base */
X
X    /*
X     * Make sure that the numbers are nonzero and the base is greater than one.
X     */
X    if (isneg(z1) || iszero(z1) || isneg(z2) || isleone(z2))
X        error("Bad arguments for log");
X    /*
X     * Reject trivial cases.
X     */
X    if (z1.len < z2.len)
X        return 0;
X    if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))
X        return 0;
X    power = zrel(z1, z2);
X    if (power <= 0)
X        return (power + 1);
X    /*
X     * Handle any power of two special.
X     */
X    if (zisonebit(z2))
X        return (zhighbit(z1) / zlowbit(z2));
X    /*
X     * Handle base 10 special
X     */
X    if ((z2.len == 1) && (*z2.v == 10))
X        return zlog10(z1);
X    /*
X     * Now loop by squaring the base each time, and see whether or
X     * not each successive square is still smaller than the number.
X     */
X    worth = 1;
X    zp = &squares[0];
X    *zp = z2;
X    while (((zp->len * 2) - 1) <= z1.len) {    /* while square not too large */
X        zsquare(*zp, zp + 1);
X        zp++;
X        worth *= 2;
X    }
X    /*
X     * Now back down the squares, and multiply them together to see
X     * exactly how many times the base can be raised by.
X     */
X    val = _one_;
X    power = 0;
X    for (; zp >= squares; zp--, worth /= 2) {
X        if ((val.len + zp->len - 1) <= z1.len) {
X            zmul(val, *zp, &temp);
X            if (zrel(z1, temp) >= 0) {
X                freeh(val.v);
X                val = temp;
X                power += worth;
X            } else
X                freeh(temp.v);
X        }
X        if (zp != squares)
X            freeh(zp->v);
X    }
X    return power;
X}
X
X
X/*
X * Return the integral log base 10 of a number.
X */
Xlong
Xzlog10(z)
X    ZVALUE z;
X{
X    register ZVALUE *zp;        /* current square */
X    long power;            /* current power */
X    long worth;            /* worth of current square */
X    ZVALUE val;            /* current value of power */
X    ZVALUE temp;            /* temporary */
X
X    if (!ispos(z))
X        error("Non-positive number for log10");
X    /*
X     * Loop by squaring the base each time, and see whether or
X     * not each successive square is still smaller than the number.
X     */
X    worth = 1;
X    zp = &_tenpowers_[0];
X    *zp = _ten_;
X    while (((zp->len * 2) - 1) <= z.len) {    /* while square not too large */
X        if (zp[1].len == 0)
X            zsquare(*zp, zp + 1);
X        zp++;
X        worth *= 2;
X    }
X    /*
X     * Now back down the squares, and multiply them together to see
X     * exactly how many times the base can be raised by.
X     */
X    val = _one_;
X    power = 0;
X    for (; zp >= _tenpowers_; zp--, worth /= 2) {
X        if ((val.len + zp->len - 1) <= z.len) {
X            zmul(val, *zp, &temp);
X            if (zrel(z, temp) >= 0) {
X                freeh(val.v);
X                val = temp;
X                power += worth;
X            } else
X                freeh(temp.v);
X        }
X    }
X    return power;
X}
X
X
X/*
X * Return the number of times that one number will divide another.
X * This works similarly to zlog, except that divisions must be exact.
X * For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
X */
Xlong
Xzdivcount(z1, z2)
X    ZVALUE z1, z2;
X{
X    long count;        /* number of factors removed */
X    ZVALUE tmp;        /* ignored return value */
X
X    count = zfacrem(z1, z2, &tmp);
X    freeh(tmp.v);
X    return count;
X}
X
X
X/*
X * Remove all occurances of the specified factor from a number.
X * Also returns the number of factors removed as a function return value.
X * Example:  zfacrem(540, 3, &x) returns 3 and sets x to 20.
X */
Xlong
Xzfacrem(z1, z2, rem)
X    ZVALUE z1, z2, *rem;
X{
X    register ZVALUE *zp;        /* current square */
X    long count;            /* total count of divisions */
X    long worth;            /* worth of current square */
X    ZVALUE temp1, temp2, temp3;    /* temporaries */
X    ZVALUE squares[32];        /* table of squares of factor */
X
X    z1.sign = 0;
X    z2.sign = 0;
X    /*
X     * Make sure factor isn't 0 or 1.
X     */
X    if (isleone(z2))
X        error("Bad argument for facrem");
X    /*
X     * Reject trivial cases.
X     */
X    if ((z1.len < z2.len) || (isodd(z1) && iseven(z2)) ||
X        ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
X        rem->v = alloc(z1.len);
X        rem->len = z1.len;
X        rem->sign = 0;
X        copyval(z1, *rem);
X        return 0;
X    }
X    /*
X     * Handle any power of two special.
X     */
X    if (zisonebit(z2)) {
X        count = zlowbit(z1) / zlowbit(z2);
X        rem->v = alloc(z1.len);
X        rem->len = z1.len;
X        rem->sign = 0;
X        copyval(z1, *rem);
X        shiftr(*rem, count);
X        trim(rem);
X        return count;
X    }
X    /*
X     * See if the factor goes in even once.
X     */
X    zdiv(z1, z2, &temp1, &temp2);
X    if (!iszero(temp2)) {
X        freeh(temp1.v);
X        freeh(temp2.v);
X        rem->v = alloc(z1.len);
X        rem->len = z1.len;
X        rem->sign = 0;
X        copyval(z1, *rem);
X        return 0;
X    }
X    freeh(temp2.v);
X    z1 = temp1;
X    /*
X     * Now loop by squaring the factor each time, and see whether
X     * or not each successive square will still divide the number.
X     */
X    count = 1;
X    worth = 1;
X    zp = &squares[0];
X    *zp = z2;
X    while (((zp->len * 2) - 1) <= z1.len) {    /* while square not too large */
X        zsquare(*zp, &temp1);
X        zdiv(z1, temp1, &temp2, &temp3);
X        if (!iszero(temp3)) {
X            freeh(temp1.v);
X            freeh(temp2.v);
X            freeh(temp3.v);
X            break;
X        }
X        freeh(temp3.v);
X        freeh(z1.v);
X        z1 = temp2;
X        *++zp = temp1;
X        worth *= 2;
X        count += worth;
X    }
X    /*
X     * Now back down the list of squares, and see if the lower powers
X     * will divide any more times.
X     */
X    for (; zp >= squares; zp--, worth /= 2) {
X        if (zp->len <= z1.len) {
X            zdiv(z1, *zp, &temp1, &temp2);
X            if (iszero(temp2)) {
X                temp3 = z1;
X                z1 = temp1;
X                temp1 = temp3;
X                count += worth;
X            }
X            freeh(temp1.v);
X            freeh(temp2.v);
X        }
X        if (zp != squares)
X            freeh(zp->v);
X    }
X    *rem = z1;
X    return count;
X}
X
X
X/*
X * Keep dividing a number by the gcd of it with another number until the
X * result is relatively prime to the second number.
X */
Xvoid
Xzgcdrem(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE tmp1, tmp2;
X
X    /*
X     * Begin by taking the gcd for the first time.
X     * If the number is already relatively prime, then we are done.
X     */
X    zgcd(z1, z2, &tmp1);
X    if (isunit(tmp1) || iszero(tmp1)) {
X        res->len = z1.len;
X        res->v = alloc(z1.len);
X        res->sign = z1.sign;
X        copyval(z1, *res);
X        return;
X    }
X    zquo(z1, tmp1, &tmp2);
X    z1 = tmp2;
X    z2 = tmp1;
X    /*
X     * Now keep alternately taking the gcd and removing factors until
X     * the gcd becomes one.
X     */
X    while (!isunit(z2)) {
X        (void) zfacrem(z1, z2, &tmp1);
X        freeh(z1.v);
X        z1 = tmp1;
X        zgcd(z1, z2, &tmp1);
X        freeh(z2.v);
X        z2 = tmp1;
X    }
X    *res = z1;
X}
X
X
X/*
X * Find the lowest prime factor of a number if one can be found.
X * Search is conducted for the first count primes.
X * Returns 1 if no factor was found.
X */
Xlong
Xzlowfactor(z, count)
X    ZVALUE z;
X    long count;
X{
X    FULL p, d;
X    ZVALUE div, tmp;
X    HALF divval[2];
X
X    if ((--count < 0) || iszero(z))
X        return 1;
X    if (iseven(z))
X        return 2;
X    div.sign = 0;
X    div.v = divval;
X    for (p = 3; (count > 0); p += 2) {
X        for (d = 3; (d * d) <= p; d += 2)
X            if ((p % d) == 0)
X                goto next;
X        divval[0] = (p & BASE1);
X        divval[1] = (p / BASE);
X        div.len = 1 + (p >= BASE);
X        zmod(z, div, &tmp);
X        if (iszero(tmp))
X            return p;
X        freeh(tmp.v);
X        count--;
Xnext:;
X    }
X    return 1;
X}
X
X
X/*
X * Return the number of digits (base 10) in a number, ignoring the sign.
X */
Xlong
Xzdigits(z1)
X    ZVALUE z1;
X{
X    long count, val;
X
X    z1.sign = 0;
X    if (istiny(z1)) {    /* do small numbers ourself */
X        count = 1;
X        val = 10;
X        while (*z1.v >= (HALF)val) {
X            count++;
X            val *= 10;
X        }
X        return count;
X    }
X    return (zlog10(z1) + 1);
X}
X
X
X/*
X * Return the single digit at the specified decimal place of a number,
X * where 0 means the rightmost digit.  Example:  zdigit(1234, 1) = 3.
X */
XFLAG
Xzdigit(z1, n)
X    ZVALUE z1;
X    long n;
X{
X    ZVALUE tmp1, tmp2;
X    FLAG res;
X
X    z1.sign = 0;
X    if (iszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
X        return 0;
X    if (n == 0)
X        return zmodi(z1, 10L);
X    if (n == 1)
X        return zmodi(z1, 100L) / 10;
X    if (n == 2)
X        return zmodi(z1, 1000L) / 100;
X    if (n == 3)
X        return zmodi(z1, 10000L) / 1000;
X    ztenpow(n, &tmp1);
X    zquo(z1, tmp1, &tmp2);
X    res = zmodi(tmp2, 10L);
X    freeh(tmp1.v);
X    freeh(tmp2.v);
X    return res;
X}
X
X
X/*
X * Find the square root of a number.  This is the greatest integer whose
X * square is less than or equal to the number. Returns TRUE if the
X * square root is exact.
X */
XBOOL
Xzsqrt(z1, dest)
X    ZVALUE z1, *dest;
X{
X    ZVALUE try, quo, rem, old, temp;
X    FULL iquo, val;
X    long i,j;
X
X    if (z1.sign)
X        error("Square root of negative number");
X    if (iszero(z1)) {
X        *dest = _zero_;
X        return TRUE;
X    }
X    if ((*z1.v < 4) && istiny(z1)) {
X        *dest = _one_;
X        return (*z1.v == 1);
X    }
X    /*
X     * Pick the square root of the leading one or two digits as a first guess.
X     */
X    val = z1.v[z1.len-1];
X    if ((z1.len & 0x1) == 0)
X        val = (val * BASE) + z1.v[z1.len-2];
X
X    /*
X     * Find the largest power of 2 that when squared
X     * is <= val > 0.  Avoid multiply overflow by doing
X     * a careful check at the BASE boundary.
X     */
X    j = 1<<(BASEB+BASEB-2);
X    if (val > j) {
X        iquo = BASE;
X    } else {
X        i = BASEB-1;
X        while (j > val) {
X            --i;
X            j >>= 2;
X        }
X        iquo = bitmask[i];
X    }
X
X    for (i = 8; i > 0; i--)
X        iquo = (iquo + (val / iquo)) / 2;
X    if (iquo > BASE1)
X        iquo = BASE1;
X    /*
X     * Allocate the numbers to use for the main loop.
X     * The size and high bits of the final result are correctly set here.
X     * Notice that the remainder of the test value is rubbish, but this
X     * is unimportant.
X     */
X    try.sign = 0;
X    try.len = (z1.len + 1) / 2;
X    try.v = alloc(try.len);
X    try.v[try.len-1] = (HALF)iquo;
X    old.sign = 0;
X    old.v = alloc(try.len);
X    old.len = 1;
X    *old.v = 0;
X    /*
X     * Main divide and average loop
X     */
X    for (;;) {
X        zdiv(z1, try, &quo, &rem);
X        i = zrel(try, quo);
X        if ((i == 0) && iszero(rem)) {    /* exact square root */
X            freeh(rem.v);
X            freeh(quo.v);
X            freeh(old.v);
X            *dest = try;
X            return TRUE;
X        }
X        freeh(rem.v);
X        if (i <= 0) {
X            /*
X            * Current try is less than or equal to the square root since
X            * it is less than the quotient.  If the quotient is equal to
X            * the try, we are all done.  Also, if the try is equal to the
X            * old try value, we are done since no improvement occurred.
X            * If not, save the improved value and loop some more.
X            */
X            if ((i == 0) || (zcmp(old, try) == 0)) {
X                freeh(quo.v);
X                freeh(old.v);
X                *dest = try;
X                return FALSE;
X            }
X            old.len = try.len;
X            copyval(try, old);
X        }
X        /* average current try and quotent for the new try */
X        zadd(try, quo, &temp);
X        freeh(quo.v);
X        freeh(try.v);
X        try = temp;
X        shiftr(try, 1L);
X        quicktrim(try);
X    }
X}
X
X
X/*
X * Take an arbitrary root of a number (to the greatest integer).
X * This uses the following iteration to get the Kth root of N:
X *    x = ((K-1) * x + N / x^(K-1)) / K
X */
Xvoid
Xzroot(z1, z2, dest)
X    ZVALUE z1, z2, *dest;
X{
X    ZVALUE try, quo, old, temp, temp2;
X    ZVALUE k1;            /* holds k - 1 */
X    int sign;
X    long i, k, highbit;
X    SIUNION sival;
X
X    sign = z1.sign;
X    if (sign && iseven(z2))
X        error("Even root of negative number");
X    if (iszero(z2) || isneg(z2))
X        error("Non-positive root");
X    if (iszero(z1)) {    /* root of zero */
X        *dest = _zero_;
X        return;
X    }
X    if (isunit(z2)) {    /* first root */
X        zcopy(z1, dest);
X        return;
X    }
X    if (isbig(z2)) {    /* humongous root */
X        *dest = _one_;
X        dest->sign = (HALF)sign;
X        return;
X    }
X    k = (istiny(z2) ? z1tol(z2) : z2tol(z2));
X    highbit = zhighbit(z1);
X    if (highbit < k) {    /* too high a root */
X        *dest = _one_;
X        dest->sign = (HALF)sign;
X        return;
X    }
X    sival.ivalue = k - 1;
X    k1.v = &sival.silow;
X    k1.len = 1 + (sival.sihigh != 0);
X    k1.sign = 0;
X    z1.sign = 0;
X    /*
X     * Allocate the numbers to use for the main loop.
X     * The size and high bits of the final result are correctly set here.
X     * Notice that the remainder of the test value is rubbish, but this
X     * is unimportant.
X     */
X    highbit = (highbit + k - 1) / k;
X    try.len = (highbit / BASEB) + 1;
X    try.v = alloc(try.len);
X    try.v[try.len-1] = ((HALF)1 << (highbit % BASEB));
X    try.sign = 0;
X    old.v = alloc(try.len);
X    old.len = 1;
X    *old.v = 0;
X    old.sign = 0;
X    /*
X     * Main divide and average loop
X     */
X    for (;;) {
X        zpowi(try, k1, &temp);
X        zquo(z1, temp, &quo);
X        freeh(temp.v);
X        i = zrel(try, quo);
X        if (i <= 0) {
X            /*
X             * Current try is less than or equal to the root since it is
X             * less than the quotient. If the quotient is equal to the try,
X             * we are all done.  Also, if the try is equal to the old value,
X             * we are done since no improvement occurred.
X             * If not, save the improved value and loop some more.
X             */
X            if ((i == 0) || (zcmp(old, try) == 0)) {
X                freeh(quo.v);
X                freeh(old.v);
X                try.sign = (HALF)sign;
X                quicktrim(try);
X                *dest = try;
X                return;
X            }
X            old.len = try.len;
X            copyval(try, old);
X        }
X        /* average current try and quotent for the new try */
X        zmul(try, k1, &temp);
X        freeh(try.v);
X        zadd(quo, temp, &temp2);
X        freeh(temp.v);
X        freeh(quo.v);
X        zquo(temp2, z2, &try);
X        freeh(temp2.v);
X    }
X}
X
X
X/*
X * Test to see if a number is an exact square or not.
X */
XBOOL
Xzissquare(z)
X    ZVALUE z;        /* number to be tested */
X{
X    long n, i;
X    ZVALUE tmp;
X
X    if (z.sign)        /* negative */
X        return FALSE;
X    while ((z.len > 1) && (*z.v == 0)) {    /* ignore trailing zero words */
X        z.len--;
X        z.v++;
X    }
X    if (isleone(z))        /* zero or one */
X        return TRUE;
X    n = *z.v & 0xf;        /* check mod 16 values */
X    if ((n != 0) && (n != 1) && (n != 4) && (n != 9))
X        return FALSE;
X    n = *z.v & 0xff;    /* check mod 256 values */
X    i = 0x80;
X    while (((i * i) & 0xff) != n)
X        if (--i <= 0)
X            return FALSE;
X    n = zsqrt(z, &tmp);    /* must do full square root test now */
X    freeh(tmp.v);
X    return n;
X}
X
X/* END CODE */
END_OF_FILE
if test 34345 -ne `wc -c <'zfunc.c'`; then
    echo shar: \"'zfunc.c'\" unpacked with wrong size!
fi
# end of 'zfunc.c'
fi
echo shar: End of archive 18 \(of 21\).
cp /dev/null ark18isdone
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