StdLib/LibC/gdtoa/misc.c

/* $NetBSD: misc.c,v 1.3.12.1 2008/04/08 21:10:55 jdc Exp $ */

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

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  <LibConfig.h>

#include "gdtoaimp.h"

#if defined(_MSC_VER)
  // Disable warnings about assignment within conditional expressions.
  #pragma warning ( disable : 4706 )
#endif

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
{
  int x;
  Bigint *rv;
#ifndef Omit_Private_Memory
  unsigned int len;
#endif

  ACQUIRE_DTOA_LOCK(0);
  if ( (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));
#else
    len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
      /sizeof(double);
    if (pmem_next - private_mem + len <= PRIVATE_mem) {
      rv = (Bigint*)(void *)pmem_next;
      pmem_next += len;
      }
    else
      rv = (Bigint*)MALLOC(len*sizeof(double));
#endif
    if (rv == NULL)
      return NULL;
    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
{
  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;
    /* LINTED conversion */
    *x++ = (uint32_t)(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) {
        Bfree(b);
        return NULL;
        }
      Bcopy(b1, b);
      Bfree(b);
      b = b1;
      }
    /* LINTED conversion */
      b->x[wds++] = (uint32_t)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;
        /* LINTED conversion */
        *xc++ = (uint32_t)(z & 0xffffffffUL);
        }
        while(x < xae);
      /* LINTED conversion */
        *xc = (uint32_t)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 CONST 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 = (unsigned int)k >> 2) == 0)
    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;
      b = b1;
      }
    if ((k = (unsigned int)k >> 1) == 0)
      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 = (unsigned int)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;
    /* LINTED conversion */
    *xc++ = (uint32_t)(y & 0xffffffffUL);
    }
    while(xb < xbe);
  while(xa < xae) {
    y = *xa++ - borrow;
    borrow = y >> 32 & 1UL;
    /* LINTED conversion */
    *xc++ = (uint32_t)(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 = (UINT32)(Exp_1 | y >> (Ebits - k));
    w = xa > xa0 ? *--xa : 0;
    d1 = (UINT32)(y << ((32-Ebits) + k) | w >> (Ebits - k));
    goto ret_d;
    }
  z = xa > xa0 ? *--xa : 0;
  if (k -= Ebits) {
    d0 = (UINT32)(Exp_1 | y << k | z >> (32 - k));
    y = xa > xa0 ? *--xa : 0;
    d1 = (UINT32)(z << k | y >> (32 - k));
    }
  else {
    d0 = (UINT32)(Exp_1 | y);
    d1 = (UINT32)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
  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 {
#ifdef DEBUG
    if (!z)
      Bug("Zero passed to d2b");
#endif
    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
{
  char *a = (char*)a1, *ae = a + len;
  char *b = (char*)b1, *a0 = a;
  while(a < ae)
    *a++ = *b++;
  return a0;
  }

#endif /* NO_STRING_H */