Logo Search packages:      
Sourcecode: mingw-w64 version File versions  Download package

misc.c

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

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 ".").  */


#if defined(__MINGW32__) || defined(__MINGW64__)
/* we have to include windows.h before gdtoa
   headers, otherwise defines cause conflicts. */
#ifndef WIN32_LEAN_AND_MEAN
#define WIN32_LEAN_AND_MEAN
#endif
#include <windows.h>

#define NLOCKS 2

#ifdef USE_WIN32_SL
/* Use spin locks. */
static long dtoa_sl[NLOCKS];

#define ACQUIRE_DTOA_LOCK(n) \
  while (InterlockedCompareExchange (&dtoa_sl[n], 1, 0) != 0) \
     Sleep (0);
#define FREE_DTOA_LOCK(n) InterlockedExchange (&dtoa_sl[n], 0);

#else /* USE_WIN32_SL */

#include <stdlib.h>
static CRITICAL_SECTION dtoa_CritSec[NLOCKS];
static long dtoa_CS_init = 0;
/*
   1 = initializing
   2 = initialized
   3 = deleted
*/
static void dtoa_lock_cleanup (void)
{
      long last_CS_init = InterlockedExchange (&dtoa_CS_init,3);
      if (2 == last_CS_init) {
            int i;
            for (i = 0; i < NLOCKS; i++)
                  DeleteCriticalSection (&dtoa_CritSec[i]);
      }
}

static void dtoa_lock (int n)
{
      if (2 == dtoa_CS_init) {
            EnterCriticalSection (&dtoa_CritSec[n]);
            return;
      }
      else if (0 == dtoa_CS_init) {
            long last_CS_init = InterlockedExchange (&dtoa_CS_init, 1);
            if (0 == last_CS_init) {
                  int i;
                  for (i = 0; i < NLOCKS;  i++)
                        InitializeCriticalSection (&dtoa_CritSec[i]);
                  atexit (dtoa_lock_cleanup);
                  dtoa_CS_init = 2;
            }
            else if (2 == last_CS_init)
                  dtoa_CS_init = 2;
      }
      /*  Another thread is initializing. Wait. */
      while (1 == dtoa_CS_init)
            Sleep (1);

      /* It had better be initialized now. */
      if (2 == dtoa_CS_init)
            EnterCriticalSection(&dtoa_CritSec[n]);
}

static void dtoa_unlock (int n)
{
      if (2 == dtoa_CS_init)
            LeaveCriticalSection (&dtoa_CritSec[n]);
}

#define ACQUIRE_DTOA_LOCK(n) dtoa_lock(n)
#define FREE_DTOA_LOCK(n) dtoa_unlock(n)
#endif      /* USE_WIN32_SL */

#endif      /* __MINGW32__ / __MINGW64__ */

#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 (int k)
{
      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
                && (size_t) (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 (Bigint *v)
{
      if (v) {
            if (v->k > Kmax)
                  free((void*)v);
            else {
                  ACQUIRE_DTOA_LOCK(0);
                  v->next = freelist[v->k];
                  freelist[v->k] = v;
                  FREE_DTOA_LOCK(0);
            }
      }
}

/* lo0bits():  Shift y so lowest bit is 1 and return the
 *           number of bits y was shifted.
 * With GCC, we use an inline wrapper for __builtin_clz()
 */
#ifndef __GNUC__
int lo0bits (ULong *y)
{
      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;
}
#endif      /* __GNUC__ */

Bigint *multadd (Bigint *b, int m, int a) /* multiply by m and add a */
{
      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;
}

/* hi0bits(); 
 * With GCC, we use an inline wrapper for __builtin_clz()
 */
#ifndef __GNUC__
int hi0bits_D2A (ULong x)
{
      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;
}
#endif      /* __GNUC__ */

Bigint *i2b (int i)
{
      Bigint *b;

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

Bigint *mult (Bigint *a, Bigint *b)
{
      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 (Bigint *b, int k)
{
      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 (Bigint *b, int k)
{
      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 (Bigint *a, Bigint *b)
{
      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 (Bigint *a, Bigint *b)
{
      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 (Bigint *a, int *e)
{
      ULong *xa, *xa0, w, y, z;
      int k;
      union _dbl_union d;
#define d0 word0(&d)
#define d1 word1(&d)

      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:
      return dval(&d);
#undef d0
#undef d1
}

Bigint *d2b (double dd, int *e, int *bits)
{
      Bigint *b;
      union _dbl_union d;
#ifndef Sudden_Underflow
      int i;
#endif
      int de, k;
      ULong *x, y, z;
#define d0 word0(&d)
#define d1 word1(&d)
      d.d = dd;

#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);
      z |= Exp_msk11;
#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
            *e = de - Bias - (P-1) + k;
            *bits = P - k;
#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
bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };

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

char *strcp_D2A (char *a, const char *b)
{
      while((*a = *b++))
            a++;
      return a;
}

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


Generated by  Doxygen 1.6.0   Back to index