mprec.c

标准c库代码,可以应用于各个系统提供了大量的基本函数

	{
	  *x1++ = *x << k | z;
	  z = *x++ >> k1;
	}
      while (x < xe);
      if (*x1 = z)
	++n1;
    }
#else
  if (k &= 0xf)
    {
      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 (ptr, b);
  return b1;
}

int
_DEFUN (cmp, (a, b), _Bigint * a _AND _Bigint * b)
{
  unsigned long *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 *
_DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND
	_Bigint * a _AND _Bigint * b)
{
  _Bigint *c;
  int i, wa, wb;
  long borrow, y;		/* We need signed shifts here. */
  unsigned long *xa, *xae, *xb, *xbe, *xc;
#ifdef Pack_32
  long z;
#endif

  i = cmp (a, b);
  if (!i)
    {
      c = Balloc (ptr, 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 (ptr, 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 Pack_32
  do
    {
      y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
      borrow = y >> 16;
      Sign_Extend (borrow, y);
      z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
      borrow = z >> 16;
      Sign_Extend (borrow, z);
      Storeinc (xc, z, y);
    }
  while (xb < xbe);
  while (xa < xae)
    {
      y = (*xa & 0xffff) + borrow;
      borrow = y >> 16;
      Sign_Extend (borrow, y);
      z = (*xa++ >> 16) + borrow;
      borrow = z >> 16;
      Sign_Extend (borrow, z);
      Storeinc (xc, z, y);
    }
#else
  do
    {
      y = *xa++ - *xb++ + borrow;
      borrow = y >> 16;
      Sign_Extend (borrow, y);
      *xc++ = y & 0xffff;
    }
  while (xb < xbe);
  while (xa < xae)
    {
      y = *xa++ + borrow;
      borrow = y >> 16;
      Sign_Extend (borrow, y);
      *xc++ = y & 0xffff;
    }
#endif
  while (!*--xc)
    wa--;
  c->_wds = wa;
  return c;
}

double
_DEFUN (ulp, (x), double x)
{
  register long L;
  double a;

  L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
#ifndef Sudden_Underflow
  if (L > 0)
    {
#endif
#ifdef IBM
      L |= Exp_msk1 >> 4;
#endif
      word0 (a) = L;
      word1 (a) = 0;
#ifndef Sudden_Underflow
    }
  else
    {
      L = -L >> Exp_shift;
      if (L < Exp_shift)
	{
	  word0 (a) = 0x80000 >> L;
	  word1 (a) = 0;
	}
      else
	{
	  word0 (a) = 0;
	  L -= Exp_shift;
	  word1 (a) = L >= 31 ? 1 : 1 << 31 - L;
	}
    }
#endif
  return a;
}

double
_DEFUN (b2d, (a, e),
	_Bigint * a _AND int *e)
{
  unsigned long *xa, *xa0, w, y, z;
  int k;
  double d;
#ifdef VAX
  unsigned long 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;
#else
#undef d0
#undef d1
#endif
  return d;
}

_Bigint *
_DEFUN (d2b,
	(ptr, d, e, bits),
	struct _reent * ptr _AND
	double d _AND
	int *e _AND
	int *bits)

{
  _Bigint *b;
  int de, i, k;
  unsigned long *x, y, z;
#ifdef VAX
  unsigned long 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 (ptr, 1);
#else
  b = Balloc (ptr, 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))
    z |= Exp_msk1;
#endif
#ifdef Pack_32
  if (y = d1)
    {
      if (k = lo0bits (&y))
	{
	  x[0] = y | z << 32 - k;
	  z >>= k;
	}
      else
	x[0] = y;
      i = b->_wds = (x[1] = z) ? 2 : 1;
    }
  else
    {
#ifdef DEBUG
      if (!z)
	Bug ("Zero passed to d2b");
#endif
      k = lo0bits (&z);
      x[0] = z;
      i = b->_wds = 1;
      k += 32;
    }
#else
  if (y = d1)
    {
      if (k = lo0bits (&y))
	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

double
_DEFUN (ratio, (a, b), _Bigint * a _AND _Bigint * b)

{
  double da, db;
  int k, ka, kb;

  da = b2d (a, &ka);
  db = b2d (b, &kb);
#ifdef Pack_32
  k = ka - kb + 32 * (a->_wds - b->_wds);
#else
  k = ka - kb + 16 * (a->_wds - b->_wds);
#endif
#ifdef IBM
  if (k > 0)
    {
      word0 (da) += (k >> 2) * Exp_msk1;
      if (k &= 3)
	da *= 1 << k;
    }
  else
    {
      k = -k;
      word0 (db) += (k >> 2) * Exp_msk1;
      if (k &= 3)
	db *= 1 << k;
    }
#else
  if (k > 0)
    word0 (da) += k * Exp_msk1;
  else
    {
      k = -k;
      word0 (db) += k * Exp_msk1;
    }
#endif
  return da / db;
}


_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, 1e23, 1e24

};

#if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
_CONST double bigtens[] =
{1e16, 1e32, 1e64, 1e128, 1e256};

_CONST double tinytens[] =
{1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
#else
_CONST double bigtens[] =
{1e16, 1e32};

_CONST double tinytens[] =
{1e-16, 1e-32};
#endif


double
_DEFUN (_mprec_log10, (dig),
	int dig)
{
  double v = 1.0;
  if (dig < 24)
    return tens[dig];
  while (dig > 0)
    {
      v *= 10;
      dig--;
    }
  return v;
}