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

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

#ifdef MULTIPLE_THREADS
static void *__dtoa_locks[] = { NULL, NULL };
#endif

 Bigint *
Balloc
#ifdef KR_headers
        (k) int k;
#else
        (int k)
#endif
{
        int x;
        Bigint *rv;
#ifndef Omit_Private_Memory
        unsigned int len;
#endif

        ACQUIRE_DTOA_LOCK(0);
        /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
        /* but this case seems very unlikely. */
        if (k <= Kmax && (rv = freelist[k]) !=0) {
                freelist[k] = rv->next;
                }
        else {
                x = 1 << k;
#ifdef Omit_Private_Memory
                rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
                if (rv == NULL)
                        return (NULL);
#else
                len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
                        /sizeof(double);
                if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
                        rv = (Bigint*)pmem_next;
                        pmem_next += len;
                        }
                else {
                        rv = (Bigint*)MALLOC(len*sizeof(double));
                        if (rv == NULL)
                                return (NULL);
                }
#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) {
                if (v->k > Kmax)
#ifdef FREE
                        FREE(v);
#else
                        free(v);
#endif
                else {
                        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
{
        int k;
        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 = 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);
                        if (b1 == NULL)
                                return (NULL);
                        Bcopy(b1, b);
                        Bfree(b);
                        b = b1;
                        }
                b->x[wds++] = carry;
                b->wds = wds;
                }
        return b;
        }

 int
hi0bits_D2A
#ifdef KR_headers
        (x) ULong x;
#else
        (ULong x)
#endif
{
        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);
        if (b == NULL)
                return (NULL);
        b->x[0] = 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);
        if (c == NULL)
                return (NULL);
        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 = 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;
        static int p05[3] = { 5, 25, 125 };

        if ( (i = k & 3) !=0) {
                b = multadd(b, p05[i-1], 0);
                if (b == NULL)
                        return (NULL);
                }

        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);
                        if (p5 == NULL)
                                return (NULL);
                        p5->next = 0;
                        }
                FREE_DTOA_LOCK(1);
#else
                p5 = p5s = i2b(625);
                if (p5 == NULL)
                        return (NULL);
                p5->next = 0;
#endif
                }
        for(;;) {
                if (k & 1) {
                        b1 = mult(b, p5);
                        if (b1 == NULL)
                                return (NULL);
                        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);
                                if (p51 == NULL)
                                        return (NULL);
                                p51->next = 0;
                                }
                        FREE_DTOA_LOCK(1);
#else
                        p51 = p5->next = mult(p5,p5);
                        if (p51 == NULL)
                                return (NULL);
                        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);
        if (b1 == NULL)
                return (NULL);
        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);
                if (c == NULL)
                        return (NULL);
                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);
        if (c == NULL)
                return (NULL);
        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;
        U 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
        (dd, e, bits) double dd; int *e, *bits;
#else
        (double dd, int *e, int *bits)
#endif
{
        Bigint *b;
        U d;
#ifndef Sudden_Underflow
        int i;
#endif
        int de, k;
        ULong *x, y, z;
#ifdef VAX
        ULong d0, d1;
#else
#define d0 word0(&d)
#define d1 word1(&d)
#endif
        d.d = dd;
#ifdef VAX
        d0 = word0(&d) >> 16 | word0(&d) << 16;
        d1 = word1(&d) >> 16 | word1(&d) << 16;
#endif

#ifdef Pack_32
        b = Balloc(1);
#else
        b = Balloc(2);
#endif
        if (b == NULL)
                return (NULL);
        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
                };

#ifdef NO_STRING_H

 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;
        }

 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
{
        char *a = (char*)a1, *ae = a + len;
        char *b = (char*)b1, *a0 = a;
        while(a < ae)
                *a++ = *b++;
        return a0;
        }

#endif /* NO_STRING_H */