home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch1_4 / rat.c next >
C/C++ Source or Header  |  1995-04-04  |  6KB  |  135 lines
  1. /****** rat.c ******/
  2. /* Ken Shoemake, 1994 */
  3.  
  4. #include <math.h>
  5. #include "rat.h"
  6.  
  7. static void Mul32(UINT32 x, UINT32 y, UINT32 *hi, UINT32 *lo)
  8. {
  9.     UINT32 xlo = x&0xFFFF, xhi = (x>>16)&0xFFFF;
  10.     UINT32 ylo = y&0xFFFF, yhi = (y>>16)&0xFFFF;
  11.     UINT32 t1, t2, t3;
  12.     UINT32 lolo, lohi, t1lo, t1hi, t2lo, t2hi, carry;
  13.     *lo = xlo * ylo; *hi = xhi * yhi;
  14.     t1 = xhi * ylo; t2 = xlo * yhi;
  15.     lolo = *lo&0xFFFF; lohi = (*lo>>16)&0xFFFF;
  16.     t1lo = t1&0xFFFF; t1hi = (t1>>16)&0xFFFF;
  17.     t2lo = t2&0xFFFF; t2hi = (t2>>16)&0xFFFF;
  18.     t3 = lohi + t1lo + t2lo;
  19.     carry = t3&0xFFFF; lohi = (t3<<16)&0xFFFF;
  20.     *hi += t1hi + t2hi + carry; *lo = (lohi<<16) + lolo;
  21. }
  22.  
  23. /* ratapprox(x,n) returns the best rational approximation to x whose numerator
  24.     and denominator are less than or equal to n in absolute value. The denominator
  25.     will be positive, and the numerator and denominator will be in lowest terms.
  26.     IEEE 32-bit floating point and 32-bit integers are required.
  27.    All best rational approximations of a real x may be obtained from x's
  28.     continued fraction representation, x = c0 + 1/(c1 + 1/(c2 + 1/(...)))
  29.     by truncation to k terms and possibly "interpolation" of the last term.
  30.     The continued fraction expansion itself is obtained by a variant of the
  31.     standard GCD algorithm, which is folded into the recursions generating
  32.     successive numerators and denominators. These recursions both have the
  33.     same form: f[k] = c[k]*f[k-1] + f[k-2]. For further information, see
  34.     Fundamentals of Mathematics, Volume I, MIT Press, 1983.
  35.  */
  36. Rational ratapprox(float x, INT32 limit)
  37. {
  38.     float tooLargeToFix = ldexp(1.0, BITS);         /* 0x4f000000=2147483648.0 */
  39.     float tooSmallToFix = ldexp(1.0, -BITS);        /* 0x30000000=4.6566e-10 */
  40.     float halfTooSmallToFix = ldexp(1.0, -BITS-1);  /* 0x2f800000=2.3283e-10 */
  41.     int expForInt = 24;     /* This exponent in float makes mantissa an INT32 */
  42.     static Rational ratZero = {0, 1};
  43.     INT32 sign = 1;
  44.     BOOL flip = FALSE;      /* If TRUE, nk and dk are swapped */
  45.     int scale;              /* Power of 2 to get x into integer domain */
  46.     UINT32 ak2, ak1, ak;    /* GCD arguments, initially 1 and x */
  47.     UINT32 ck, climit;      /* ck is GCD quotient and c.f. term k */
  48.     INT32 nk, dk;           /* Result num. and den., recursively found */
  49.     INT32 nk1 = 0, dk2 = 0; /* History terms for recursion */
  50.     INT32 nk2 = 1, dk1 = 1;
  51.     BOOL hard = FALSE;
  52.     Rational val;
  53.     
  54.     if (limit <= 0) return (ratZero);       /* Insist limit > 0 */
  55.     if (x < 0.0) {x = -x; sign = -1;}
  56.     val.numer = sign; val.denom = limit;
  57.     /* Handle first non-zero term of continued fraction,
  58.         rest prepared for integer GCD, sure to fit.
  59.      */
  60.     if (x >= 1.0) {/* First continued fraction term is non-zero */
  61.             float rest;
  62.             if (x >= tooLargeToFix || (ck = x) >= limit)
  63.                 {val.numer = sign*limit; val.denom = 1; return (val);}
  64.             flip = TRUE;        /* Keep denominator larger, for fast loop test */
  65.             nk = 1;  dk = ck;   /* Make new numerator and denominator */
  66.             rest = x - ck;
  67.             frexp(1.0,&scale);
  68.             scale = expForInt - scale;
  69.             ak = ldexp(rest, scale);
  70.             ak1 = ldexp(1.0, scale);
  71.         } else {/* First continued fraction term is zero */
  72.             int n;
  73.             UINT32 num = 1;
  74.             if (x <= tooSmallToFix) {       /* Is x too tiny to be 1/INT32 ? */
  75.                 if (x <= halfTooSmallToFix) return (ratZero);
  76.                 if (limit > (UINT32)(0.5/x)) return (val);
  77.                 else return (ratZero);
  78.             }
  79.             /* Treating 1.0 and x as integers, divide 1/x in a peculiar way
  80.                 to get accurate remainder
  81.              */
  82.             frexp(x,&scale);
  83.             scale = expForInt - scale;
  84.             ak1 = ldexp(x, scale);
  85.             n = (scale<BITS)?scale:BITS;    /* Stay within UINT32 arithmetic */
  86.             num <<= n;
  87.             ck = num/ak1;                   /* First attempt at 1/x */
  88.             ak = num%ak1;                   /* First attempt at remainder */
  89.             while ((scale -= n) > 0) {/* Shift quotient, remainder until done */
  90.                 n = (scale<8)?scale:8;      /* The 8 is 24 bits of x in 32 bits */
  91.                 num = ak<<n;
  92.                 ck = ck<<n + num/ak1;
  93.                 ak = num%ak1;               /* Reduce remainder */
  94.             }
  95.             /* All done with divide */
  96.             if (ck >= limit) {              /* Is x too tiny to be 1/limit ? */
  97.                 if (2*limit > ck)
  98.                     return (val);
  99.                 else return (ratZero);
  100.             }
  101.             nk = 1;  dk = ck;               /* Make new numer and denom */
  102.         }
  103.     while (ak != 0) {                   /* If possible, quit when have exact result */
  104.         ak2 = ak1;  ak1 = ak;           /* Prepare for next term */
  105.         nk2 = nk1;  nk1 = nk;           /* (This loop does almost all the work) */
  106.         dk2 = dk1;  dk1 = dk;
  107.         ck = ak2/ak1;                   /* Get next term of continued fraction */
  108.         ak = ak2 - ck*ak1;              /* Get remainder (GCD step) */
  109.         climit = (limit - dk2)/dk1;     /* Anticipate result of recursion on denom */
  110.         if (climit <= ck) {hard = TRUE; break;} /* Do not let denom exceed limit */
  111.         nk = ck*nk1 + nk2;              /* Make new result numer and denom */
  112.         dk = ck*dk1 + dk2;
  113.     }
  114.     if (hard) {
  115.         UINT32 twoClimit = 2*climit;
  116.         if (twoClimit >= ck) {      /* If climit < ck/2 no improvement possible */
  117.             nk = climit*nk1 + nk2;  /* Make limited numerator and denominator */
  118.             dk = climit*dk1 + dk2;
  119.             if (twoClimit == ck) {  /* If climit == ck improvement not sure */
  120.                 /* Using climit is better only when dk2/dk1 > ak/ak1 */
  121.                 /* For full precision, test dk2*ak1 > dk1*ak */
  122.                 UINT32 dk2ak1Hi, dk2ak1Lo, dk1akHi, dk1akLo;
  123.                 Mul32(flip?nk2:dk2, ak1, &dk2ak1Hi, &dk2ak1Lo);
  124.                 Mul32(flip?nk1:dk1, ak, &dk1akHi, &dk1akLo);
  125.                 if ((dk2ak1Hi < dk1akHi)
  126.                 || ((dk2ak1Hi == dk1akHi) && (dk2ak1Lo <= dk1akLo)))
  127.                     { nk = nk1;  dk = dk1; }    /* Not an improvement, so undo step */
  128.             }
  129.         }
  130.     }
  131.     if (flip) {val.numer = sign*dk;  val.denom = nk;}
  132.     else      {val.numer = sign*nk;  val.denom = dk;}
  133.     return (val);
  134. }
  135.