limit fstBC to 30bp in Python3 ver.

[GalaxyCodeBases.git] / c_cpp / etc / calc / zmath.c

/*
 * zmath - extended precision integral arithmetic primitives
 *
 * Copyright (C) 1999-2007  David I. Bell and Ernest Bowen
 *
 * Primary author:  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 * @(#) $Revision: 30.1 $
 * @(#) $Id: zmath.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
 * @(#) $Source: /usr/local/src/bin/calc/RCS/zmath.c,v $
 *
 * Under source code control:   1990/02/15 01:48:28
 * File existed as early as:    before 1990
 *
 * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
 */


#include <stdio.h>
#include "zmath.h"

HALF _zeroval_[] = { 0 };
HALF _oneval_[] = { 1 };
HALF _twoval_[] = { 2 };
HALF _threeval_[] = { 3 };
HALF _fourval_[] = { 4 };
HALF _fiveval_[] = { 5 };
HALF _sixval_[] = { 6 };
HALF _sevenval_[] = { 7 };
HALF _eightval_[] = { 8 };
HALF _nineval_[] = { 9 };
HALF _tenval_[] = { 10 };
HALF _elevenval_[] = { 11 };
HALF _twelveval_[] = { 12 };
HALF _thirteenval_[] = { 13 };
HALF _fourteenval_[] = { 14 };
HALF _fifteenval_[] = { 15 };
HALF _sixteenval_[] = { 16 };
HALF _seventeenval_[] = { 17 };
HALF _eightteenval_[] = { 18 };
HALF _nineteenval_[] = { 19 };
HALF _twentyval_[] = { 20 };
HALF _sqbaseval_[] = { 0, 1 };
HALF _pow4baseval_[] = { 0, 0, 1 };
HALF _pow8baseval_[] = { 0, 0, 0, 0, 1 };

ZVALUE zconst[] = {
    { _zeroval_, 1, 0}, { _oneval_, 1, 0}, { _twoval_, 1, 0},
    { _threeval_, 1, 0}, { _fourval_, 1, 0}, { _fiveval_, 1, 0},
    { _sixval_, 1, 0}, { _sevenval_, 1, 0}, { _eightval_, 1, 0},
    { _nineval_, 1, 0}, { _tenval_, 1, 0}, { _elevenval_, 1, 0},
    { _twelveval_, 1, 0}, { _thirteenval_, 1, 0}, { _fourteenval_, 1, 0},
    { _fifteenval_, 1, 0}, { _sixteenval_, 1, 0}, { _seventeenval_, 1, 0},
    { _eightteenval_, 1, 0}, { _nineteenval_, 1, 0}, { _twentyval_, 1, 0}
};

ZVALUE _zero_ = { _zeroval_, 1, 0};
ZVALUE _one_ = { _oneval_, 1, 0 };
ZVALUE _two_ = { _twoval_, 1, 0 };
ZVALUE _ten_ = { _tenval_, 1, 0 };
ZVALUE _sqbase_ = { _sqbaseval_, 2, 0 };
ZVALUE _pow4base_ = { _pow4baseval_, 4, 0 };
ZVALUE _pow8base_ = { _pow8baseval_, 4, 0 };
ZVALUE _neg_one_ = { _oneval_, 1, 1 };

/*
 * 2^64 as a ZVALUE
 */
#if BASEB == 32
ZVALUE _b32_ = { _sqbaseval_, 2, 0 };
ZVALUE _b64_ = { _pow4baseval_, 3, 0 };
#elif BASEB == 16
ZVALUE _b32_ = { _pow4baseval_, 3, 0 };
ZVALUE _b64_ = { _pow8baseval_, 5, 0 };
#else
        -=@=- BASEB not 16 or 32 -=@=-
#endif


/*
 * highhalf[i] - masks off the upper i bits of a HALF
 * rhighhalf[i] - masks off the upper BASEB-i bits of a HALF
 * lowhalf[i] - masks off the upper i bits of a HALF
 * rlowhalf[i] - masks off the upper BASEB-i bits of a HALF
 * bitmask[i] - (1 << i) for  0 <= i <= BASEB*2
 */
