Floating-Point Printer Sample Code

• fp-free.ss: Scheme code for the fully general free-format algorithm for positive floating-point numbers

  ;;; Copyright (c) 1996 Robert G. Burger
  ;;;
  ;;; Permission is hereby granted, free of charge, to any person
  ;;; obtaining a copy of this software and associated documentation
  ;;; files (the "Software"), to deal in the Software without
  ;;; restriction, including without limitation the rights to use, copy,
  ;;; modify, merge, publish, distribute, sublicense, and/or sell copies
  ;;; of the Software, and to permit persons to whom the Software is
  ;;; furnished to do so, subject to the following conditions:
  ;;;
  ;;; The above copyright notice and this permission notice shall be
  ;;; included in all copies or substantial portions of the Software.
  ;;;
  ;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  ;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  ;;; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  ;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  ;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  ;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  ;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  ;;; DEALINGS IN THE SOFTWARE.
  
  ;;; Fully general free-format algorithm for positive floating-point numbers
  
  ;;; It uses the floating-point logarithm to estimate the scaling factor
  ;;; and a table to look up powers of ten.
  
  ;;; Input to flonum->digits:
  ;;;       v -- a positive floating-point number, f x b^e
  ;;;       f -- mantissa of v
  ;;;       e -- exponent of v
  ;;;   min-e -- smallest representable exponent
  ;;;       p -- fixed size of mantissa
  ;;;       b -- input base (usually two)
  ;;;      ob -- output base (usually ten)
  
  ;;; Output: (k d_1 d_2 ... d_n),
  ;;;   where 0.d_1...d_n x ob^k is the shortest correctly rounded base-ob
  ;;;   number that rounds to v when input (it does not assume any particular
  ;;;   input rounding algorithm)
  
  ;;; See also "Printing Floating-Point Numbers Quickly and Accurately"
  ;;; in Proceedings of the SIGPLAN '96 Conference on Programming Language
  ;;; Design and Implementation.
  
  ;;; Author: Bob Burger  Date: March 1996
  
  (define flonum->digits
    (lambda (v f e min-e p b ob)
      (if (>= e 0)
          (if (not (= f (expt b (- p 1))))
              (let ([be (expt b e)])
                (scale (* f be 2) 2 be be 0 ob #f #f v))
              (let* ([be (expt b e)] [be1 (* be b)])
                (scale (* f be1 2) (* b 2) be1 be 0 ob #f #f v)))
          (if (or (= e min-e) (not (= f (expt b (- p 1)))))
              (scale (* f 2) (* (expt b (- e)) 2) 1 1 0 ob #f #f v)
              (scale (* f b 2) (* (expt b (- 1 e)) 2)
                     b 1 0 ob #f #f v)))))
  
  (define scale
    (lambda (r s m+ m- k ob low-ok? high-ok? v)
      (let ([est (inexact->exact (ceiling (- (logB ob v) 1e-10)))])
        (if (>= est 0)
            (fixup r (* s (exptt ob est)) m+ m- est ob low-ok? high-ok?)
            (let ([scale (exptt ob (- est))])
              (fixup (* r scale) s (* m+ scale) (* m- scale)
                     est ob low-ok? high-ok?))))))
  
  (define fixup
    (lambda (r s m+ m- k ob low-ok? high-ok?)
      (if ((if high-ok? >= >) (+ r m+) s) ; too low?
          (cons (+ k 1) (generate r s m+ m- ob low-ok? high-ok?))
          (cons k
                (generate (* r ob) s (* m+ ob) (* m- ob) ob low-ok? high-ok?)))))
  
  (define generate
    (lambda (r s m+ m- ob low-ok? high-ok?)
      (let ([d (quotient r s)]
            [r (remainder r s)])
        (let ([tc1 ((if low-ok? <= <) r m-)]
              [tc2 ((if high-ok? >= >) (+ r m+) s)])
          (if (not tc1)
              (if (not tc2)
                  (cons d (generate (* r ob) s (* m+ ob) (* m- ob)
                                    ob low-ok? high-ok?))
                  (list (+ d 1)))
              (if (not tc2)
                  (list d)
                  (if (< (* r 2) s)
                      (list d)
                      (list (+ d 1)))))))))
  
  (define exptt
    (let ([table (make-vector 326)])
      (do ([k 0 (+ k 1)] [v 1 (* v 10)])
          ((= k 326))
        (vector-set! table k v))
      (lambda (B k)
        (if (and (= B 10) (<= 0 k 325))
            (vector-ref table k)
            (expt B k)))))
  
  (define logB
    (let ([table (make-vector 37)])
      (do ([B 2 (+ B 1)])
          ((= B 37))
        (vector-set! table B (/ (log B))))
      (lambda (B x)
        (if (<= 2 B 36)
            (* (log x) (vector-ref table B))
            (/ (log x) (log B))))))
  

• fp-free10.ss: Scheme code for the free-format algorithm for printing IEEE double-precision positive floating-point numbers in base 10

  ;;; Copyright (c) 1996 Robert G. Burger
  ;;;
  ;;; Permission is hereby granted, free of charge, to any person
  ;;; obtaining a copy of this software and associated documentation
  ;;; files (the "Software"), to deal in the Software without
  ;;; restriction, including without limitation the rights to use, copy,
  ;;; modify, merge, publish, distribute, sublicense, and/or sell copies
  ;;; of the Software, and to permit persons to whom the Software is
  ;;; furnished to do so, subject to the following conditions:
  ;;;
  ;;; The above copyright notice and this permission notice shall be
  ;;; included in all copies or substantial portions of the Software.
  ;;;
  ;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  ;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  ;;; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  ;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  ;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  ;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  ;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  ;;; DEALINGS IN THE SOFTWARE.
  
  ;;; Free-format algorithm for printing IEEE double-precision positive
  ;;; floating-point numbers in base 10
  
  ;;; It uses the floating-point logarithm to estimate the scaling factor
  ;;; and a table to look up powers of ten.
  
  ;;; Input to flonum->digits:
  ;;;       v -- a positive floating-point number, f x 2^e
  ;;;       f -- mantissa of v
  ;;;       e -- exponent of v
  
  ;;; Output: (k d_1 d_2 ... d_n),
  ;;;   where 0.d_1...d_n x 10^k is the shortest correctly rounded base-10
  ;;;   number that rounds to v when input (it assumes the input
  ;;;   routine rounds to even)
  
  ;;; See also "Printing Floating-Point Numbers Quickly and Accurately"
  ;;; in Proceedings of the SIGPLAN '96 Conference on Programming Language
  ;;; Design and Implementation.
  
  ;;; Author: Bob Burger  Date: March 1996
  
  (define flonum->digits
    (let ([min-e -1074]
          [bp-1 (expt 2 52)])
      (lambda (v f e)
        (let ([round? (even? f)])
          (if (>= e 0)
              (if (not (= f bp-1))
                  (let ([be (expt 2 e)])
                    (scale (* f be 2) 2 be be 0 round? round? v))
                  (let ([be (expt 2 e)])
                    (scale (* f be 4) 4 (* be 2) be 0 round? round? v)))
              (if (or (= e min-e) (not (= f bp-1)))
                  (scale (* f 2) (expt 2 (- 1 e)) 1 1 0 round? round? v)
                  (scale (* f 4) (expt 2 (- 2 e)) 2 1 0 round? round? v)))))))
  
  (define scale
    (lambda (r s m+ m- k low-ok? high-ok? v)
      (let ([est (inexact->exact (ceiling (- (log10 v) 1e-10)))])
        (if (>= est 0)
            (fixup r (* s (expt10 est)) m+ m- est low-ok? high-ok?)
            (let ([scale (expt10 (- est))])
              (fixup (* r scale) s (* m+ scale) (* m- scale)
                     est low-ok? high-ok?))))))
  
  (define fixup
    (lambda (r s m+ m- k low-ok? high-ok?)
      (if ((if high-ok? >= >) (+ r m+) s) ; too low?
          (cons (+ k 1) (generate r s m+ m- low-ok? high-ok?))
          (cons k
                (generate (* r 10) s (* m+ 10) (* m- 10) low-ok? high-ok?)))))
  
  (define generate
    (lambda (r s m+ m- low-ok? high-ok?)
      (let ([d (quotient r s)]
            [r (remainder r s)])
        (let ([tc1 ((if low-ok? <= <) r m-)]
              [tc2 ((if high-ok? >= >) (+ r m+) s)])
          (if (not tc1)
              (if (not tc2)
                  (cons d (generate (* r 10) s (* m+ 10) (* m- 10)
                                    low-ok? high-ok?))
                  (list (+ d 1)))
              (if (not tc2)
                  (list d)
                  (if (< (* r 2) s)
                      (list d)
                      (list (+ d 1)))))))))
  
  (define expt10
    (let ([table (make-vector 326)])
      (do ([k 0 (+ k 1)] [v 1 (* v 10)])
          ((= k 326))
        (vector-set! table k v))
      (lambda (k)
        (vector-ref table k))))
  
  (define log10
    (let ([f (/ (log 10))])
      (lambda (x)
        (* (log x) f))))
  

• free.c: C code for printing IEEE double-precision positive floating-point numbers in free-format base 10

  /* Copyright (c) 1996 Robert G. Burger
  
   * Permission is hereby granted, free of charge, to any person
   * obtaining a copy of this software and associated documentation
   * files (the "Software"), to deal in the Software without
   * restriction, including without limitation the rights to use, copy,
   * modify, merge, publish, distribute, sublicense, and/or sell copies
   * of the Software, and to permit persons to whom the Software is
   * furnished to do so, subject to the following conditions:
  
   * The above copyright notice and this permission notice shall be
   * included in all copies or substantial portions of the Software.
  
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
   * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
   * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
   * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   * SOFTWARE.
   */
  
  /* Free-format floating point 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;
  
  Bignum R, S, MP, MM, five[MAX_FIVE];
  Bignum S2, S3, S4, S5, S6, S7, S8, S9;
  int ruf, k, s_n, use_mp, qr_shift, sl, slr;
  
  void mul10 PROTO((Bignum *x));
  void big_short_mul PROTO((Bignum *x, Bigit y, Bignum *z));
  void print_big PROTO((Bignum *x));
  int estimate PROTO((int n));
  void one_shift_left PROTO((int y, Bignum *z));
  void short_shift_left PROTO((Bigit x, int y, Bignum *z));
  void big_shift_left PROTO((Bignum *x, int y, Bignum *z));
  int big_comp PROTO((Bignum *x, Bignum *y));
  int sub_big PROTO((Bignum *x, Bignum *y, Bignum *z));
  void add_big PROTO((Bignum *x, Bignum *y, Bignum *z));
  int add_cmp PROTO((void));
  int qr PROTO((void));
  
  int dragon PROTO((char *buf, double v));
  void free_init PROTO((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));\
  }
  
  void mul10(x) 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;
  }
  
  void big_short_mul(x, y, z) Bignum *x, *z; Bigit y; {
     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;
  }
  
  void print_big(x) 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));
     }
  }
  
  int estimate(n) int n; {
     if (n < 0)
        return (int)(n*0.3010299956639812);
     else
        return 1+(int)(n*0.3010299956639811);
  }
  
  void one_shift_left(y, z) 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;
  }
  
  void short_shift_left(x, y, z) 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;
  }
  
  void big_shift_left(x, y, z) Bignum *x, *z; int y; {
     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;
  }
  
  
  int big_comp(x, y) Bignum *x, *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;
  }
  
  int sub_big(x, y, z) Bignum *x, *y, *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 (zl > 0 && *--zp == 0) zl--;
    z->l = zl;
    return 0;
  }
  
  void add_big(x, y, z) Bignum *x, *y, *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;
  }
  
  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);
  }
  
  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 dragon(buf, v) char *buf; double v; {
     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;
  
     if (sign) *buf++ = '-';
     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 free_init() {
     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);
     */
  }
  

• fixed.c: C code for a quick hack of free.c for printing fixed-format IEEE double-precision floating-point numbers in base 10

  /* Copyright (c) 1996 Robert G. Burger
  
   * Permission is hereby granted, free of charge, to any person
   * obtaining a copy of this software and associated documentation
   * files (the "Software"), to deal in the Software without
   * restriction, including without limitation the rights to use, copy,
   * modify, merge, publish, distribute, sublicense, and/or sell copies
   * of the Software, and to permit persons to whom the Software is
   * furnished to do so, subject to the following conditions:
  
   * The above copyright notice and this permission notice shall be
   * included in all copies or substantial portions of the Software.
  
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
   * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
   * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
   * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   * SOFTWARE.
   */
  
  /* 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;
  }
  

• fp.h: C header file for free.c and fixed.c

  /* Copyright (c) 1996 Robert G. Burger
  
   * Permission is hereby granted, free of charge, to any person
   * obtaining a copy of this software and associated documentation
   * files (the "Software"), to deal in the Software without
   * restriction, including without limitation the rights to use, copy,
   * modify, merge, publish, distribute, sublicense, and/or sell copies
   * of the Software, and to permit persons to whom the Software is
   * furnished to do so, subject to the following conditions:
  
   * The above copyright notice and this permission notice shall be
   * included in all copies or substantial portions of the Software.
  
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
   * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
   * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
   * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   * SOFTWARE.
   */
  
  #if defined(__alpha) || defined(__i386) || defined(_M_IX86) || defined(__x86_64)
  #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) || defined(__x86_64)
  #define UNSIGNED64 unsigned long
  #else
  #define UNSIGNED64 unsigned long long
  #endif
  
  #if defined(sparc) || defined(hpux)
  #define PROTO(x) ()
  #else
  #define PROTO(x) x
  #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