real == frac

real == frac

[Kevin Ryde]

I noticed that

	(= 0.5 (+ 1/2 (/ 1 (ash 1 1000))))

gives #t, but I think with exact fractions it should now be #f, since
the second arg is 1/2+1/2^1000, not an actual 0.5.

I guess fraction2double is rounding, making the comparison not quite
right.  Going the other way, converting double to fraction and
comparing fractions would avoid that.

I suspect the same applies to all uses of fraction2double in <, min
and max too.


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: real == frac

[Marius Vollmer]

Kevin Ryde <user42@zip.com.au> writes:

> I noticed that
>
> 	(= 0.5 (+ 1/2 (/ 1 (ash 1 1000))))
>
> gives #t, but I think with exact fractions it should now be #f, since
> the second arg is 1/2+1/2^1000, not an actual 0.5.
>
> I guess fraction2double is rounding, making the comparison not quite
> right.  Going the other way, converting double to fraction and
> comparing fractions would avoid that.

Yes, that would be an improvement.  Could you implement it?

-- 
GPG: D5D4E405 - 2F9B BCCC 8527 692A 04E3  331E FAF8 226A D5D4 E405


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: real == frac

[Kevin Ryde]

[-- Attachment #1: Type: text/plain, Size: 452 bytes --]

Marius Vollmer <mvo@zagadka.de> writes:
>
> Yes, that would be an improvement.  Could you implement it?

Starting with less_p,

        * numbers.c (scm_less_p): Don't convert frac to float for compares,
        that can give wrong results through rounding.

        * tests/numbers.test (<): Add tests inum/bignum/flonum/frac with frac.

For min and max, I'd be inclined to have them call scm_less_p for
their comparison, to avoid duplicating code.



[-- Attachment #2: numbers.c.less-frac.diff --]
[-- Type: text/plain, Size: 3725 bytes --]

--- numbers.c.~1.219.~	2003-12-03 07:37:10.000000000 +1000
+++ numbers.c	2003-12-09 14:50:46.000000000 +1000
@@ -3074,6 +3074,12 @@
 }
 
 
+/* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
+   done are good for inums, but for bignums an answer can almost always be
+   had by just examining a few high bits of the operands, as done in GMP by
+   mpq_cmp.  flonum/frac compares likewise, but with the slight complication
+   of the float exponent to take into account.  */
+
 SCM_GPROC1 (s_less_p, "<", scm_tc7_rpsubr, scm_less_p, g_less_p);
 /* "Return @code{#t} if the list of parameters is monotonically\n"
  * "increasing."
@@ -3081,6 +3087,7 @@
 SCM
 scm_less_p (SCM x, SCM y)
 {
+ again:
   if (SCM_INUMP (x))
     {
       long xx = SCM_INUM (x);
@@ -3098,7 +3105,13 @@
       else if (SCM_REALP (y))
 	return SCM_BOOL ((double) xx < SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
-	return SCM_BOOL ((double) xx < scm_i_fraction2double (y));
+        {
+          /* "x < a/b" becomes "x*b < a" */
+        int_frac:
+          x = scm_product (x, SCM_FRACTION_DENOMINATOR (y));
+          y = SCM_FRACTION_NUMERATOR (y);
+          goto again;
+        }
       else
 	SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
     }
@@ -3126,12 +3139,7 @@
 	  return SCM_BOOL (cmp < 0);
 	}
       else if (SCM_FRACTIONP (y))
-	{
-	  int cmp;
-	  cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), scm_i_fraction2double (y));
-	  scm_remember_upto_here_1 (x);
-	  return SCM_BOOL (cmp < 0);
-	}
+        goto int_frac;
       else
 	SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
     }
@@ -3151,25 +3159,48 @@
       else if (SCM_REALP (y))
 	return SCM_BOOL (SCM_REAL_VALUE (x) < SCM_REAL_VALUE (y));
       else if (SCM_FRACTIONP (y))
