/* Fixed-format floating point printer
   A quick hack of the free-format printer */

#include "fp.h"

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;

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

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

int fixed PROTO((char *buf, double v, int prec));

int fixed(buf, v, prec) char *buf; double v; int prec; {
   struct dblflt *x;
   int sign, e, f_n, i, d, n;
   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;

   if (sign) *buf++ = '-';
   if (f == 0) {
     for (i = prec; i > 0; i--) *buf++ = '0';
     *buf = 0;
     return 0;
   }
   
   /* Compute the scaling factor estimate, k */
   if (e > MIN_E)
      k = estimate(e+52);
   else {
      Bigit y;
      
      for (n = e+52, y = (Bigit)1 << 52; f < y; n--) y >>= 1;
      k = estimate(n);
   }

   if (e >= 0)
      f_n = e, s_n = 0;
   else
      f_n = 0, s_n = -e;
   
   /* Scale it! */
   if (k == 0) {
      short_shift_left(f, f_n, &R);
      one_shift_left(s_n, &S);
      qr_shift = 1;
   }
   else if (k > 0) {
      s_n += k;
      if (f_n >= s_n)
         f_n -= s_n, s_n = 0;
      else 
         s_n -= f_n, f_n = 0;
      short_shift_left(f, f_n, &R);
      big_shift_left(&five[k-1], s_n, &S);
      qr_shift = 0;
   }
   else {
      s_n += k;
      big_short_mul(&five[-k-1], f, &S);
      big_shift_left(&S, f_n, &R);
      one_shift_left(s_n, &S);
      qr_shift = 1;
   }

   /* fixup */
   if (big_comp(&R, &S) < 0) k--, mul10(&R);

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

   for (n = prec; ;) {

      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();
      
      *buf++ = d + '0';
      if (--n == 0) break;
      mul10(&R);
   }

   big_shift_left(&R, 1, &MM);
   switch (big_comp(&MM, &S)) {
    case -1: /* No rounding needed */
      *buf = 0;
      return k;
    case 0: /* Exactly in the middle */
      *buf++ = '5';
      *buf = 0;
      return k;
    default:  /* Round up */
      *buf-- = 0;
      break;
   }

   for (n = prec; n > 0; n--) {
      char c;
      c = *buf;
      if (c != '9') {
         *buf = c+1;
         return k;
      }
      *buf-- = '0';
   }
   *++buf = '1';
   return k+1;
}