Objects/doubledigits.c

/* Free-format floating point printer
 * 
 * Based on "Floating-Point Printer Sample Code", by Robert G. Burger,
 * http://www.cs.indiana.edu/~burger/fp/index.html
 */

#include "Python.h"

#if defined(__alpha) || defined(__i386) || defined(_M_IX86) || defined(_M_X64) || defined(_M_IA64)
#define LITTLE_ENDIAN_IEEE_DOUBLE
#elif !(defined(__ppc__) || defined(sparc) || defined(__sgi) || defined(_IBMR2) || defined(hpux))
#error unknown machine type
#endif

#if defined(_M_IX86)
#define UNSIGNED64 unsigned __int64
#elif defined(__alpha)
#define UNSIGNED64 unsigned long
#else
#define UNSIGNED64 unsigned long long
#endif

#ifndef U32
#define U32 unsigned int
#endif

/* exponent stored + 1024, hidden bit to left of decimal point */
#define bias 1023
#define bitstoright 52
#define m1mask 0xf
#define hidden_bit 0x100000
#ifdef LITTLE_ENDIAN_IEEE_DOUBLE
struct dblflt {
    unsigned int m4: 16;
    unsigned int m3: 16;
    unsigned int m2: 16;
    unsigned int m1: 4;
    unsigned int e: 11;
    unsigned int s: 1;
};
#else
/* Big Endian IEEE Double Floats */
struct dblflt {
    unsigned int s: 1;
    unsigned int e: 11;
    unsigned int m1: 4;
    unsigned int m2: 16;
    unsigned int m3: 16;
    unsigned int m4: 16;
};
#endif
#define float_radix 2.147483648e9


typedef UNSIGNED64 Bigit;
#define BIGSIZE 24
#define MIN_E -1074
#define MAX_FIVE 325
#define B_P1 ((Bigit)1 << 52)

typedef struct {
   int l;
   Bigit d[BIGSIZE];
} Bignum;

static Bignum R, S, MP, MM, five[MAX_FIVE];
static Bignum S2, S3, S4, S5, S6, S7, S8, S9;
static int ruf, k, s_n, use_mp, qr_shift, sl, slr;

static void mul10(Bignum *x);
static void big_short_mul(Bignum *x, Bigit y, Bignum *z);
/*
static void print_big(Bignum *x);
*/
static int estimate(int n);
static void one_shift_left(int y, Bignum *z);
static void short_shift_left(Bigit x, int y, Bignum *z);
static void big_shift_left(Bignum *x, int y, Bignum *z);
static int big_comp(Bignum *x, Bignum *y);
static int sub_big(Bignum *x, Bignum *y, Bignum *z);
static void add_big(Bignum *x, Bignum *y, Bignum *z);
static int add_cmp(void);
static int qr(void);

/*static int _PyFloat_Digits(char *buf, double v, int *signum);*/
/*static void _PyFloat_DigitsInit(void);*/

#define ADD(x, y, z, k) {\
   Bigit x_add, z_add;\
   x_add = (x);\
   if ((k))\
     z_add = x_add + (y) + 1, (k) = (z_add <= x_add);\
   else\
     z_add = x_add + (y), (k) = (z_add < x_add);\
   (z) = z_add;\
}

#define SUB(x, y, z, b) {\
   Bigit x_sub, y_sub;\
   x_sub = (x); y_sub = (y);\
   if ((b))\
     (z) = x_sub - y_sub - 1, b = (y_sub >= x_sub);\
   else\
     (z) = x_sub - y_sub, b = (y_sub > x_sub);\
}

#define MUL(x, y, z, k) {\
   Bigit x_mul, low, high;\
   x_mul = (x);\
   low = (x_mul & 0xffffffff) * (y) + (k);\
   high = (x_mul >> 32) * (y) + (low >> 32);\
   (k) = high >> 32;\
   (z) = (low & 0xffffffff) | (high << 32);\
}

#define SLL(x, y, z, k) {\
   Bigit x_sll = (x);\
   (z) = (x_sll << (y)) | (k);\
   (k) = x_sll >> (64 - (y));\
}

static void
mul10(Bignum *x)
{
   int i, l;
   Bigit *p, k;

   l = x->l;
   for (i = l, p = &x->d[0], k = 0; i >= 0; i--)
      MUL(*p, 10, *p++, k);
   if (k != 0)
      *p = k, x->l = l+1;
}