HALF highhalf[BASEB+1] = {
#if BASEB == 32
        0x00000000,
        0x80000000, 0xC0000000, 0xE0000000, 0xF0000000,
        0xF8000000, 0xFC000000, 0xFE000000, 0xFF000000,
        0xFF800000, 0xFFC00000, 0xFFE00000, 0xFFF00000,
        0xFFF80000, 0xFFFC0000, 0xFFFE0000, 0xFFFF0000,
        0xFFFF8000, 0xFFFFC000, 0xFFFFE000, 0xFFFFF000,
        0xFFFFF800, 0xFFFFFC00, 0xFFFFFE00, 0xFFFFFF00,
        0xFFFFFF80, 0xFFFFFFC0, 0xFFFFFFE0, 0xFFFFFFF0,
        0xFFFFFFF8, 0xFFFFFFFC, 0xFFFFFFFE, 0xFFFFFFFF
#elif BASEB == 16
        0x0000,
        0x8000, 0xC000, 0xE000, 0xF000,
        0xF800, 0xFC00, 0xFE00, 0xFF00,
        0xFF80, 0xFFC0, 0xFFE0, 0xFFF0,
        0xFFF8, 0xFFFC, 0xFFFE, 0xFFFF
#else
        -=@=- BASEB not 16 or 32 -=@=-
#endif
};
HALF rhighhalf[BASEB+1] = {
#if BASEB == 32
        0xFFFFFFFF, 0xFFFFFFFE, 0xFFFFFFFC, 0xFFFFFFF8,
        0xFFFFFFF0, 0xFFFFFFE0, 0xFFFFFFC0, 0xFFFFFF80,
        0xFFFFFF00, 0xFFFFFE00, 0xFFFFFC00, 0xFFFFF800,
        0xFFFFF000, 0xFFFFE000, 0xFFFFC000, 0xFFFF8000,
        0xFFFF0000, 0xFFFE0000, 0xFFFC0000, 0xFFF80000,
        0xFFF00000, 0xFFE00000, 0xFFC00000, 0xFF800000,
        0xFF000000, 0xFE000000, 0xFC000000, 0xF8000000,
        0xF0000000, 0xE0000000, 0xC0000000, 0x80000000,
        0x00000000
#elif BASEB == 16
        0xFFFF, 0xFFFE, 0xFFFC, 0xFFF8,
        0xFFF0, 0xFFE0, 0xFFC0, 0xFF80,
        0xFF00, 0xFE00, 0xFC00, 0xF800,
        0xF000, 0xE000, 0xC000, 0x8000,
        0x0000
#else
        -=@=- BASEB not 16 or 32 -=@=-
#endif
};
HALF lowhalf[BASEB+1] = {
        0x0,
        0x1,  0x3, 0x7, 0xF,
        0x1F,  0x3F, 0x7F, 0xFF,
        0x1FF,  0x3FF, 0x7FF, 0xFFF,
        0x1FFF,  0x3FFF, 0x7FFF, 0xFFFF
#if BASEB == 32
       ,0x1FFFF,  0x3FFFF, 0x7FFFF, 0xFFFFF,
        0x1FFFFF,  0x3FFFFF, 0x7FFFFF, 0xFFFFFF,
        0x1FFFFFF,  0x3FFFFFF, 0x7FFFFFF, 0xFFFFFFF,
        0x1FFFFFFF,  0x3FFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF
#endif
};
HALF rlowhalf[BASEB+1] = {
#if BASEB == 32
        0xFFFFFFFF, 0x7FFFFFFF, 0x3FFFFFFF, 0x1FFFFFFF,
        0xFFFFFFF, 0x7FFFFFF, 0x3FFFFFF, 0x1FFFFFF,
        0xFFFFFF, 0x7FFFFF, 0x3FFFFF, 0x1FFFFF,
        0xFFFFF, 0x7FFFF, 0x3FFFF, 0x1FFFF,
#endif
        0xFFFF, 0x7FFF, 0x3FFF, 0x1FFF,
        0xFFF, 0x7FF, 0x3FF, 0x1FF,
        0xFF, 0x7F, 0x3F, 0x1F,
        0xF, 0x7, 0x3, 0x1,
        0x0
};
HALF bitmask[(2*BASEB)+1] = {
#if BASEB == 32
        0x00000001, 0x00000002, 0x00000004, 0x00000008,
        0x00000010, 0x00000020, 0x00000040, 0x00000080,
        0x00000100, 0x00000200, 0x00000400, 0x00000800,
        0x00001000, 0x00002000, 0x00004000, 0x00008000,
        0x00010000, 0x00020000, 0x00040000, 0x00080000,
        0x00100000, 0x00200000, 0x00400000, 0x00800000,
        0x01000000, 0x02000000, 0x04000000, 0x08000000,
        0x10000000, 0x20000000, 0x40000000, 0x80000000,
        0x00000001, 0x00000002, 0x00000004, 0x00000008,
        0x00000010, 0x00000020, 0x00000040, 0x00000080,
        0x00000100, 0x00000200, 0x00000400, 0x00000800,
        0x00001000, 0x00002000, 0x00004000, 0x00008000,
        0x00010000, 0x00020000, 0x00040000, 0x00080000,
        0x00100000, 0x00200000, 0x00400000, 0x00800000,
        0x01000000, 0x02000000, 0x04000000, 0x08000000,
        0x10000000, 0x20000000, 0x40000000, 0x80000000,
        0x00000001
#elif BASEB == 16
        0x0001, 0x0002, 0x0004, 0x0008,
        0x0010, 0x0020, 0x0040, 0x0080,
        0x0100, 0x0200, 0x0400, 0x0800,
        0x1000, 0x2000, 0x4000, 0x8000,
        0x0001, 0x0002, 0x0004, 0x0008,
        0x0010, 0x0020, 0x0040, 0x0080,
        0x0100, 0x0200, 0x0400, 0x0800,
        0x1000, 0x2000, 0x4000, 0x8000,
        0x0001
#else
        -=@=- BASEB not 16 or 32 -=@=-
#endif
};              /* actual rotation thru 8 cycles */

BOOL _math_abort_;              /* nonzero to abort calculations */


/*
 * popcnt - popcnt[x] number of 1 bits in 0 <= x < 256
 */
char popcnt[256] = {
    0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
};



#ifdef ALLOCTEST
STATIC long nalloc = 0;
STATIC long nfree = 0;
#endif


HALF *
alloc(LEN len)
{
        HALF *hp;

        if (_math_abort_) {
                math_error("Calculation aborted");
                /*NOTREACHED*/
        }
        hp = (HALF *) malloc((len+1) * sizeof(HALF));
        if (hp == 0) {
                math_error("Not enough memory");
                /*NOTREACHED*/
        }
#ifdef ALLOCTEST
        ++nalloc;
#endif
        return hp;
}


#ifdef ALLOCTEST
void
freeh(HALF *h)
{
        if ((h != _zeroval_) && (h != _oneval_)) {
                free(h);
                ++nfree;
        }
}


void
allocStat(void)
{
        fprintf(stderr, "nalloc: %ld nfree: %ld kept: %ld\n",
                nalloc, nfree, nalloc - nfree);
}
#endif


/*
 * Convert a normal integer to a number.
 */
void
itoz(long i, ZVALUE *res)
{
        long diddle, len;

        res->len = 1;
        res->sign = 0;
        diddle = 0;
        if (i == 0) {
                res->v = _zeroval_;
                return;
        }
        if (i < 0) {
                res->sign = 1;
                i = -i;
                if (i < 0) {    /* fix most negative number */
                        diddle = 1;
                        i--;
                }
        }
        if (i == 1) {
                res->v = _oneval_;
                return;
        }
        len = 1 + (((FULL) i) >= BASE);
        res->len = (LEN)len;
        res->v = alloc((LEN)len);
        res->v[0] = (HALF) (i + diddle);
        if (len == 2)
                res->v[1] = (HALF) (i / BASE);
}


/*
 * Convert a number to a normal integer, as far as possible.
 * If the number is out of range, the largest number is returned.
 */
long
ztoi(ZVALUE z)
{
        long i;

        if (zgtmaxlong(z)) {
                i = MAXLONG;
                return (z.sign ? -i : i);
        }
        i = ztolong(z);
        return (z.sign ? -i : i);
}


/*
 * Convert a normal unsigned integer to a number.
 */
void
utoz(FULL i, ZVALUE *res)
{
        long len;

        res->len = 1;
        res->sign = 0;
        if (i == 0) {
                res->v = _zeroval_;
                return;
        }
        if (i == 1) {
                res->v = _oneval_;
                return;
        }
        len = 1 + (((FULL) i) >= BASE);
        res->len = (LEN)len;
        res->v = alloc((LEN)len);
        res->v[0] = (HALF)i;
        if (len == 2)
                res->v[1] = (HALF) (i / BASE);
}


/*
 * Convert a normal signed integer to a number.
 */
