text: /****************************************************************
   *
   * The author of this software is David M. Gay.
   *
   * Copyright (c) 1991 by AT&T.
   *
   * Permission to use, copy, modify, and distribute this software for any
   * purpose without fee is hereby granted, provided that this entire notice
   * is included in all copies of any software which is or includes a copy
   * or modification of this software and in all copies of the supporting
   * documentation for such software.
   *
   * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
   * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
   * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
   * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
   *
   ***************************************************************/
  
  /* Please send bug reports to
  	David M. Gay
  	AT&T Bell Laboratories, Room 2C-463
  	600 Mountain Avenue
  	Murray Hill, NJ 07974-2070
  	U.S.A.
  	dmg@research.att.com or research!dmg
   */
  
  /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
   *
   * This strtod returns a nearest machine number to the input decimal
   * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
   * broken by the IEEE round-even rule.  Otherwise ties are broken by
   * biased rounding (add half and chop).
   *
   * Inspired loosely by William D. Clinger's paper "How to Read Floating
   * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
   *
   * Modifications:
   *
   *	1. We only require IEEE, IBM, or VAX double-precision
   *		arithmetic (not IEEE double-extended).
   *	2. We get by with floating-point arithmetic in a case that
   *		Clinger missed -- when we're computing d * 10^n
   *		for a small integer d and the integer n is not too
   *		much larger than 22 (the maximum integer k for which
   *		we can represent 10^k exactly), we may be able to
   *		compute (d*10^k) * 10^(e-k) with just one roundoff.
   *	3. Rather than a bit-at-a-time adjustment of the binary
   *		result in the hard case, we use floating-point
   *		arithmetic to determine the adjustment to within
   *		one bit; only in really hard cases do we need to
   *		compute a second residual.
   *	4. Because of 3., we don't need a large table of powers of 10
   *		for ten-to-e (just some small tables, e.g. of 10^k
   *		for 0 <= k <= 22).
   */
  
  /*
   * #define IEEE_8087 for IEEE-arithmetic machines where the least
   *	significant byte has the lowest address.
   * #define IEEE_MC68k for IEEE-arithmetic machines where the most
   *	significant byte has the lowest address.
   * #define Sudden_Underflow for IEEE-format machines without gradual
   *	underflow (i.e., that flush to zero on underflow).
   * #define IBM for IBM mainframe-style floating-point arithmetic.
   * #define VAX for VAX-style floating-point arithmetic.
   * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
   * #define No_leftright to omit left-right logic in fast floating-point
   *	computation of dtoa.
   * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
   * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
   *	that use extended-precision instructions to compute rounded
   *	products and quotients) with IBM.
   * #define ROUND_BIASED for IEEE-format with biased rounding.
   * #define Inaccurate_Divide for IEEE-format with correctly rounded
   *	products but inaccurate quotients, e.g., for Intel i860.
   * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
   *	integer arithmetic.  Whether this speeds things up or slows things
   *	down depends on the machine and the number being converted.
   */
  
  #include <stdlib.h>
  #include <string.h>
  #include <java-assert.h>
  #include "mprec.h"
  
  /* reent.c knows this value */
  #define _Kmax 15
  #include <stdio.h>
  
  _Jv_Bigint *
  _DEFUN (Balloc, (ptr, k), struct _Jv_reent *ptr _AND int k)
  {
    _Jv_Bigint *rv = NULL;
  
    int i = 0;
    int j = 1;
  
    JvAssert ((1 << k) < MAX_BIGNUM_WDS);
  
    while ((ptr->_allocation_map & j) && i < MAX_BIGNUMS)
      i++, j <<= 1;
  
    JvAssert (i < MAX_BIGNUMS);
  
    if (i >= MAX_BIGNUMS) 
      return NULL;
  
    ptr->_allocation_map |= j;
    rv = &ptr->_freelist[i];
        
    rv->_k = k;
    rv->_maxwds = 32;
  
    return rv;
  }
  
  
  void
  _DEFUN (Bfree, (ptr, v), struct _Jv_reent *ptr _AND _Jv_Bigint * v)
  {
    long i;
  
    i = v - ptr->_freelist;
  
    JvAssert (i >= 0 && i < MAX_BIGNUMS);
  
    if (i >= 0 && i < MAX_BIGNUMS)
      ptr->_allocation_map &= ~ (1 << i);
  }
  
  
  _Jv_Bigint *
  _DEFUN (multadd, (ptr, b, m, a),
  	struct _Jv_reent *ptr _AND
  	_Jv_Bigint * b _AND
  	int m _AND
  	int a)
  {
    int i, wds;
    unsigned long *x, y;
  #ifdef Pack_32
    unsigned long xi, z;
  #endif
    _Jv_Bigint *b1;
  
    wds = b->_wds;
    x = b->_x;
    i = 0;
    do
      {
  #ifdef Pack_32
        xi = *x;
        y = (xi & 0xffff) * m + a;
        z = (xi >> 16) * m + (y >> 16);
        a = (int) (z >> 16);
        *x++ = (z << 16) + (y & 0xffff);
  #else
        y = *x * m + a;
        a = (int) (y >> 16);
        *x++ = y & 0xffff;
  #endif
      }
    while (++i < wds);
    if (a)
      {
        if (wds >= b->_maxwds)
  	{
  	  b1 = Balloc (ptr, b->_k + 1);
  	  Bcopy (b1, b);
  	  Bfree (ptr, b);
  	  b = b1;
  	}
        b->_x[wds++] = a;
        b->_wds = wds;
      }
    return b;
  }
  
  _Jv_Bigint *
  _DEFUN (s2b, (ptr, s, nd0, nd, y9),
  	struct _Jv_reent * ptr _AND
  	_CONST char *s _AND
  	int nd0 _AND
  	int nd _AND
  	unsigned long y9)
  {
    _Jv_Bigint *b;
    int i, k;
    long x, y;
  
    x = (nd + 8) / 9;
    for (k = 0, y = 1; x > y; y <<= 1, k++);
  #ifdef Pack_32
    b = Balloc (ptr, k);
    b->_x[0] = y9;
    b->_wds = 1;
  #else
    b = Balloc (ptr, k + 1);
    b->_x[0] = y9 & 0xffff;
    b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
  #endif
  
    i = 9;
    if (9 < nd0)
      {
        s += 9;
        do
  	b = multadd (ptr, b, 10, *s++ - '0');
        while (++i < nd0);
        s++;
      }
    else
      s += 10;
    for (; i < nd; i++)
      b = multadd (ptr, b, 10, *s++ - '0');
    return b;
  }
  
  int
  _DEFUN (hi0bits,
  	(x), register unsigned long x)
  {
    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;
  }
  
  int
  _DEFUN (lo0bits, (y), unsigned long *y)
  {
    register int k;
    register unsigned long 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 & 1))
  	return 32;
      }
    *y = x;
    return k;
  }
  
  _Jv_Bigint *
  _DEFUN (i2b, (ptr, i), struct _Jv_reent * ptr _AND int i)
  {
    _Jv_Bigint *b;
  
    b = Balloc (ptr, 1);
    b->_x[0] = i;
    b->_wds = 1;
    return b;
  }
  
  _Jv_Bigint *
  _DEFUN (mult, (ptr, a, b), struct _Jv_reent * ptr _AND _Jv_Bigint * a _AND _Jv_Bigint * b)
  {
    _Jv_Bigint *c;
    int k, wa, wb, wc;
    unsigned long carry, y, z;
    unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  #ifdef Pack_32
    unsigned long z2;
  #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 (ptr, 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 Pack_32
    for (; xb < xbe; xb++, xc0++)
      {
        if ((y = *xb & 0xffff))
  	{
  	  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))
  	{
  	  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++)
  	{
  	  x = xa;
  	  xc = xc0;
  	  carry = 0;
  	  do
  	    {
  	      z = *x++ * y + *xc + carry;
  	      carry = z >> 16;
  	      *xc++ = z & 0xffff;
  	    }
  	  while (x < xae);
  	  *xc = carry;
  	}
      }
  #endif
    for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
    c->_wds = wc;
    return c;
  }
  
  _Jv_Bigint *
  _DEFUN (pow5mult,
  	(ptr, b, k), struct _Jv_reent * ptr _AND _Jv_Bigint * b _AND int k)
  {
    _Jv_Bigint *b1, *p5, *p51;
    int i;
    static _CONST int p05[3] = {5, 25, 125};
  
    if ((i = k & 3))
      b = multadd (ptr, b, p05[i - 1], 0);
  
    if (!(k >>= 2))
      return b;
    if (!(p5 = ptr->_p5s))
      {
        /* first time */
        p5 = ptr->_p5s = i2b (ptr, 625);
        p5->_next = 0;
      }
    for (;;)
      {
        if (k & 1)
  	{
  	  b1 = mult (ptr, b, p5);
  	  Bfree (ptr, b);
  	  b = b1;
  	}
        if (!(k >>= 1))
  	break;
        if (!(p51 = p5->_next))
  	{
  	  p51 = p5->_next = mult (ptr, p5, p5);
  	  p51->_next = 0;
  	}
        p5 = p51;
      }
    return b;
  }
  
  _Jv_Bigint *
  _DEFUN (lshift, (ptr, b, k), struct _Jv_reent * ptr _AND _Jv_Bigint * b _AND int k)
  {
    int i, k1, n, n1;
    _Jv_Bigint *b1;
    unsigned long *x, *x1, *xe, z;
  
  #ifdef Pack_32
    n = k >> 5;
  #else
    n = k >> 4;
  #endif
    k1 = b->_k;
    n1 = n + b->_wds + 1;
    for (i = b->_maxwds; n1 > i; i <<= 1)
      k1++;
    b1 = Balloc (ptr, k1);
    x1 = b1->_x;
    for (i = 0; i < n; i++)
      *x1++ = 0;
    x = b->_x;
    xe = x + b->_wds;
  #ifdef Pack_32
    if (k &= 0x1f)
      {
        k1 = 32 - k;
        z = 0;
        do
  	{
  	  *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), _Jv_Bigint * a _AND _Jv_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;
  }
  
  _Jv_Bigint *
  _DEFUN (diff, (ptr, a, b), struct _Jv_reent * ptr _AND
  	_Jv_Bigint * a _AND _Jv_Bigint * b)
  {
    _Jv_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)
  {
    union double_union x, a;
    register long L;
  
    x.d = _x;
  
    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;
  #ifndef _DOUBLE_IS_32BITS
        word1 (a) = 0;
  #endif
  
  #ifndef Sudden_Underflow
      }
    else
      {
        L = -L >> Exp_shift;
        if (L < Exp_shift)
  	{
  	  word0 (a) = 0x80000 >> L;
  #ifndef _DOUBLE_IS_32BITS
  	  word1 (a) = 0;
  #endif
  	}
        else
  	{
  	  word0 (a) = 0;
  	  L -= Exp_shift;
  #ifndef _DOUBLE_IS_32BITS
  	  word1 (a) = L >= 31 ? 1 : 1 << (31 - L);
  #endif
  	}
      }
  #endif
    return a.d;
  }
  
  double
  _DEFUN (b2d, (a, e),
  	_Jv_Bigint * a _AND int *e)
  {
    unsigned long *xa, *xa0, w, y, z;
    int k;
    union double_union 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;
  #ifndef _DOUBLE_IS_32BITS
        d1 = y << (32 - Ebits + k) | w >> (Ebits - k);
  #endif
        goto ret_d;
      }
    z = xa > xa0 ? *--xa : 0;
    if (k -= Ebits)
      {
        d0 = Exp_1 | y << k | z >> (32 - k);
        y = xa > xa0 ? *--xa : 0;
  #ifndef _DOUBLE_IS_32BITS
        d1 = z << k | y >> (32 - k);
  #endif
      }
    else
      {
        d0 = Exp_1 | y;
  #ifndef _DOUBLE_IS_32BITS
        d1 = z;
  #endif
      }
  #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.d;
  }
  
  _Jv_Bigint *
  _DEFUN (d2b,
  	(ptr, _d, e, bits),
  	struct _Jv_reent * ptr _AND
  	double _d _AND
  	int *e _AND
  	int *bits)
  
  {
    union double_union d;
    _Jv_Bigint *b;
    int de, i, k;
    unsigned long *x, y, z;
  #ifdef VAX
    unsigned long d0, d1;
    d.d = _d;
    d0 = word0 (d) >> 16 | word0 (d) << 16;
    d1 = word1 (d) >> 16 | word1 (d) << 16;
  #else
  #define d0 word0(d)
  #define d1 word1(d)
    d.d = _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
  #ifndef _DOUBLE_IS_32BITS
    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
  #endif
      {
  #ifdef DEBUG
        if (!z)
  	Bug ("Zero passed to d2b");
  #endif
        k = lo0bits (&z);
        x[0] = z;
        i = b->_wds = 1;
  #ifndef _DOUBLE_IS_32BITS
        k += 32;
  #endif
      }
  #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), _Jv_Bigint * a _AND _Jv_Bigint * b)
  
  {
    union double_union da, db;
    int k, ka, kb;
  
    da.d = b2d (a, &ka);
    db.d = 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.d *= 1 << k;
      }
    else
      {
        k = -k;
        word0 (db) += (k >> 2) * Exp_msk1;
        if (k &= 3)
  	db.d *= 1 << k;
      }
  #else
    if (k > 0)
      word0 (da) += k * Exp_msk1;
    else
      {
        k = -k;
        word0 (db) += k * Exp_msk1;
      }
  #endif
    return da.d / db.d;
  }
  
  
  _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
  
  
  