misc.c

Go to the documentation of this file.

/****************************************************************

The author of this software is David M. Gay.

Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved

Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.

LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.

****************************************************************/

/* Please send bug reports to David M. Gay (dmg at acm dot org,
 * with " at " changed at "@" and " dot " changed to ".").  */

#include "gdtoaimp.h"

static Bigint* freelist[Kmax + 1];
#ifndef Omit_Private_Memory
#ifndef PRIVATE_MEM
#define PRIVATE_MEM 2304
#endif
#define PRIVATE_mem ((PRIVATE_MEM + sizeof(double) - 1) / sizeof(double))
static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
#endif

Bigint* Balloc
#ifdef KR_headers
    (k) int k;
#else
    (int k)
#endif
{
    Bigint* rv;

    ACQUIRE_DTOA_LOCK(0);
    if((rv = freelist[k]) != 0)
    {
        freelist[k] = rv->next;
    }
    else
    {
        int x = 1 << k;
#ifdef Omit_Private_Memory
        rv = (Bigint*)MALLOC(sizeof(Bigint) + (x - 1) * sizeof(ULong));
#else
#ifndef Omit_Private_Memory
        unsigned int len;
#endif
        len = (sizeof(Bigint) + ((unsigned)x - 1) * sizeof(ULong) + sizeof(double) - 1) /
              sizeof(double);
        if((unsigned)(pmem_next - private_mem + len) <= PRIVATE_mem)
        {
            rv = (Bigint*)pmem_next;
            pmem_next += len;
        }
        else
        {
            {
                rv = (Bigint*)MALLOC(len * sizeof(double));
            }
        }
#endif
        rv->k = k;
        rv->maxwds = x;
    }
    FREE_DTOA_LOCK(0);
    rv->sign = rv->wds = 0;
    return rv;
}

void Bfree
#ifdef KR_headers
    (v) Bigint* v;
#else
    (Bigint* v)
#endif
{
    if(v)
    {
        ACQUIRE_DTOA_LOCK(0);
        v->next = freelist[v->k];
        freelist[v->k] = v;
        FREE_DTOA_LOCK(0);
    }
}

int lo0bits
#ifdef KR_headers
    (y) ULong* y;
#else
    (ULong* y)
#endif
{
    register int k;
    register ULong x = *y;

    if(x & 7)
    {
        if(x & 1)
        {
            {
                return 0;
            }
        }
        if(x & 2)
        {
            *y = x >> 1;
            return 1;
        }
        *y = x >> 2;
        return 2;
    }
    k = 0;
    if(!(x & 0xffff))
    {
        k = 16;
        x >>= 16;
    }
    if(!(x & 0xff))
    {
        k += 8;
        x >>= 8;
    }
    if(!(x & 0xf))
    {
        k += 4;
        x >>= 4;
    }
    if(!(x & 0x3))
    {
        k += 2;
        x >>= 2;
    }
    if(!(x & 1))
    {
        k++;
        x >>= 1;
        if(!x)
        {
            {
                return 32;
            }
        }
    }
    *y = x;
    return k;
}

Bigint* multadd
#ifdef KR_headers
    (b, m, a) Bigint* b;
int m, a;
#else
    (Bigint* b, int m, int a) /* multiply by m and add a */
#endif
{
    int i, wds;
#ifdef ULLong
    ULong* x;
    ULLong carry, y;
#else
    ULong carry, *x, y;
#ifdef Pack_32
    ULong xi, z;
#endif
#endif
    Bigint* b1;

    wds = b->wds;
    x = b->x;
    i = 0;
    carry = (unsigned long long)a;
    do
    {
#ifdef ULLong
        y = *x * (ULLong)m + carry;
        carry = y >> 32;
        *x++ = y & 0xffffffffUL;
#else
#ifdef Pack_32
        xi = *x;
        y = (xi & 0xffff) * m + carry;
        z = (xi >> 16) * m + (y >> 16);
        carry = z >> 16;
        *x++ = (z << 16) + (y & 0xffff);
#else
        y = *x * m + carry;
        carry = y >> 16;
        *x++ = y & 0xffff;
#endif
#endif
    } while(++i < wds);
    if(carry)
    {
        if(wds >= b->maxwds)
        {
            b1 = Balloc(b->k + 1);
            Bcopy(b1, b);
            Bfree(b);
            b = b1;
        }
        b->x[wds++] = (ULong)carry;
        b->wds = wds;
    }
    return b;
}