-	return SCM_BOOL (SCM_REAL_VALUE (x) < scm_i_fraction2double (y));
+        {
+          double  xx = SCM_REAL_VALUE (x);
+	  if (xisnan (xx))
+	    return SCM_BOOL_F;
+          if (xisinf (xx))
+            return SCM_BOOL (xx < 0.0);
+          x = scm_inexact_to_exact (x);  /* with x as frac or int */
+          goto again;
+        }
       else
 	SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
     }
   else if (SCM_FRACTIONP (x))
     {
-      if (SCM_INUMP (y))
-	return SCM_BOOL (scm_i_fraction2double (x) < (double) SCM_INUM (y));
-      else if (SCM_BIGP (y))
-	{
-	  int cmp;
-	  cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), scm_i_fraction2double (x));
-	  scm_remember_upto_here_1 (y);
-	  return SCM_BOOL (cmp > 0);
-	}
+      if (SCM_INUMP (y) || SCM_BIGP (y))
+        {
+          /* "a/b < y" becomes "a < y*b" */
+          y = scm_product (y, SCM_FRACTION_DENOMINATOR (x));
+          x = SCM_FRACTION_NUMERATOR (x);
+          goto again;
+        }
       else if (SCM_REALP (y))
-	return SCM_BOOL (scm_i_fraction2double (x) < SCM_REAL_VALUE (y));
+        {
+          double yy = SCM_REAL_VALUE (y);
+          if (xisnan (yy))
+            return SCM_BOOL_F;
+          if (xisinf (yy))
+            return SCM_BOOL (0.0 < yy);
+          y = scm_inexact_to_exact (y);  /* with y as frac or int */
+          goto again;
+        }
       else if (SCM_FRACTIONP (y))
-	return SCM_BOOL (scm_i_fraction2double (x) < scm_i_fraction2double (y));
+        {
+          /* "a/b < c/d" becomes "a*d < c*b" */
+          SCM new_x = scm_product (SCM_FRACTION_NUMERATOR (x),
+                                   SCM_FRACTION_DENOMINATOR (y));
+          SCM new_y = scm_product (SCM_FRACTION_NUMERATOR (y),
+                                   SCM_FRACTION_DENOMINATOR (x));
+          x = new_x;
+          y = new_y;
+          goto again;
+        }
       else
 	SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
     }