static void
big_short_mul(Bignum *x, Bigit y, Bignum *z)
{
   int i, xl, zl;
   Bigit *xp, *zp, k;
   U32 high, low;

   xl = x->l;
   xp = &x->d[0];
   zl = xl;
   zp = &z->d[0];
   high = y >> 32;
   low = y & 0xffffffff;
   for (i = xl, k = 0; i >= 0; i--, xp++, zp++) {
      Bigit xlow, xhigh, z0, t, c, z1;
      xlow = *xp & 0xffffffff;
      xhigh = *xp >> 32;
      z0 = (xlow * low) + k; /* Cout is (z0 < k) */
      t = xhigh * low;
      z1 = (xlow * high) + t;
      c = (z1 < t);
      t = z0 >> 32;
      z1 += t;
      c += (z1 < t);
      *zp = (z1 << 32) | (z0 & 0xffffffff);
      k = (xhigh * high) + (c << 32) + (z1 >> 32) + (z0 < k);
   }
   if (k != 0)
      *zp = k, zl++;
   z->l = zl;
}

/*
static void
print_big(Bignum *x)
{
   int i;
   Bigit *p;

   printf("#x");
   i = x->l;
   p = &x->d[i];
   for (p = &x->d[i]; i >= 0; i--) {
      Bigit b = *p--;
      printf("%08x%08x", (int)(b >> 32), (int)(b & 0xffffffff));
   }
}
*/

static int
estimate(int n)
{
   if (n < 0)
      return (int)(n*0.3010299956639812);
   else
      return 1+(int)(n*0.3010299956639811);
}

static void
one_shift_left(int y, Bignum *z)
{
   int n, m, i;
   Bigit *zp;

   n = y / 64;
   m = y % 64;
   zp = &z->d[0];
   for (i = n; i > 0; i--) *zp++ = 0;
   *zp = (Bigit)1 << m;
   z->l = n;
}

static void
short_shift_left(Bigit x, int y, Bignum *z)
{
   int n, m, i, zl;
   Bigit *zp;
   
   n = y / 64;
   m = y % 64;
   zl = n;
   zp = &(z->d[0]);
   for (i = n; i > 0; i--) *zp++ = 0;
   if (m == 0)
      *zp = x;
   else {
      Bigit high = x >> (64 - m);
      *zp = x << m;
      if (high != 0)
         *++zp = high, zl++;
   }
   z->l = zl;
}

static void
big_shift_left(Bignum *x, int y, Bignum *z)
{
   int n, m, i, xl, zl;
   Bigit *xp, *zp, k;
   
   n = y / 64;
   m = y % 64;
   xl = x->l;
   xp = &(x->d[0]);
   zl = xl + n;
   zp = &(z->d[0]);
   for (i = n; i > 0; i--) *zp++ = 0;
   if (m == 0)
      for (i = xl; i >= 0; i--) *zp++ = *xp++;
   else {
      for (i = xl, k = 0; i >= 0; i--)
         SLL(*xp++, m, *zp++, k);
      if (k != 0)
         *zp = k, zl++;
   }
   z->l = zl;
}


static int
big_comp(Bignum *x, Bignum *y)
{
   int i, xl, yl;
   Bigit *xp, *yp;

   xl = x->l;
   yl = y->l;
   if (xl > yl) return 1;
   if (xl < yl) return -1;
   xp = &x->d[xl];
   yp = &y->d[xl];
   for (i = xl; i >= 0; i--, xp--, yp--) {
      Bigit a = *xp;
      Bigit b = *yp;

      if (a > b) return 1;
      else if (a < b) return -1;
   }
   return 0;
}

static int
sub_big(Bignum *x, Bignum *y, Bignum *z)
{
  int xl, yl, zl, b, i;
  Bigit *xp, *yp, *zp;

  xl = x->l;
  yl = y->l;
  if (yl > xl) return 1;
  xp = &x->d[0];
  yp = &y->d[0];
  zp = &z->d[0];
  
  for (i = yl, b = 0; i >= 0; i--)
    SUB(*xp++, *yp++, *zp++, b);
  for (i = xl-yl; b && i > 0; i--) {
    Bigit x_sub;
    x_sub = *xp++;
    *zp++ = x_sub - 1;
    b = (x_sub == 0);
  }
  for (; i > 0; i--) *zp++ = *xp++;
  if (b) return 1;
  zl = xl;
  while (*--zp == 0) zl--;
  z->l = zl;
  return 0;
}