int hi0bits_D2A
#ifdef KR_headers
    (x) register ULong x;
#else
    (register ULong x)
#endif
{
    register int k = 0;

    if(!(x & 0xffff0000))
    {
        k = 16;
        x <<= 16;
    }
    if(!(x & 0xff000000))
    {
        k += 8;
        x <<= 8;
    }
    if(!(x & 0xf0000000))
    {
        k += 4;
        x <<= 4;
    }
    if(!(x & 0xc0000000))
    {
        k += 2;
        x <<= 2;
    }
    if(!(x & 0x80000000))
    {
        k++;
        if(!(x & 0x40000000))
        {
            {
                return 32;
            }
        }
    }
    return k;
}

Bigint* i2b
#ifdef KR_headers
    (i) int i;
#else
    (int i)
#endif
{
    Bigint* b;

    b = Balloc(1);
    b->x[0] = (ULong)i;
    b->wds = 1;
    return b;
}

Bigint *mult
#ifdef KR_headers
    (a, b) Bigint *a,
    *b;
#else
    (Bigint* a, Bigint* b)
#endif
{
    Bigint* c;
    int k, wa, wb, wc;
    ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
    ULong y;
#ifdef ULLong
    ULLong carry, z;
#else
    ULong carry, z;
#ifdef Pack_32
    ULong z2;
#endif
#endif

    if(a->wds < b->wds)
    {
        c = a;
        a = b;
        b = c;
    }
    k = a->k;
    wa = a->wds;
    wb = b->wds;
    wc = wa + wb;
    if(wc > a->maxwds)
    {
        {
            k++;
        }
    }
    c = Balloc(k);
    for(x = c->x, xa = x + wc; x < xa; x++)
    {
        {
            *x = 0;
        }
    }
    xa = a->x;
    xae = xa + wa;
    xb = b->x;
    xbe = xb + wb;
    xc0 = c->x;
#ifdef ULLong
    for(; xb < xbe; xc0++)
    {
        if((y = *xb++) != 0)
        {
            x = xa;
            xc = xc0;
            carry = 0;
            do
            {
                z = *x++ * (ULLong)y + *xc + carry;
                carry = z >> 32;
                *xc++ = z & 0xffffffffUL;
            } while(x < xae);
            *xc = (ULong)carry;
        }
    }
#else
#ifdef Pack_32
    for(; xb < xbe; xb++, xc0++)
    {
        if((y = *xb & 0xffff) != 0)
        {
            x = xa;
            xc = xc0;
            carry = 0;
            do
            {
                z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
                carry = z >> 16;
                z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
                carry = z2 >> 16;
                Storeinc(xc, z2, z);
            } while(x < xae);
            *xc = carry;
        }
        if((y = *xb >> 16) != 0)
        {
            x = xa;
            xc = xc0;
            carry = 0;
            z2 = *xc;
            do
            {
                z = (*x & 0xffff) * y + (*xc >> 16) + carry;
                carry = z >> 16;
                Storeinc(xc, z, z2);
                z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
                carry = z2 >> 16;
            } while(x < xae);
            *xc = z2;
        }
    }
#else
    for(; xb < xbe; xc0++)
    {
        if((y = *xb++) != 0)
        {
            x = xa;
            xc = xc0;
            carry = 0;
            do
            {
                z = *x++ * y + *xc + carry;
                carry = z >> 16;
                *xc++ = z & 0xffff;
            } while(x < xae);
            *xc = carry;
        }
    }
#endif
#endif
    for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc)
    {
        {
            ;
        }
    }
    c->wds = wc;
    return c;
}

static Bigint* p5s;

Bigint* pow5mult
#ifdef KR_headers
    (b, k) Bigint* b;
int k;
#else
    (Bigint* b, int k)
