~ chicken-core (chicken-5) 18ce467a28f4f180b5efebdaeb167498c1b19508

commit 18ce467a28f4f180b5efebdaeb167498c1b19508
Author:     Peter Bex <peter.bex@xs4all.nl>
AuthorDate: Sun Sep 22 16:30:56 2013 +0200
Commit:     Christian Kellermann <ckeen@pestilenz.org>
CommitDate: Wed Oct 2 20:20:39 2013 +0200

    Fix #1051: use C99 isnormal() and return canned values.
    
    Instead of relying on unusable code from numbers egg, we return
    +1.0/+inf.0 or -1.0/+inf.0 when given a subnormal floating-point number.
    
    Signed-off-by: Christian Kellermann <ckeen@pestilenz.org>

diff --git a/runtime.c b/runtime.c
index 4ad05445..9bef9400 100644
--- a/runtime.c
+++ b/runtime.c
@@ -7371,46 +7371,27 @@ void C_ccall C_flonum_fraction(C_word c, C_word closure, C_word k, C_word n)
 
 void C_ccall C_flonum_rat(C_word c, C_word closure, C_word k, C_word n)
 {
-  double frac, tmp, numer, denom, factor, fn = C_flonum_magnitude(n);
-  double r1a, r1b;
+  double frac, tmp, numer, denom, fn = C_flonum_magnitude(n);
   double ga, gb;
   C_word ab[WORDS_PER_FLONUM * 2], *ap = ab;
   int i = 0;
 
-  if (n < 1 && n > -1) {
-    factor = pow(2, DBL_MANT_DIG);
-    fn *= factor;
-  } else {
-    factor = 1;
-  }
-
-  /* Calculate bit-length of the fractional part (ie, after decimal point) */
-  frac = fn;
-  while(!C_isnan(frac) && !C_isinf(frac) && C_modf(frac, &tmp) != 0.0) {
-    frac *= 2;
-    if (i++ > 3000) /* should this be flonum-maximum-exponent? */
-      barf(C_CANT_REPRESENT_INEXACT_ERROR, "fprat", n);
-  }
-
-  /* r1a and r1b are integral and form the rational number r1 = r1a/r1b. */
-  r1b = pow(2, i);
-  r1a = fn*r1b;
-
-  /*
-   * We "multiply" r1 with r2 given that r2 = 1/factor.
-   * result = (r1a * (factor / g)) / abs(r1b / g)   | g = gcd(r1b, factor)
-   */
-  ga = r1b;
-  gb = factor;
-  while(gb != 0.0) {
-    tmp = fmod(ga, gb);
-    ga = gb;
-    gb = tmp;
+  if (isnormal(fn)) {
+    /* Calculate bit-length of the fractional part (ie, after decimal point) */
+    frac = fn;
+    while(!C_isnan(frac) && !C_isinf(frac) && C_modf(frac, &tmp) != 0.0) {
+      frac *= 2;
+      if (i++ > 3000) /* should this be flonum-maximum-exponent? */
+        barf(C_CANT_REPRESENT_INEXACT_ERROR, "fprat", n);
+    }
+    
+    /* Now we can compute the rational number r = 2^i/X*2^i = numer/denom. */
+    denom = pow(2, i);
+    numer = fn*denom;
+  } else { /* denormalised/subnormal number: [+-]1.0/+inf.0 */
+    numer = fn > 0.0 ? 1.0 : -1.0;
+    denom = 1.0/0.0; /* +inf */
   }
-  /* ga now holds gcd(r1b, factor), and r1b and ga are absolute already */
-  numer = r1a * (factor / ga);
-  denom = r1b / ga;
-
   C_values(4, C_SCHEME_UNDEFINED, k, C_flonum(&ap, numer), C_flonum(&ap, denom));
 }
 
diff --git a/tests/library-tests.scm b/tests/library-tests.scm
index aaef6fd1..711341b4 100644
--- a/tests/library-tests.scm
+++ b/tests/library-tests.scm
@@ -86,6 +86,15 @@
 (assert (equal? 5.0 (numerator 1.25)))
 (assert (equal? 4.0 (denominator 1.25)))
 (assert (equal? -5.0 (numerator -1.25)))
+
+;; A few denormalised numbers, cribbed from NetBSD ATF tests for ldexp():
+(assert (equal? 1.0 (numerator 1.1125369292536006915451e-308)))
+(assert (equal? +inf.0 (denominator 1.1125369292536006915451e-308)))
+(assert (equal? -1.0 (numerator -5.5626846462680034577256e-309)))
+(assert (equal? +inf.0 (denominator -5.5626846462680034577256e-309)))
+(assert (equal? 1.0 (numerator 4.9406564584124654417657e-324)))
+(assert (equal? +inf.0 (denominator 4.9406564584124654417657e-324)))
+
 (assert (equal? 4.0 (denominator -1.25)))
 (assert (equal? 1e10 (numerator 1e10)))
 (assert (equal? 1.0 (denominator 1e10)))