static void
add_big(Bignum *x, Bignum *y, Bignum *z)
{
  int xl, yl, k, i;
  Bigit *xp, *yp, *zp;

  xl = x->l;
  yl = y->l;
  if (yl > xl) {
    int tl;
    Bignum *tn;
    tl = xl; xl = yl; yl = tl;
    tn = x; x = y; y = tn;
  }

  xp = &x->d[0];
  yp = &y->d[0];
  zp = &z->d[0];
  
  for (i = yl, k = 0; i >= 0; i--)
    ADD(*xp++, *yp++, *zp++, k);
  for (i = xl-yl; k && i > 0; i--) {
    Bigit z_add;
    z_add = *xp++ + 1;
    k = (z_add == 0);
    *zp++ = z_add;
  }
  for (; i > 0; i--) *zp++ = *xp++;
  if (k)
    *zp = 1, z->l = xl+1;
  else
    z->l = xl;
}

static int
add_cmp()
{
   int rl, ml, sl, suml;
   static Bignum sum;

   rl = R.l;
   ml = (use_mp ? MP.l : MM.l);
   sl = S.l;
   
   suml = rl >= ml ? rl : ml;
   if ((sl > suml+1) || ((sl == suml+1) && (S.d[sl] > 1))) return -1;
   if (sl < suml) return 1;

   add_big(&R, (use_mp ? &MP : &MM), &sum);
   return big_comp(&sum, &S);
}

static int
qr()
{
  if (big_comp(&R, &S5) < 0)
    if (big_comp(&R, &S2) < 0)
      if (big_comp(&R, &S) < 0)
        return 0;
      else {
        sub_big(&R, &S, &R);
        return 1;
      }
    else if (big_comp(&R, &S3) < 0) {
      sub_big(&R, &S2, &R);
      return 2;
    }
    else if (big_comp(&R, &S4) < 0) {
      sub_big(&R, &S3, &R);
      return 3;
    }
    else {
      sub_big(&R, &S4, &R);
      return 4;
    }
  else if (big_comp(&R, &S7) < 0)
    if (big_comp(&R, &S6) < 0) {
      sub_big(&R, &S5, &R);
      return 5;
    }
    else {
      sub_big(&R, &S6, &R);
      return 6;
    }
  else if (big_comp(&R, &S9) < 0)
    if (big_comp(&R, &S8) < 0) {
      sub_big(&R, &S7, &R);
      return 7;
    }
    else {
      sub_big(&R, &S8, &R);
      return 8;
    }
  else {
    sub_big(&R, &S9, &R);
    return 9;
  }
}

#define OUTDIG(d) { *buf++ = (d) + '0'; *buf = 0; return k; }