#endif
{
    Bigint *b1, *p5, *p51;
    int i;

    if((i = k & 3) != 0)
    {
        {
            static int p05[3] = {5, 25, 125};
            b = multadd(b, p05[i - 1], 0);
        }
    }

    if(!(k >>= 2))
    {
        {
            return b;
        }
    }
    if((p5 = p5s) == 0)
    {
        /* first time */
#ifdef MULTIPLE_THREADS
        ACQUIRE_DTOA_LOCK(1);
        if(!(p5 = p5s))
        {
            p5 = p5s = i2b(625);
            p5->next = 0;
        }
        FREE_DTOA_LOCK(1);
#else
        p5 = p5s = i2b(625);
        p5->next = 0;
#endif
    }
    for(;;)
    {
        if(k & 1)
        {
            b1 = mult(b, p5);
            Bfree(b);
            b = b1;
        }
        if(!(k >>= 1))
        {
            {
                break;
            }
        }
        if((p51 = p5->next) == 0)
        {
#ifdef MULTIPLE_THREADS
            ACQUIRE_DTOA_LOCK(1);
            if(!(p51 = p5->next))
            {
                p51 = p5->next = mult(p5, p5);
                p51->next = 0;
            }
            FREE_DTOA_LOCK(1);
#else
            p51 = p5->next = mult(p5, p5);
            p51->next = 0;
#endif
        }
        p5 = p51;
    }
    return b;
}

Bigint* lshift
#ifdef KR_headers
    (b, k) Bigint* b;
int k;
#else
    (Bigint* b, int k)
#endif
{
    int i, k1, n, n1;
    Bigint* b1;
    ULong *x, *x1, *xe, z;

    n = k >> kshift;
    k1 = b->k;
    n1 = n + b->wds + 1;
    for(i = b->maxwds; n1 > i; i <<= 1)
    {
        {
            k1++;
        }
    }
    b1 = Balloc(k1);
    x1 = b1->x;
    for(i = 0; i < n; i++)
    {
        {
            *x1++ = 0;
        }
    }
    x = b->x;
    xe = x + b->wds;
    if(k &= kmask)
    {
#ifdef Pack_32
        k1 = 32 - k;
        z = 0;
        do
        {
            *x1++ = *x << k | z;
            z = *x++ >> k1;
        } while(x < xe);
        if((*x1 = z) != 0)
        {
            {
                ++n1;
            }
        }
#else
        k1 = 16 - k;
        z = 0;
        do
        {
            *x1++ = *x << k & 0xffff | z;
            z = *x++ >> k1;
        } while(x < xe);
        if(*x1 = z)
            ++n1;
#endif
    }
    else
    {
        {
            do
            {
                {
                    *x1++ = *x++;
                }
            } while(x < xe);
        }
    }
    b1->wds = n1 - 1;
    Bfree(b);
    return b1;
}

int cmp
#ifdef KR_headers
    (a, b) Bigint *a,
    *b;
#else
    (Bigint* a, Bigint* b)
#endif
{
    ULong *xa, *xa0, *xb, *xb0;
    int i, j;

    i = a->wds;
    j = b->wds;
#ifdef DEBUG
    if(i > 1 && !a->x[i - 1])
        Bug("cmp called with a->x[a->wds-1] == 0");
    if(j > 1 && !b->x[j - 1])
        Bug("cmp called with b->x[b->wds-1] == 0");
#endif
    if(i -= j)
    {
        {
            return i;
        }
    }
    xa0 = a->x;
    xa = xa0 + j;
    xb0 = b->x;
    xb = xb0 + j;
    for(;;)
    {
        if(*--xa != *--xb)
        {
            {
                return *xa < *xb ? -1 : 1;
            }
        }
        if(xa <= xa0)
        {
            {
                break;
            }
        }
    }
    return 0;
}

Bigint *diff
#ifdef KR_headers
    (a, b) Bigint *a,
    *b;
#else
    (Bigint* a, Bigint* b)