void
stoz(SFULL i, ZVALUE *res)
{
        long diddle, len;

        res->len = 1;
        res->sign = 0;
        diddle = 0;
        if (i == 0) {
                res->v = _zeroval_;
                return;
        }
        if (i < 0) {
                res->sign = 1;
                i = -i;
                if (i < 0) {    /* fix most negative number */
                        diddle = 1;
                        i--;
                }
        }
        if (i == 1) {
                res->v = _oneval_;
                return;
        }
        len = 1 + (((FULL) i) >= BASE);
        res->len = (LEN)len;
        res->v = alloc((LEN)len);
        res->v[0] = (HALF) (i + diddle);
        if (len == 2)
                res->v[1] = (HALF) (i / BASE);
}


/*
 * Convert a number to a unsigned normal integer, as far as possible.
 * If the number is out of range, the largest number is returned.
 * The absolute value of z is converted.
 */
FULL
ztou(ZVALUE z)
{
        if (z.len > 2) {
                return MAXUFULL;
        }
        return ztofull(z);
}


/*
 * Convert a number to a signed normal integer, as far as possible.
 *
 * If the number is too large to fit into an integer, than the largest
 * positive or largest negative integer is returned depending on the sign.
 */
SFULL
ztos(ZVALUE z)
{
        FULL val;       /* absolute value of the return value */

        /* negative value processing */
        if (z.sign) {
            if (z.len > 2) {
                return MINSFULL;
            }
            val = ztofull(z);
            if (val > TOPFULL) {
                return MINSFULL;
            }
            return (SFULL)(-val);
        }

        /* positive value processing */
        if (z.len > 2) {
            return (SFULL)MAXFULL;
        }
        val = ztofull(z);
        if (val > MAXFULL) {
            return (SFULL)MAXFULL;
        }
        return (SFULL)val;
}


/*
 * Make a copy of an integer value
 */
void
zcopy(ZVALUE z, ZVALUE *res)
{
        res->sign = z.sign;
        res->len = z.len;
        if (zisabsleone(z)) {   /* zero or plus or minus one are easy */
                res->v = (z.v[0] ? _oneval_ : _zeroval_);
                return;
        }
        res->v = alloc(z.len);
        zcopyval(z, *res);
}


/*
 * Add together two integers.
 */
void
zadd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
        ZVALUE dest;
        HALF *p1, *p2, *pd;
        long len;
        FULL carry;
        SIUNION sival;

        if (z1.sign && !z2.sign) {
                z1.sign = 0;
                zsub(z2, z1, res);
                return;
        }
        if (z2.sign && !z1.sign) {
                z2.sign = 0;
                zsub(z1, z2, res);
                return;
        }
        if (z2.len > z1.len) {
                pd = z1.v; z1.v = z2.v; z2.v = pd;
                len = z1.len; z1.len = z2.len; z2.len = (LEN)len;
        }
        dest.len = z1.len + 1;
        dest.v = alloc(dest.len);
        dest.sign = z1.sign;
        carry = 0;
        pd = dest.v;
        p1 = z1.v;
        p2 = z2.v;
        len = z2.len;
        while (len--) {
                sival.ivalue = ((FULL) *p1++) + ((FULL) *p2++) + carry;
                /* ignore Saber-C warning #112 - get ushort from uint */
                /*        ok to ignore on name zadd`sival */
                *pd++ = sival.silow;
                carry = sival.sihigh;
        }
        len = z1.len - z2.len;
        while (len--) {
                sival.ivalue = ((FULL) *p1++) + carry;
                *pd++ = sival.silow;
                carry = sival.sihigh;
        }
        *pd = (HALF)carry;
        zquicktrim(dest);
        *res = dest;
}


/*
 * Subtract two integers.
 */
void
zsub(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
        register HALF *h1, *h2, *hd;
        long len1, len2;
        FULL carry;
        SIUNION sival;
        ZVALUE dest;

        if (z1.sign != z2.sign) {
                z2.sign = z1.sign;
                zadd(z1, z2, res);
                return;
        }
        len1 = z1.len;
        len2 = z2.len;
        if (len1 == len2) {
                h1 = z1.v + len1;
                h2 = z2.v + len2;
                while ((len1 > 0) && ((FULL)*--h1 == (FULL)*--h2)) {
                        len1--;
                }
                if (len1 == 0) {
                        *res = _zero_;
                        return;
                }
                len2 = len1;
                carry = ((FULL)*h1 < (FULL)*h2);
        } else {
                carry = (len1 < len2);
        }
        dest.sign = z1.sign;
        h1 = z1.v;
        h2 = z2.v;
        if (carry) {
                carry = len1;
                len1 = len2;
                len2 = (long)carry;
                h1 = z2.v;
                h2 = z1.v;
                dest.sign = !dest.sign;
        }
        hd = alloc((LEN)len1);
        dest.v = hd;
        dest.len = (LEN)len1;
        len1 -= len2;
        carry = 0;
        while (--len2 >= 0) {
                /* ignore Saber-C warning #112 - get ushort from uint */
                /*        ok to ignore on name zsub`sival */
                sival.ivalue = (BASE1 - ((FULL) *h1++)) + *h2++ + carry;
                *hd++ = (HALF)(BASE1 - sival.silow);
                carry = sival.sihigh;
        }
        while (--len1 >= 0) {
                sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
                *hd++ = (HALF)(BASE1 - sival.silow);
                carry = sival.sihigh;
        }
        if (hd[-1] == 0)
                ztrim(&dest);
        *res = dest;
}


/*
 * Multiply an integer by a small number.
 */
void
zmuli(ZVALUE z, long n, ZVALUE *res)
{
        register HALF *h1, *sd;
        FULL low;
        FULL high;
        FULL carry;
        long len;
        SIUNION sival;
        ZVALUE dest;

        if ((n == 0) || ziszero(z)) {
                *res = _zero_;
                return;
        }
        if (n < 0) {
                n = -n;
                z.sign = !z.sign;
        }
        if (n == 1) {
                zcopy(z, res);
                return;
        }
#if LONG_BITS > BASEB
        low = ((FULL) n) & BASE1;
        high = ((FULL) n) >> BASEB;
#else
        low = (FULL)n;
        high = 0;
#endif
        dest.len = z.len + 2;
        dest.v = alloc(dest.len);
        dest.sign = z.sign;
        /*
         * Multiply by the low digit.
         */
        h1 = z.v;
        sd = dest.v;
        len = z.len;
        carry = 0;
        while (len--) {
                /* ignore Saber-C warning #112 - get ushort from uint */
                /*        ok to ignore on name zmuli`sival */
                sival.ivalue = ((FULL) *h1++) * low + carry;
                *sd++ = sival.silow;
                carry = sival.sihigh;
        }
        *sd = (HALF)carry;
        /*
         * If there was only one digit, then we are all done except
         * for trimming the number if there was no last carry.
         */
        if (high == 0) {
                dest.len--;
                if (carry == 0)
                        dest.len--;
                *res = dest;
                return;
        }
        /*
         * Need to multiply by the high digit and add it into the
         * previous value.  Clear the final word of rubbish first.
         */
        *(++sd) = 0;
        h1 = z.v;
        sd = dest.v + 1;
        len = z.len;
        carry = 0;
        while (len--) {
                sival.ivalue = ((FULL) *h1++) * high + ((FULL) *sd) + carry;
                *sd++ = sival.silow;
                carry = sival.sihigh;
        }
        *sd = (HALF)carry;
        zquicktrim(dest);
        *res = dest;
}


/*
 * Divide two numbers by their greatest common divisor.
 * This is useful for reducing the numerator and denominator of
 * a fraction to its lowest terms.
 */
void
zreduce(ZVALUE z1, ZVALUE z2, ZVALUE *z1res, ZVALUE *z2res)
{
        ZVALUE tmp;

        if (zisabsleone(z1) || zisabsleone(z2))
                tmp = _one_;
        else
                zgcd(z1, z2, &tmp);
        if (zisunit(tmp)) {
                zcopy(z1, z1res);
                zcopy(z2, z2res);
        } else {
                zequo(z1, tmp, z1res);
                zequo(z2, tmp, z2res);
        }
        zfree(tmp);
}



/*
 * Compute the quotient and remainder for division of an integer by an
 * integer, rounding criteria determined by rnd.  Returns the sign of
 * the remainder.
 */
long
zdiv(ZVALUE z1, ZVALUE z2, ZVALUE *quo, ZVALUE *rem, long rnd)
{
        register HALF *a, *b;
        HALF s, u;
        HALF *A, *B, *a1, *b0;
        FULL f, g, h, x;
        BOOL adjust, onebit;
        LEN m, n, len, i, p, j1, j2, k;
        long t, val;

        if (ziszero(z2)) {
                math_error("Division by zero in zdiv");
                /*NOTREACHED*/
        }
        m = z1.len;
        n = z2.len;
        B = z2.v;
        s = 0;
        if (m < n) {
                A = alloc(n + 1);
                memcpy(A, z1.v, m * sizeof(HALF));
                for (i = m; i <= n; i++)
                        A[i] = 0;
                a1 = A + n;
                len = 1;
                goto done;
        }
        A = alloc(m + 2);
        memcpy(A, z1.v, m * sizeof(HALF));
        A[m] = 0;
        A[m + 1] = 0;
        len = m - n + 1;        /* quotient length will be len or len +/- 1 */
        a1 = A + n;             /* start of digits for quotient */
        b0 = B;
        p = n;
        while (!*b0) {          /* b0: working start for divisor */
                ++b0;
                --p;
        }
        if (p == 1) {
                u = *b0;
                if (u == 1) {
                        for (; m >= n; m--)
                                A[m] = A[m - 1];
                        A[m] = 0;
                        goto done;
                }
                f = 0;
                a = A + m;
                i = len;
                while (i--) {
                        f = f << BASEB | *--a;
                        a[1] = (HALF)(f / u);
                        f = f % u;
                }
                *a = (HALF)f;
                m = n;
                goto done;
        }
        f = B[n - 1];
        k = 1;
        while (f >>= 1)
                k++;            /* k: number of bits in top divisor digit */
        j1 = BASEB - k;
        j2 = BASEB + j1;
        h = j1 ? ((FULL) B[n - 1] << j1 | B[n - 2] >> k) : B[n-1];
        onebit = (BOOL)((B[n - 2] >> (k - 1)) & 1);
        m++;
        while (m > n) {
                m--;
                f = (FULL) A[m] << j2 | (FULL) A[m - 1] << j1;
                if (j1) f |= A[m - 2] >> k;
                if (s) f = ~f;
                x = f / h;
                if (x) {
                        if (onebit && x > TOPHALF + f % h)
                                x--;
                        a = A + m - p;
                        b = b0;
                        u = 0;
                        i = p;
                        if (s) {
                                while (i--) {
                                        f = (FULL) *a + u + x * *b++;
                                        *a++ = (HALF) f;
                                        u = (HALF) (f >> BASEB);
                                }
                                s = *a + u;
                                A[m] = (HALF) (~x + !s);
                        } else {
                                while (i--) {
                                        f = (FULL) *a - u - x * *b++;
                                        *a++ = (HALF) f;
                                        u = -(HALF)(f >> BASEB);
                                }
                                s = *a - u;
                                A[m] = (HALF)(x + s);
                        }
                }
                else
                        A[m] = s;
        }
done:   while (m > 0 && A[m - 1] == 0)
                m--;
        if (m == 0 && s == 0) {
                *rem = _zero_;
                val = 0;
                if (a1[len - 1] == 0)
                        len--;
                if (len == 0) {
                        *quo = _zero_;
                } else {
                        quo->len = len;
                        quo->v = alloc(len);
                        memcpy(quo->v, a1, len * sizeof(HALF));
                        quo->sign = z1.sign ^ z2.sign;
                }
                freeh(A);
                return val;
        }
        if (rnd & 8)
                adjust = (((*a1 ^ rnd) & 1) ? TRUE : FALSE);
        else
                adjust = (((rnd & 1) ^ z1.sign ^ z2.sign) ? TRUE : FALSE);
        if (rnd & 2)
                adjust ^= z1.sign ^ z2.sign;
        if (rnd & 4)
                adjust ^= z2.sign;
        if (rnd & 16) {                 /* round-off case */
                a = A + n;
                b = B + n;
                i = n + 1;
                f = g = 0;
                t = -1;
                if (s) {
                        while (--i > 0) {
                                g = (FULL) *--a + (*--b >> 1 | f);
                                f = *b & 1 ? TOPHALF : 0;
                                if (g != BASE1)
                                        break;
                        }
                        if (g == BASE && f == 0) {
                                while ((--i > 0) && ((*--a | *--b) == 0));
                                t = (i > 0);
                        } else if (g >= BASE) {
                                t = 1;
                        }
                } else {
                        while (--i > 0) {
                                g = (FULL) *--a - (*--b >> 1 | f);
                                f = *b & 1 ? TOPHALF : 0;
                                if (g != 0)
                                        break;
                        }
                        if (g > 0 && g < BASE)
                                t = 1;
                        else if (g == 0 && f == 0)
                                t = 0;
                }
                if (t)
                        adjust = (t > 0);
        }
        if (adjust) {
                i = len;
                a = a1;
                while (i > 0 && *a == BASE1) {
                        i--;
                        *a++ = 0;
                }
                (*a)++;
                if (i == 0)
                        len++;
        }
        if (s && adjust) {
                i = 0;
                while (A[i] == 0)
                        i++;
                A[i] = -A[i];
                while (++i < n)
                        A[i] = ~A[i];
                m = n;
                while (A[m - 1] == 0)
                        m--;
        }
        if (!s && adjust) {
                a = A;
                b = B;
                i = n;
                u = 0;
                while (i--) {
                        f = (FULL) *b++ - *a - (FULL) u;
                        *a++ = (HALF) f;
                        u = -(HALF)(f >> BASEB);
                }
                m = n;
                while (A[m - 1] == 0)
                        m--;
        }
        if (s && !adjust) {
                a = A;
                b = B;
                i = n;
                f = 0;
                while (i--) {
                        f = (FULL) *b++ + *a + (f >> BASEB);
                        *a++ = (HALF) f;
                }
                m = n;
                while (A[m-1] == 0)
                        m--;
        }
        rem->len = m;
        rem->v = alloc(m);
        memcpy(rem->v, A, m * sizeof(HALF));
        rem->sign = z1.sign ^ adjust;
        val = rem->sign ? -1 : 1;
        if (a1[len - 1] == 0)
                len--;
        if (len == 0) {
                *quo = _zero_;
        } else {
                quo->len = len;
                quo->v = alloc(len);
                memcpy(quo->v, a1, len * sizeof(HALF));
                quo->sign = z1.sign ^ z2.sign;
        }
        freeh(A);
        return val;
}


/*
 * Compute and store at a specified location the integer quotient
 * of two integers, the type of rounding being determined by rnd.
 * Returns the sign of z1/z2 - calculated quotient.
 */
long
zquo(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
{
        ZVALUE tmp;
        long val;

        val = zdiv(z1, z2, res, &tmp, rnd);
        if (z2.sign)
                val = -val;
        zfree(tmp);
        return val;
}


/*
 * Compute and store at a specified location the remainder for
 * division of an integer by an integer, the type of rounding
 * used being determined by rnd.  Returns the sign of the remainder.
 */
long
zmod(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
{
        ZVALUE tmp;
        long val;

        val = zdiv(z1, z2, &tmp, res, rnd);
        zfree(tmp);
        return val;
}


/*
 * Computes the quotient z1/z2 on the assumption that this is exact.
 * There is no thorough check on the exactness of the division
 * so this should not be called unless z1/z2 is an integer.
 */
void
zequo(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
        LEN i, j, k, len, m, n, o, p;
        HALF *a, *a0, *A, *b, *B, u, v, w, x;
        FULL f, g;

        if (ziszero(z1)) {
                *res = _zero_;
                return;
        }
        if (ziszero(z2)) {
                math_error("Division by zero");
                /*NOTREACHED*/
        }
        if (zisone(z2)) {
                zcopy(z1, res);
                return;
        }
        if (zhighbit(z1) < zhighbit(z2)) {
                math_error("Bad call to zequo");
                /*NOTREACHED*/
        }
        B = z2.v;
        o = 0;
        while (!*B) {
                ++B;
                ++o;
        }
        m = z1.len - o;
        n = z2.len - o;
        len = m - n + 1;                /* Maximum length of quotient */
        v = *B;
        A = alloc(len+1);
        memcpy(A, z1.v + o, len * sizeof(HALF));
        A[len] = 0;
        if (n == 1) {
                if (v > 1) {
                        a = A + len;
                        i = len;
                        f = 0;
                        while (i--) {
                                f = f << BASEB | *--a;
                                *a = (HALF)(f / v);
                                f %= v;
                        }
                }
        } else {
                k = 0;
                while (!(v & 1)) {
                        k++;
                        v >>= 1;
                }
                j = BASEB - k;
                if (n > 1 && k > 0)
                        v |= B[1] << j;
                u = v - 1;
                w = x = 1;
                while (u) {     /* To find w = inverse of v modulo BASE */
                        do {
                                v <<= 1;
                                x <<= 1;
                        }
                        while (!(u & x));
                        u += v;
                        w |= x;
                }
                a0 = A;
                p = len;
                while (p > 1) {
                        if (!*a0) {
                                while (!*++a0 && p > 1)
                                        p--;
                                --a0;
                        }
                        if (p == 1)
                                break;
                        x = k ? w * (*a0 >> k | a0[1] << j) : w * *a0;
                        g = x;
                        if (x) {
                                a = a0;
                                b = B;
                                u = 0;
                                i = n;
                                if (i > p)
                                        i = p;
                                while (i--) {
                                        f = (FULL) *a - g * *b++ - (FULL) u;
                                        *a++ = (HALF)f;
                                        u = -(HALF)(f >> BASEB);
                                }
                                if (u && p > n) {
                                        i = p - n;
                                        while (u && i--) {
                                                f = (FULL) *a - u;
                                                *a++ = (HALF) f;
                                                u = -(HALF)(f >> BASEB);
                                        }
                                }
                        }
                        *a0++ = x;
                        p--;
                }
                if (k == 0) {
                        *a0 = w * *a0;
                } else {
                        u = (HALF)(w * *a0) >> k;
                        x = (HALF)(((FULL) z1.v[z1.len - 1] << BASEB
                                | z1.v[z1.len - 2])
                                /((FULL) B[n-1] << BASEB | B[n-2]));
                        if ((x ^ u) & 1) x--;
                        *a0 = x;
                }
        }
        if (!A[len - 1]) len--;
        res->v = A;
        res->len = len;
        res->sign = z1.sign != z2.sign;
}



/*
 * Return the quotient and remainder of an integer divided by a small
 * number.  A nonzero remainder is only meaningful when both numbers
 * are positive.
 */
long
zdivi(ZVALUE z, long n, ZVALUE *res)
{
        HALF *h1, *sd;
        FULL val;
        HALF divval[2];
        ZVALUE div;
        ZVALUE dest;
        LEN len;

        if (n == 0) {
                math_error("Division by zero");
                /*NOTREACHED*/
        }
        if (ziszero(z)) {
                *res = _zero_;
                return 0;
        }
        if (n < 0) {
                n = -n;
                z.sign = !z.sign;
        }
        if (n == 1) {
                zcopy(z, res);
                return 0;
        }
        /*
         * If the division is by a large number, then call the normal
         * divide routine.
         */
        if (n & ~BASE1) {
                div.sign = 0;
                div.v = divval;
                divval[0] = (HALF) n;
#if LONG_BITS > BASEB
                divval[1] = (HALF)(((FULL) n) >> BASEB);
                div.len = 2;
#else
                div.len = 1;
#endif
                zdiv(z, div, res, &dest, 0);
                n = ztolong(dest);
                zfree(dest);
                return n;
        }
        /*
         * Division is by a small number, so we can be quick about it.
         */
        len = z.len;
        dest.sign = z.sign;
        dest.len = len;
        dest.v = alloc(len);
        h1 = z.v + len;
        sd = dest.v + len;
        val = 0;
        while (len--) {
                val = ((val << BASEB) + ((FULL) *--h1));
                *--sd = (HALF)(val / n);
                val %= n;
        }
        zquicktrim(dest);
        *res = dest;
        return (long)val;
}



/*
 * Calculate the mod of an integer by a small number.
 * This is only defined for positive moduli.
 */
long
zmodi(ZVALUE z, long n)
{
        register HALF *h1;
        FULL val;
        HALF divval[2];
        ZVALUE div;
        ZVALUE temp;
        long len;

        if (n == 0) {
                math_error("Division by zero");
                /*NOTREACHED*/
        }
        if (n < 0) {
                math_error("Non-positive modulus");
                /*NOTREACHED*/
        }
        if (ziszero(z) || (n == 1))
                return 0;
        if (zisone(z))
                return 1;
        /*
         * If the modulus is by a large number, then call the normal
         * modulo routine.
         */
        if (n & ~BASE1) {
                div.sign = 0;
                div.v = divval;
                divval[0] = (HALF) n;
#if LONG_BITS > BASEB
                divval[1] = (HALF)(((FULL) n) >> BASEB);
                div.len = 2;
#else
                div.len = 1;
#endif
                zmod(z, div, &temp, 0);
                n = ztolong(temp);
                zfree(temp);
                return n;
        }
        /*
         * The modulus is by a small number, so we can do this quickly.
         */
        len = z.len;
        h1 = z.v + len;
        val = 0;
        while (len-- > 0)
                val = ((val << BASEB) + ((FULL) *--h1)) % n;
        if (val && z.sign)
                val = n - val;
        return (long)val;
}


/*
 * Return whether or not one number exactly divides another one.
 * Returns TRUE if division occurs with no remainder.
 * z1 is the number to be divided by z2.
 */

BOOL
zdivides(ZVALUE z1, ZVALUE z2)
{
        LEN i, j, k, m, n;
        HALF u, v, w, x;
        HALF *a, *a0, *A, *b, *B, *c, *d;
        FULL f;
        BOOL ans;

        if (zisunit(z2) || ziszero(z1)) return TRUE;
        if (ziszero(z2)) return FALSE;

        m = z1.len;
        n = z2.len;
        if (m < n) return FALSE;

        c = z1.v;
        d = z2.v;

        while (!*d) {
                if (*c++) return FALSE;
                d++;
                m--;
                n--;
        }

        j = 0;
        u = *c;
        v = *d;
        while (!(v & 1)) {      /* Counting and checking zero bits */
                if (u & 1) return FALSE;
                u >>= 1;
                v >>= 1;
                j++;
        }

        if (n == 1 && v == 1) return TRUE;
        if (!*c) {                      /* Removing any further zeros */
                while(!*++c)
                        m--;
                c--;
        }

        if (m < n) return FALSE;

        if (j) {
                B = alloc((LEN)n);      /* Array for shifted z2 */
                d += n;
                b = B + n;
                i = n;
                f = 0;
                while(i--) {
                        f = f << BASEB | *--d;
                        *--b = (HALF)(f >> j);
                }
                if (!B[n - 1]) n--;
        }
        else B = d;
        u = *B;
        v = x = 1;
        w = 0;
        while (x) {                     /* Finding minv(*B, BASE) */
                if (v & x) {
                        v -= x * u;
                        w |= x;
                }
                x <<= 1;
        }

        A = alloc((LEN)(m + 2));        /* Main working array */
        memcpy(A, c, m * sizeof(HALF));
        A[m + 1] = 1;
        A[m] = 0;

        k = m - n + 1;                  /* Length of presumed quotient */

        a0 = A;

        while (k--) {
                if (*a0) {
                        x = w * *a0;
                        a = a0;
                        b = B;
                        i = n;
                        u = 0;
                        while (i--) {
                                f = (FULL) *a - (FULL) x * *b++ - u;
                                *a++ = (HALF)f;
                                u = -(HALF)(f >> BASEB);
                        }
                        f = (FULL) *a - u;
                        *a++ = (HALF)f;
                        u = -(HALF)(f >> BASEB);
                        if (u) {
                                while (*a == 0) *a++ = BASE1;
                                (*a)--;
                        }
                }
                a0++;
        }
        ans = FALSE;
        if (A[m + 1]) {
                a = A + m;
                while (--n && !*--a);
                if (!n) ans = TRUE;
        }
        freeh(A);
        if (j) freeh(B);
        return ans;
}


/*
 * Compute the bitwise OR of two integers
 */
void
zor(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
        register HALF *sp, *dp;
        long len;
        ZVALUE bz, lz, dest;

        if (z1.len >= z2.len) {
                bz = z1;
                lz = z2;
        } else {
                bz = z2;
                lz = z1;
        }
        dest.len = bz.len;
        dest.v = alloc(dest.len);
        dest.sign = 0;
        zcopyval(bz, dest);
        len = lz.len;
        sp = lz.v;
        dp = dest.v;
        while (len--)
                *dp++ |= *sp++;
        *res = dest;
}


/*
 * Compute the bitwise AND of two integers
 */
void
zand(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
        HALF *h1, *h2, *hd;
        LEN len;
        ZVALUE dest;

        len = ((z1.len <= z2.len) ? z1.len : z2.len);
        h1 = &z1.v[len-1];
        h2 = &z2.v[len-1];
        while ((len > 1) && ((*h1 & *h2) == 0)) {
                h1--;
                h2--;
                len--;
        }
        dest.len = len;
        dest.v = alloc(len);
        dest.sign = 0;
        h1 = z1.v;
        h2 = z2.v;
        hd = dest.v;
        while (len--)
                *hd++ = (*h1++ & *h2++);
        *res = dest;
}


/*
 * Compute the bitwise XOR of two integers.
 */
void
zxor(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
        HALF *dp, *h1, *h2;
        LEN len, j, k;
        ZVALUE dest;

        h1 = z1.v;
        h2 = z2.v;
        len = z1.len;
        j = z2.len;
        if (z1.len < z2.len) {
                len = z2.len;
                j = z1.len;
                h1 = z2.v;
                h2 = z1.v;
        } else if (z1.len == z2.len) {
                while (len > 1 && z1.v[len-1] == z2.v[len-1])
                        len--;
                j = len;
        }
        k = len - j;
        dest.len = len;
        dest.v = alloc(len);
        dest.sign = 0;
        dp = dest.v;
        while (j-- > 0)
                *dp++ = *h1++ ^ *h2++;
        while (k-- > 0)
                *dp++ = *h1++;
        *res = dest;
}


/*
 * Compute the bitwise ANDNOT of two integers.
 */
void
zandnot(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
        HALF *dp, *h1, *h2;
        LEN len, j, k;
        ZVALUE dest;

        len = z1.len;
        if (z2.len >= len) {
                while (len > 1 && (z1.v[len-1] & ~z2.v[len-1]) == 0)
                        len--;
                j = len;
                k = 0;
        } else {
                j = z2.len;
                k = len - z2.len;
        }
        dest.len = len;
        dest.v = alloc(len);
        dest.sign = 0;
        dp = dest.v;
        h1 = z1.v;
        h2 = z2.v;
        while (j-- > 0)
                *dp++ = *h1++ & ~*h2++;
        while (k-- > 0)
                *dp++ = *h1++;
        *res = dest;
}


/*
 * Shift a number left (or right) by the specified number of bits.
 * Positive shift means to the left.  When shifting right, rightmost
 * bits are lost.  The sign of the number is preserved.
 */
void
zshift(ZVALUE z, long n, ZVALUE *res)
{
        ZVALUE ans;
        LEN hc;         /* number of halfwords shift is by */

        if (ziszero(z)) {
                *res = _zero_;
                return;
        }
        if (n == 0) {
                zcopy(z, res);
                return;
        }
        /*
         * If shift value is negative, then shift right.
         * Check for large shifts, and handle word-sized shifts quickly.
         */
        if (n < 0) {
                n = -n;
                if ((n < 0) || (n >= (z.len * BASEB))) {
                        *res = _zero_;
                        return;
                }
                hc = (LEN)(n / BASEB);
                n %= BASEB;
                z.v += hc;
                z.len -= hc;
                ans.len = z.len;
                ans.v = alloc(ans.len);
                ans.sign = z.sign;
                zcopyval(z, ans);
                if (n > 0) {
                        zshiftr(ans, n);
                        ztrim(&ans);
                }
                if (ziszero(ans)) {
                        zfree(ans);
                        ans = _zero_;
                }
                *res = ans;
                return;
        }
        /*
         * Shift value is positive, so shift leftwards.
         * Check specially for a shift of the value 1, since this is common.
         * Also handle word-sized shifts quickly.
         */
        if (zisunit(z)) {
                zbitvalue(n, res);
                res->sign = z.sign;
                return;
        }
        hc = (LEN)(n / BASEB);
        n %= BASEB;
        ans.len = z.len + hc + 1;
        ans.v = alloc(ans.len);
        ans.sign = z.sign;
        if (hc > 0)
                memset((char *) ans.v, 0, hc * sizeof(HALF));
        memcpy((char *) (ans.v + hc),
            (char *) z.v, z.len * sizeof(HALF));
        ans.v[ans.len - 1] = 0;
        if (n > 0) {
                ans.v += hc;
                ans.len -= hc;
                zshiftl(ans, n);
                ans.v -= hc;
                ans.len += hc;
        }
        ztrim(&ans);
        *res = ans;
}


/*
 * Return the position of the lowest bit which is set in the binary
 * representation of a number (counting from zero).  This is the highest
 * power of two which evenly divides the number.
 */
long
zlowbit(ZVALUE z)
{
        register HALF *zp;
        long n;
        HALF dataval;
        HALF *bitval;

        n = 0;
        for (zp = z.v; *zp == 0; zp++)
                if (++n >= z.len)
                        return 0;
        dataval = *zp;
        bitval = bitmask;
        /* ignore Saber-C warning #530 about empty while statement */
        /*        ok to ignore in proc zlowbit */
        while ((*(bitval++) & dataval) == 0) {
        }
        return (n*BASEB)+(bitval-bitmask-1);
}


/*
 * Return the position of the highest bit which is set in the binary
 * representation of a number (counting from zero).  This is the highest power
 * of two which is less than or equal to the number (which is assumed nonzero).
 */
LEN
zhighbit(ZVALUE z)
{
        HALF dataval;
        HALF *bitval;

        dataval = z.v[z.len-1];
        if (dataval == 0) {
                return 0;
        }
        bitval = bitmask+BASEB;
        if (dataval) {
                /* ignore Saber-C warning #530 about empty while statement */
                /*        ok to ignore in proc zhighbit */
                while ((*(--bitval) & dataval) == 0) {
                }
        }
        return (z.len*BASEB)+(LEN)(bitval-bitmask-BASEB);
}


/*
 * Return whether or not the specifed bit number is set in a number.
 * Rightmost bit of a number is bit 0.
 */
BOOL
zisset(ZVALUE z, long n)
{
        if ((n < 0) || ((n / BASEB) >= z.len))
                return FALSE;
        return ((z.v[n / BASEB] & (((HALF) 1) << (n % BASEB))) != 0);
}


/*
 * Check whether or not a number has exactly one bit set, and
 * thus is an exact power of two.  Returns TRUE if so.
 */
BOOL
zisonebit(ZVALUE z)
{
        register HALF *hp;
        register LEN len;

        if (ziszero(z) || zisneg(z))
                return FALSE;
        hp = z.v;
        len = z.len;
        while (len > 4) {
                len -= 4;
                if (*hp++ || *hp++ || *hp++ || *hp++)
                        return FALSE;
        }
        while (--len > 0) {
                if (*hp++)
                        return FALSE;
        }
        return ((*hp & -*hp) == *hp);           /* NEEDS 2'S COMPLEMENT */
}


/*
 * Check whether or not a number has all of its bits set below some
 * bit position, and thus is one less than an exact power of two.
 * Returns TRUE if so.
 */
BOOL
zisallbits(ZVALUE z)
{
        register HALF *hp;
        register LEN len;
        HALF digit;

        if (ziszero(z) || zisneg(z))
                return FALSE;
        hp = z.v;
        len = z.len;
        while (len > 4) {
                len -= 4;
                if ((*hp++ != BASE1) || (*hp++ != BASE1) ||
                        (*hp++ != BASE1) || (*hp++ != BASE1))
                                return FALSE;
        }
        while (--len > 0) {
                if (*hp++ != BASE1)
                        return FALSE;
        }
        digit = (HALF)(*hp + 1);
        return ((digit & -digit) == digit);     /* NEEDS 2'S COMPLEMENT */
}


/*
 * Return the number whose binary representation contains only one bit which
 * is in the specified position (counting from zero).  This is equivilant
 * to raising two to the given power.
 */
void
zbitvalue(long n, ZVALUE *res)
{
        ZVALUE z;

        if (n < 0) n = 0;
        z.sign = 0;
        z.len = (LEN)((n / BASEB) + 1);
        z.v = alloc(z.len);
        zclearval(z);
        z.v[z.len-1] = (((HALF) 1) << (n % BASEB));
        *res = z;
}


/*
 * Compare a number against zero.
 * Returns the sgn function of the number (-1, 0, or 1).
 */
FLAG
ztest(ZVALUE z)
{
        register int sign;
        register HALF *h;
        register long len;

        sign = 1;
        if (z.sign)
                sign = -sign;
        h = z.v;
        len = z.len;
        while (len--)
                if (*h++)
                        return sign;
        return 0;
}


/*
 * Return the sign of z1 - z2, i.e. 1 if the first integer is greater,
 * 0 if they are equal, -1 otherwise.
 */
FLAG
zrel(ZVALUE z1, ZVALUE z2)
{
        HALF *h1, *h2;
        LEN len;
        int sign;

        sign = 1;
        if (z1.sign < z2.sign)
                return 1;
        if (z2.sign < z1.sign)
                return -1;
        if (z2.sign)
                sign = -1;
        if (z1.len != z2.len)
                return (z1.len > z2.len) ? sign : -sign;
        len = z1.len;
        h1 = z1.v + len;
        h2 = z2.v + len;

        while (len > 0) {
                if (*--h1 != *--h2)
                        break;
                len--;
        }
        if (len > 0)
                return (*h1 > *h2) ? sign : -sign;
        return 0;
}


/*
 * Return the sign of abs(z1) - abs(z2), i.e. 1 if the first integer
 * has greater absolute value, 0 is they have equal absolute value,
 * -1 otherwise.
 */
FLAG
zabsrel(ZVALUE z1, ZVALUE z2)
{
        HALF *h1, *h2;
        LEN len;

        if (z1.len != z2.len)
                return (z1.len > z2.len) ? 1 : -1;

        len = z1.len;
        h1 = z1.v + len;
        h2 = z2.v + len;
        while (len > 0) {
                if (*--h1 != *--h2)
                        break;
                len--;
        }
        if (len > 0)
                return (*h1 > *h2) ? 1 : -1;
        return 0;
}


/*
 * Compare two numbers to see if they are equal or not.
 * Returns TRUE if they differ.
 */
BOOL
zcmp(ZVALUE z1, ZVALUE z2)
{
        register HALF *h1, *h2;
        register long len;

        if ((z1.sign != z2.sign) || (z1.len != z2.len) || (*z1.v != *z2.v))
                return TRUE;
        len = z1.len;
        h1 = z1.v;
        h2 = z2.v;
        while (--len > 0) {
                if (*++h1 != *++h2)
                        return TRUE;
        }
        return FALSE;
}


/*
 * Utility to calculate the gcd of two FULL integers.
 */
FULL
uugcd(FULL f1, FULL f2)
{
        FULL temp;

        while (f1) {
                temp = f2 % f1;
                f2 = f1;
                f1 = temp;
        }
        return (FULL) f2;
}


/*
 * Utility to calculate the gcd of two small integers.
 */
long
iigcd(long i1, long i2)
{
        FULL f1, f2, temp;

        f1 = (FULL) ((i1 >= 0) ? i1 : -i1);
        f2 = (FULL) ((i2 >= 0) ? i2 : -i2);
        while (f1) {
                temp = f2 % f1;
                f2 = f1;
                f1 = temp;
        }
        return (long) f2;
}


void
ztrim(ZVALUE *z)
{
        HALF *h;
        LEN len;

        h = z->v + z->len - 1;
        len = z->len;
        while (*h == 0 && len > 1) {
                --h;
                --len;
        }
        z->len = len;
}


/*
 * Utility routine to shift right.
 *
 * NOTE: The ZVALUE length is not adjusted instead, the value is
 *       zero padded from the left.  One may need to call ztrim()
 *       or use zshift() instead.
 */
void
zshiftr(ZVALUE z, long n)
{
        register HALF *h, *lim;
        FULL mask, maskt;
        long len;

        if (n >= BASEB) {
                len = n / BASEB;
                h = z.v;
                lim = z.v + z.len - len;
                while (h < lim) {
                        h[0] = h[len];
                        ++h;
                }
                n -= BASEB * len;
                lim = z.v + z.len;
                while (h < lim)
                        *h++ = 0;
        }
        if (n) {
                len = z.len;
                h = z.v + len;
                mask = 0;
                while (len--) {
                        maskt = (((FULL) *--h) << (BASEB - n)) & BASE1;
                        *h = ((*h >> n) | (HALF)mask);
                        mask = maskt;
                }
        }
}


/*
 * Utility routine to shift left.
 *
 * NOTE: The ZVALUE length is not adjusted.  The bits in the upper
 *       HALF are simply tossed.  You may want to use zshift() instead.
 */
void
zshiftl(ZVALUE z, long n)
{
        register HALF *h;
        FULL mask, i;
        long len;

        if (n >= BASEB) {
                len = n / BASEB;
                h = z.v + z.len - 1;
                while (!*h)
                        --h;
                while (h >= z.v) {
                        h[len] = h[0];
                        --h;
                }
                n -= BASEB * len;
                while (len)
                        h[len--] = 0;
        }
        if (n > 0) {
                len = z.len;
                h = z.v;
                mask = 0;
                while (len--) {
                        i = (((FULL) *h) << n) | mask;
                        if (i > BASE1) {
                                mask = i >> BASEB;
                                i &= BASE1;
                        } else {
                                mask = 0;
                        }
                        *h = (HALF) i;
                        ++h;
                }
        }
}


/*
 * popcnt - count the number of 0 or 1 bits in an integer
 *
 * We ignore all 0 bits above the highest bit.
 */
long
zpopcnt(ZVALUE z, int bitval)
{
        long cnt = 0;   /* number of times found */
        HALF h;         /* HALF to count */
        int i;

        /*
         * count 1's
         */
        if (bitval) {

                /*
                 * count each HALF
                 */
                for (i=0; i < z.len; ++i) {
                        /* count each octet */
                        for (h = z.v[i]; h; h >>= 8) {
                                cnt += (long)popcnt[h & 0xff];
                        }
                }

        /*
         * count 0's
         */
        } else {

                /*
                 * count each HALF up until the last
                 */
                for (i=0; i < z.len-1; ++i) {

                        /* count each octet */
                        cnt += BASEB;
                        for (h = z.v[i]; h; h >>= 8) {
                                cnt -= (long)popcnt[h & 0xff];
                        }
                }

                /*
                 * count the last octet up until the highest 1 bit
                 */
                for (h = z.v[z.len-1]; h; h>>=1) {
                        /* count each 0 bit */
                        if ((h & 0x1) == 0) {
                                ++cnt;
                        }
                }
        }

        /*
         * return count
         */
        return cnt;
}