int
_PyFloat_Digits(char *buf, double v, int *signum)
{
   struct dblflt *x;
   int sign, e, f_n, m_n, i, d, tc1, tc2;
   Bigit f;

   /* decompose float into sign, mantissa & exponent */
   x = (struct dblflt *)&v;
   sign = x->s;
   e = x->e; 
   f = (Bigit)(x->m1 << 16 | x->m2) << 32 | (U32)(x->m3 << 16 | x->m4);
   if (e != 0) {
      e = e - bias - bitstoright;
      f |= (Bigit)hidden_bit << 32;
   }
   else if (f != 0)
      /* denormalized */
      e = 1 - bias - bitstoright;

   *signum = sign;
   if (f == 0) {
     *buf++ = '0';
     *buf = 0;
     return 0;
   }
   
   ruf = !(f & 1); /* ruf = (even? f) */

   /* Compute the scaling factor estimate, k */
   if (e > MIN_E)
      k = estimate(e+52);
   else {
      int n;
      Bigit y;
      
      for (n = e+52, y = (Bigit)1 << 52; f < y; n--) y >>= 1;
      k = estimate(n);
   }

   if (e >= 0)
      if (f != B_P1)
         use_mp = 0, f_n = e+1, s_n = 1, m_n = e;
      else
         use_mp = 1, f_n = e+2, s_n = 2, m_n = e;
   else
      if ((e == MIN_E) || (f != B_P1))
         use_mp = 0, f_n = 1, s_n = 1-e, m_n = 0;
      else
         use_mp = 1, f_n = 2, s_n = 2-e, m_n = 0;
   
   /* Scale it! */
   if (k == 0) {
      short_shift_left(f, f_n, &R);
      one_shift_left(s_n, &S);
      one_shift_left(m_n, &MM);
      if (use_mp) one_shift_left(m_n+1, &MP);
      qr_shift = 1;
   }
   else if (k > 0) {
      s_n += k;
      if (m_n >= s_n)
         f_n -= s_n, m_n -= s_n, s_n = 0;
      else 
         f_n -= m_n, s_n -= m_n, m_n = 0;
      short_shift_left(f, f_n, &R);
      big_shift_left(&five[k-1], s_n, &S);
      one_shift_left(m_n, &MM);
      if (use_mp) one_shift_left(m_n+1, &MP);
      qr_shift = 0;
   }
   else {
      Bignum *power = &five[-k-1];

      s_n += k;
      big_short_mul(power, f, &S);
      big_shift_left(&S, f_n, &R);
      one_shift_left(s_n, &S);
      big_shift_left(power, m_n, &MM);
      if (use_mp) big_shift_left(power, m_n+1, &MP);
      qr_shift = 1;
   }

   /* fixup */
   if (add_cmp() <= -ruf) {
     k--;
     mul10(&R);
     mul10(&MM);
     if (use_mp) mul10(&MP);
   }

   /*
   printf("\nk = %d\n", k);
   printf("R = "); print_big(&R);
   printf("\nS = "); print_big(&S);
   printf("\nM- = "); print_big(&MM);
   if (use_mp) printf("\nM+ = "), print_big(&MP);
   putchar('\n');
   fflush(0);
   */
   
   if (qr_shift) {
     sl = s_n / 64;
     slr = s_n % 64;
   }
   else {
     big_shift_left(&S, 1, &S2);
     add_big(&S2, &S, &S3);
     big_shift_left(&S2, 1, &S4);
     add_big(&S4, &S, &S5);
     add_big(&S4, &S2, &S6);
     add_big(&S4, &S3, &S7);
     big_shift_left(&S4, 1, &S8);
     add_big(&S8, &S, &S9);
   }

again:
   if (qr_shift) { /* Take advantage of the fact that S = (ash 1 s_n) */
      if (R.l < sl)
         d = 0;
      else if (R.l == sl) {
            Bigit *p;

            p = &R.d[sl];
            d = *p >> slr;
            *p &= ((Bigit)1 << slr) - 1;
            for (i = sl; (i > 0) && (*p == 0); i--) p--;
            R.l = i;
         }
      else {
         Bigit *p;
         
         p = &R.d[sl+1];
         d = *p << (64 - slr) | *(p-1) >> slr;
         p--;
         *p &= ((Bigit)1 << slr) - 1;
         for (i = sl; (i > 0) && (*p == 0); i--) p--;
         R.l = i;
      }
   }
   else /* We need to do quotient-remainder */
     d = qr();

   tc1 = big_comp(&R, &MM) < ruf;
   tc2 = add_cmp() > -ruf;
   if (!tc1)
      if (!tc2) {
         mul10(&R);
         mul10(&MM);
         if (use_mp) mul10(&MP);
         *buf++ = d + '0';
         goto again;
      }
      else
        OUTDIG(d+1)
   else
      if (!tc2)
         OUTDIG(d)
      else {
         big_shift_left(&R, 1, &MM);
         if (big_comp(&MM, &S) == -1)
            OUTDIG(d)
         else
            OUTDIG(d+1)
      }
}

void
_PyFloat_DigitsInit()
{
   int n, i, l;
   Bignum *b;
   Bigit *xp, *zp, k;

   five[0].l = l = 0;
   five[0].d[0] = 5;
   for (n = MAX_FIVE-1, b = &five[0]; n > 0; n--) {
      xp = &b->d[0];
      b++;
      zp = &b->d[0];
      for (i = l, k = 0; i >= 0; i--)
         MUL(*xp++, 5, *zp++, k);
      if (k != 0)
         *zp = k, l++;
      b->l = l;
   }

   /*
   for (n = 1, b = &five[0]; n <= MAX_FIVE; n++) {
      big_shift_left(b++, n, &R);
      print_big(&R);
      putchar('\n');
   }
   fflush(0);
   */
}