#endif
{
    Bigint* c;
    int i, wa, wb;
    ULong *xa, *xae, *xb, *xbe, *xc;
#ifdef ULLong
    ULLong borrow, y;
#else
    ULong borrow, y;
#ifdef Pack_32
    ULong z;
#endif
#endif

    i = cmp(a, b);
    if(!i)
    {
        c = Balloc(0);
        c->wds = 1;
        c->x[0] = 0;
        return c;
    }
    if(i < 0)
    {
        c = a;
        a = b;
        b = c;
        i = 1;
    }
    else
    {
        {
            i = 0;
        }
    }
    c = Balloc(a->k);
    c->sign = i;
    wa = a->wds;
    xa = a->x;
    xae = xa + wa;
    wb = b->wds;
    xb = b->x;
    xbe = xb + wb;
    xc = c->x;
    borrow = 0;
#ifdef ULLong
    do
    {
        y = (ULLong)*xa++ - *xb++ - borrow;
        borrow = y >> 32 & 1UL;
        *xc++ = y & 0xffffffffUL;
    } while(xb < xbe);
    while(xa < xae)
    {
        y = *xa++ - borrow;
        borrow = y >> 32 & 1UL;
        *xc++ = y & 0xffffffffUL;
    }
#else
#ifdef Pack_32
    do
    {
        y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
        borrow = (y & 0x10000) >> 16;
        z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
        borrow = (z & 0x10000) >> 16;
        Storeinc(xc, z, y);
    } while(xb < xbe);
    while(xa < xae)
    {
        y = (*xa & 0xffff) - borrow;
        borrow = (y & 0x10000) >> 16;
        z = (*xa++ >> 16) - borrow;
        borrow = (z & 0x10000) >> 16;
        Storeinc(xc, z, y);
    }
#else
    do
    {
        y = *xa++ - *xb++ - borrow;
        borrow = (y & 0x10000) >> 16;
        *xc++ = y & 0xffff;
    } while(xb < xbe);
    while(xa < xae)
    {
        y = *xa++ - borrow;
        borrow = (y & 0x10000) >> 16;
        *xc++ = y & 0xffff;
    }
#endif
#endif
    while(!*--xc)
    {
        {
            wa--;
        }
    }
    c->wds = wa;
    return c;
}

double b2d
#ifdef KR_headers
    (a, e) Bigint* a;
int* e;
#else
    (Bigint* a, int* e)
#endif
{
    ULong *xa, *xa0, w, y, z;
    int k;
    double d;
#ifdef VAX
    ULong d0, d1;
#else
#define d0 word0(d)
#define d1 word1(d)
#endif

    xa0 = a->x;
    xa = xa0 + a->wds;
    y = *--xa;
#ifdef DEBUG
    if(!y)
        Bug("zero y in b2d");
#endif
    k = hi0bits(y);
    *e = 32 - k;
#ifdef Pack_32
    if(k < Ebits)
    {
        d0 = Exp_1 | y >> (Ebits - k);
        w = xa > xa0 ? *--xa : 0;
        d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k);
        goto ret_d;
    }
    z = xa > xa0 ? *--xa : 0;
    if(k -= Ebits)
    {
        d0 = Exp_1 | y << k | z >> (32 - k);
        y = xa > xa0 ? *--xa : 0;
        d1 = z << k | y >> (32 - k);
    }
    else
    {
        d0 = Exp_1 | y;
        d1 = z;
    }
#else
    if(k < Ebits + 16)
    {
        z = xa > xa0 ? *--xa : 0;
        d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
        w = xa > xa0 ? *--xa : 0;
        y = xa > xa0 ? *--xa : 0;
        d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
        goto ret_d;
    }
    z = xa > xa0 ? *--xa : 0;
    w = xa > xa0 ? *--xa : 0;
    k -= Ebits + 16;
    d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
    y = xa > xa0 ? *--xa : 0;
    d1 = w << k + 16 | y << k;
#endif
ret_d:
#ifdef VAX
    word0(d) = d0 >> 16 | d0 << 16;
    word1(d) = d1 >> 16 | d1 << 16;
#endif
    return dval(d);
}
#undef d0
#undef d1

Bigint* d2b
#ifdef KR_headers
    (d, e, bits) double d;
int *e, *bits;
#else
    (double d, int* e, int* bits)
#endif
{
    Bigint* b;
#ifndef Sudden_Underflow
    int i;
#endif
    int de, k;
    ULong *x, y, z;
#ifdef VAX
    ULong d0, d1;
    d0 = word0(d) >> 16 | word0(d) << 16;
    d1 = word1(d) >> 16 | word1(d) << 16;
#else
#define d0 word0(d)
#define d1 word1(d)
#endif

#ifdef Pack_32
    b = Balloc(1);
#else
    b = Balloc(2);
#endif
    x = b->x;

    z = d0 & Frac_mask;
    d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
#ifdef Sudden_Underflow
    de = (int)(d0 >> Exp_shift);
#ifndef IBM
    z |= Exp_msk11;
#endif
#else
    if((de = (int)(d0 >> Exp_shift)) != 0)
    {
        {
            z |= Exp_msk1;
        }
    }
#endif
#ifdef Pack_32
    if((y = d1) != 0)
    {
        if((k = lo0bits(&y)) != 0)
        {
            x[0] = y | z << (32 - k);
            z >>= k;
        }
        else
        {
            {
                x[0] = y;
            }
        }
#ifndef Sudden_Underflow
        i =
#endif
            b->wds = (x[1] = z) != 0 ? 2 : 1;
    }
    else
    {
        k = lo0bits(&z);
        x[0] = z;
#ifndef Sudden_Underflow
        i =
#endif
            b->wds = 1;
        k += 32;
    }
#else
    if((y = d1) != 0)
    {
        if((k = lo0bits(&y)) != 0)
            if(k >= 16)
            {
                x[0] = y | z << 32 - k & 0xffff;
                x[1] = z >> k - 16 & 0xffff;
                x[2] = z >> k;
                i = 2;
            }
            else
            {
                x[0] = y & 0xffff;
                x[1] = y >> 16 | z << 16 - k & 0xffff;
                x[2] = z >> k & 0xffff;
                x[3] = z >> k + 16;
                i = 3;
            }
        else
        {
            x[0] = y & 0xffff;
            x[1] = y >> 16;
            x[2] = z & 0xffff;
            x[3] = z >> 16;
            i = 3;
        }
    }
    else
    {
#ifdef DEBUG
        if(!z)
            Bug("Zero passed to d2b");
#endif
        k = lo0bits(&z);
        if(k >= 16)
        {
            x[0] = z;
            i = 0;
        }
        else
        {
            x[0] = z & 0xffff;
            x[1] = z >> 16;
            i = 1;
        }
        k += 32;
    }
    while(!x[i])
        --i;
    b->wds = i + 1;
#endif
#ifndef Sudden_Underflow
    if(de)
    {
#endif
#ifdef IBM
        *e = (de - Bias - (P - 1) << 2) + k;
        *bits = 4 * P + 8 - k - hi0bits(word0(d) & Frac_mask);
#else
    *e = de - Bias - (P - 1) + k;
    *bits = P - k;
#endif
#ifndef Sudden_Underflow
    }
    else
    {
        *e = de - Bias - (P - 1) + 1 + k;
#ifdef Pack_32
        *bits = 32 * i - hi0bits(x[i - 1]);
#else
        *bits = (i + 2) * 16 - hi0bits(x[i]);
#endif
    }
#endif
    return b;
}
#undef d0
#undef d1

CONST double
#ifdef IEEE_Arith
    bigtens[] = {1e16, 1e32, 1e64, 1e128, 1e256};
CONST double tinytens[] = {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
#else
#ifdef IBM
    bigtens[] = {1e16, 1e32, 1e64};
CONST double tinytens[] = {1e-16, 1e-32, 1e-64};
#else
    bigtens[] = {1e16, 1e32};
CONST double tinytens[] = {1e-16, 1e-32};
#endif
#endif

CONST double tens[] = {1e0,
                       1e1,
                       1e2,
                       1e3,
                       1e4,
                       1e5,
                       1e6,
                       1e7,
                       1e8,
                       1e9,
                       1e10,
                       1e11,
                       1e12,
                       1e13,
                       1e14,
                       1e15,
                       1e16,
                       1e17,
                       1e18,
                       1e19,
                       1e20,
                       1e21,
                       1e22
#ifdef VAX
                       ,
                       1e23,
                       1e24
#endif
};

char*
#ifdef KR_headers
    strcp_D2A(a, b) char* a;
char* b;
#else
strcp_D2A(char *a, CONST char *b)
#endif
{
    while((*a = *b++))
    {
        {
            a++;
        }
    }
    return a;
}

#ifdef NO_STRING_H

Char*
#ifdef KR_headers
    memcpy_D2A(a, b, len) Char* a;
Char* b;
size_t len;
#else
memcpy_D2A(void *a1, void *b1, size_t len)
#endif
{
    register char *a = (char*)a1, *ae = a + len;
    register char *b = (char*)b1, *a0 = a;
    while(a < ae)
        *a++ = *b++;
    return a0;
}

#endif /* NO_STRING_H */