[-- Attachment #3: numbers.test.less-frac.diff --]
[-- Type: text/plain, Size: 3078 bytes --]

--- numbers.test.~1.39.~	2003-11-25 08:11:18.000000000 +1000
+++ numbers.test	2003-12-09 15:06:40.000000000 +1000
@@ -1684,7 +1684,95 @@
   (pass-if (not (< (1- (ash 3 1023)) +nan.0)))
   (pass-if (not (< +nan.0 (ash 3 1023))))
   (pass-if (not (< +nan.0 (1+ (ash 3 1023)))))
-  (pass-if (not (< +nan.0 (1- (ash 3 1023))))))
+  (pass-if (not (< +nan.0 (1- (ash 3 1023)))))
+
+  (with-test-prefix "inum/frac"
+    (pass-if (< 2 9/4))
+    (pass-if (< -2 9/4))
+    (pass-if (< -2 7/4))
+    (pass-if (< -2 -7/4))
+    (pass-if (eq? #f (< 2 7/4)))
+    (pass-if (eq? #f (< 2 -7/4)))
+    (pass-if (eq? #f (< 2 -9/4)))
+    (pass-if (eq? #f (< -2 -9/4))))
+
+  (with-test-prefix "bignum/frac"
+    (let ((x (ash 1 2048)))
+      (pass-if (< x (* 4/3 x)))
+      (pass-if (< (- x) (* 4/3 x)))
+      (pass-if (< (- x) (* 2/3 x)))
+      (pass-if (< (- x) (* -2/3 x)))
+      (pass-if (eq? #f (< x (* 2/3 x))))
+      (pass-if (eq? #f (< x (* -2/3 x))))
+      (pass-if (eq? #f (< x (* -4/3 x))))
+      (pass-if (eq? #f (< (- x) (* -4/3 x))))))
+
+  (with-test-prefix "flonum/frac"
+    (pass-if (< 0.75 4/3))
+    (pass-if (< -0.75 4/3))
+    (pass-if (< -0.75 2/3))
+    (pass-if (< -0.75 -2/3))
+    (pass-if (eq? #f (< 0.75 2/3)))
+    (pass-if (eq? #f (< 0.75 -2/3)))
+    (pass-if (eq? #f (< 0.75 -4/3)))
+    (pass-if (eq? #f (< -0.75 -4/3)))
+
+    (pass-if (< -inf.0 4/3))
+    (pass-if (< -inf.0 -4/3))
+    (pass-if (eq? #f (< +inf.0 4/3)))
+    (pass-if (eq? #f (< +inf.0 -4/3)))
+
+    (pass-if (eq? #f (< +nan.0 4/3)))
+    (pass-if (eq? #f (< +nan.0 -4/3))))
+
+  (with-test-prefix "frac/inum"
+    (pass-if (< 7/4 2))
+    (pass-if (< -7/4 2))
+    (pass-if (< -9/4 2))
+    (pass-if (< -9/4 -2))
+    (pass-if (eq? #f (< 9/4 2)))
+    (pass-if (eq? #f (< 9/4 -2)))
+    (pass-if (eq? #f (< 7/4 -2)))
+    (pass-if (eq? #f (< -7/4 -2))))
+
+  (with-test-prefix "frac/bignum"
+    (let ((x (ash 1 2048)))
+      (pass-if (< (* 2/3 x) x))
+      (pass-if (< (* -2/3 x) x))
+      (pass-if (< (* -4/3 x) x))
+      (pass-if (< (* -4/3 x) (- x)))
+      (pass-if (eq? #f (< (* 4/3 x) x)))
+      (pass-if (eq? #f (< (* 4/3 x) (- x))))
+      (pass-if (eq? #f (< (* 2/3 x) (- x))))
+      (pass-if (eq? #f (< (* -2/3 x) (- x))))))
+
+  (with-test-prefix "frac/flonum"
+    (pass-if (< 2/3 0.75))
+    (pass-if (< -2/3 0.75))
+    (pass-if (< -4/3 0.75))
+    (pass-if (< -4/3 -0.75))
+    (pass-if (eq? #f (< 4/3 0.75)))
+    (pass-if (eq? #f (< 4/3 -0.75)))
+    (pass-if (eq? #f (< 2/3 -0.75)))
+    (pass-if (eq? #f (< -2/3 -0.75)))
+
+    (pass-if (< 4/3 +inf.0))
+    (pass-if (< -4/3 +inf.0))
+    (pass-if (eq? #f (< 4/3 -inf.0)))
+    (pass-if (eq? #f (< -4/3 -inf.0)))
+
+    (pass-if (eq? #f (< 4/3 +nan.0)))
+    (pass-if (eq? #f (< -4/3 +nan.0))))
+
+  (with-test-prefix "frac/frac"
+    (pass-if (< 2/3 6/7))
+    (pass-if (< -2/3 6/7))
+    (pass-if (< -4/3 6/7))
+    (pass-if (< -4/3 -6/7))
+    (pass-if (eq? #f (< 4/3 6/7)))
+    (pass-if (eq? #f (< 4/3 -6/7)))
+    (pass-if (eq? #f (< 2/3 -6/7)))
+    (pass-if (eq? #f (< -2/3 -6/7)))))
 
 ;;;
 ;;; >

[-- Attachment #4: Type: text/plain, Size: 142 bytes --]

_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: real == frac

[Kevin Ryde]

[-- Attachment #1: Type: text/plain, Size: 392 bytes --]

        * numbers.c (scm_num_eq_p): For real==frac, complex==frac, frac==real
        and frac==complex, make an exact comparison rather than converting
        with fraction2double.

I've added code under the complex cases because it was there already.
But I might have mentioned before that I thought all those could be
false immediately, if COMPLEXP always has a non-zero imaginary part.


[-- Attachment #2: numbers.c.eq-frac.diff --]
[-- Type: text/plain, Size: 2650 bytes --]

--- numbers.c.~1.221.~	2004-01-07 07:54:59.000000000 +1000
+++ numbers.c	2004-02-12 11:58:18.000000000 +1000
@@ -2945,6 +2945,7 @@
 SCM
 scm_num_eq_p (SCM x, SCM y)
 {
+ again:
   if (SCM_INUMP (x))
     {
       long xx = SCM_INUM (x);
@@ -3019,7 +3020,15 @@
 	return SCM_BOOL ((SCM_REAL_VALUE (x) == SCM_COMPLEX_REAL (y))
 			 && (0.0 == SCM_COMPLEX_IMAG (y)));
       else if (SCM_FRACTIONP (y))
-	return SCM_BOOL (SCM_REAL_VALUE (x) == scm_i_fraction2double (y));
+        {
+          double  xx = SCM_REAL_VALUE (x);
+          if (xisnan (xx))
+            return SCM_BOOL_F;
+          if (xisinf (xx))
+            return SCM_BOOL (xx < 0.0);
+          x = scm_inexact_to_exact (x);  /* with x as frac or int */
+          goto again;
+        }
       else
 	SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
     }
@@ -3046,8 +3055,18 @@
 	return SCM_BOOL ((SCM_COMPLEX_REAL (x) == SCM_COMPLEX_REAL (y))
 			 && (SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y)));
       else if (SCM_FRACTIONP (y))
-	return SCM_BOOL ((SCM_COMPLEX_REAL (x) == scm_i_fraction2double (y))
-			 && (SCM_COMPLEX_IMAG (x) == 0.0));
+        {
+          double  xx;
+          if (SCM_COMPLEX_IMAG (x) != 0.0)
+            return SCM_BOOL_F;
+          xx = SCM_COMPLEX_REAL (x);
+          if (xisnan (xx))
+            return SCM_BOOL_F;
+          if (xisinf (xx))
+            return SCM_BOOL (xx < 0.0);
+          x = scm_inexact_to_exact (x);  /* with x as frac or int */
+          goto again;
+        }
       else
 	SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
     }
@@ -3058,10 +3077,28 @@
       else if (SCM_BIGP (y))
 	return SCM_BOOL_F;
       else if (SCM_REALP (y))
-	return SCM_BOOL (scm_i_fraction2double (x) == SCM_REAL_VALUE (y));
+        {
+          double yy = SCM_REAL_VALUE (y);
+          if (xisnan (yy))
+            return SCM_BOOL_F;
+          if (xisinf (yy))
+            return SCM_BOOL (0.0 < yy);
+          y = scm_inexact_to_exact (y);  /* with y as frac or int */
+          goto again;
+        }
       else if (SCM_COMPLEXP (y))
-	return SCM_BOOL ((scm_i_fraction2double (x) == SCM_COMPLEX_REAL (y))
-			 && (0.0 == SCM_COMPLEX_IMAG (y)));
+        {
+          double yy;
+          if (SCM_COMPLEX_IMAG (y) != 0.0)
+            return SCM_BOOL_F;
+          yy = SCM_COMPLEX_REAL (y);
+          if (xisnan (yy))
+            return SCM_BOOL_F;
+          if (xisinf (yy))
+            return SCM_BOOL (0.0 < yy);
+          y = scm_inexact_to_exact (y);  /* with y as frac or int */
+          goto again;
+        }
       else if (SCM_FRACTIONP (y))
 	return scm_i_fraction_equalp (x, y);
       else

[-- Attachment #3: numbers.test.eq-frac.diff --]
[-- Type: text/plain, Size: 935 bytes --]

--- numbers.test.~1.40.~	2003-12-09 15:06:40.000000000 +1000
+++ numbers.test	2004-02-17 09:07:42.000000000 +1000
@@ -1308,7 +1308,25 @@
   ;; in gmp prior to 4.2, mpz_cmp_d ended up treating NaN as 3*2^1023, make
   ;; sure we've avoided that
   (pass-if (not (= (ash 3 1023) +nan.0)))
-  (pass-if (not (= +nan.0 (ash 3 1023)))))
+  (pass-if (not (= +nan.0 (ash 3 1023))))
+
+  (pass-if (= 1/2 0.5))
+  (pass-if (not (= 1/3 0.333333333333333333333333333333333)))
+  (pass-if (not (= 2/3 0.5)))
+  (pass-if (not (= 0.5 (+ 1/2 (/ 1 (ash 1 1000))))))
+
+  (pass-if (= 1/2 0.5+0i))
+  (pass-if (not (= 0.333333333333333333333333333333333 1/3)))
+  (pass-if (not (= 2/3 0.5+0i)))
+  (pass-if (not (= 1/2 0+0.5i)))
+
+  (pass-if (= 0.5 1/2))
+  (pass-if (not (= 0.5 2/3)))
+  (pass-if (not (= (+ 1/2 (/ 1 (ash 1 1000))) 0.5)))
+
+  (pass-if (= 0.5+0i 1/2))
+  (pass-if (not (= 0.5+0i 2/3)))
+  (pass-if (not (= 0+0.5i 1/2))))
 
 ;;;
 ;;; <

[-- Attachment #4: Type: text/plain, Size: 142 bytes --]

_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel