[Guile-commits] 43/85: Simplify scm_product, use integer lib

[Andy Wingo]

wingo pushed a commit to branch main
in repository guile.

commit 9179525a058affd702971a0dd57926acabd5d0cb
Author: Andy Wingo <wingo@pobox.com>
AuthorDate: Tue Jan 4 20:46:50 2022 +0100

    Simplify scm_product, use integer lib
    
    * libguile/numbers.c (scm_product): Remove need for s_product defines.
    Call out to product, as appropriate.
    (product): New helper.
    * libguile/integers.h:
    * libguile/integers.c (scm_integer_mul_ii):
    (scm_integer_mul_zi):
    (scm_integer_mul_zz): New internal functions.
---
 libguile/integers.c |  55 +++++++++++
 libguile/integers.h |   4 +
 libguile/numbers.c  | 277 +++++++++++++++++-----------------------------------
 3 files changed, 147 insertions(+), 189 deletions(-)

diff --git a/libguile/integers.c b/libguile/integers.c
index 1133d2215..715c210f5 100644
--- a/libguile/integers.c
+++ b/libguile/integers.c
@@ -2603,3 +2603,58 @@ scm_integer_sub_zz (struct scm_bignum *x, struct 
scm_bignum *y)
 {
   return scm_integer_add_zz (x, negate_bignum (clone_bignum (y)));
 }
+
+SCM
+scm_integer_mul_ii (scm_t_inum x, scm_t_inum y)
+{
+#if SCM_I_FIXNUM_BIT < 32
+  int64_t k = x * (int64_t) y;
+  if (SCM_FIXABLE (k))
+    return SCM_I_MAKINUM (k);
+#else
+  if (x == 0)
+    return SCM_INUM0;
+  scm_t_inum ax = (x > 0) ? x : -x;
+  scm_t_inum ay = (y > 0) ? y : -y;
+  if (SCM_MOST_POSITIVE_FIXNUM / ax >= ay)
+    return SCM_I_MAKINUM (x * y);
+#endif
+
+  // FIXME: Use mpn_mul with two-limb result to avoid allocating.
+  return scm_integer_mul_zi (long_to_bignum (x), y);
+}
+
+SCM
+scm_integer_mul_zi (struct scm_bignum *x, scm_t_inum y)
+{
+  switch (y)
+    {
+    case -1:
+      return scm_integer_negate_z (x);
+    case 0:
+      return SCM_INUM0;
+    case 1:
+      return scm_from_bignum (x);
+    default:
+      {
+        mpz_t result, zx;
+        mpz_init (result);
+        alias_bignum_to_mpz (x, zx);
+        mpz_mul_si (result, zx, y);
+        scm_remember_upto_here_1 (x);
+        return take_mpz (result);
+      }
+    }
+}
+
+SCM
+scm_integer_mul_zz (struct scm_bignum *x, struct scm_bignum *y)
+{
+  mpz_t result, zx, zy;
+  mpz_init (result);
+  alias_bignum_to_mpz (x, zx);
+  alias_bignum_to_mpz (y, zy);
+  mpz_mul (result, zx, zy);
+  scm_remember_upto_here_2 (x, y);
+  return take_mpz (result);
+}
diff --git a/libguile/integers.h b/libguile/integers.h
index 3f3bf279e..8795a6288 100644
--- a/libguile/integers.h
+++ b/libguile/integers.h
@@ -180,6 +180,10 @@ SCM_INTERNAL SCM scm_integer_sub_iz (scm_t_inum x, struct 
scm_bignum *y);
 SCM_INTERNAL SCM scm_integer_sub_zi (struct scm_bignum *x, scm_t_inum y);
 SCM_INTERNAL SCM scm_integer_sub_zz (struct scm_bignum *x, struct scm_bignum 
*y);
 
+SCM_INTERNAL SCM scm_integer_mul_ii (scm_t_inum x, scm_t_inum y);
+SCM_INTERNAL SCM scm_integer_mul_zi (struct scm_bignum *x, scm_t_inum y);
+SCM_INTERNAL SCM scm_integer_mul_zz (struct scm_bignum *x, struct scm_bignum 
*y);
+
 
 
 #endif  /* SCM_INTEGERS_H */
diff --git a/libguile/numbers.c b/libguile/numbers.c
index 80c775f28..a2f9de14e 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5417,237 +5417,136 @@ SCM_DEFINE (scm_oneminus, "1-", 1, 0, 0,
 }
 #undef FUNC_NAME
 
-
-SCM_PRIMITIVE_GENERIC (scm_i_product, "*", 0, 2, 1,
-                       (SCM x, SCM y, SCM rest),
-                       "Return the product of all arguments.  If called 
without arguments,\n"
-                       "1 is returned.")
-#define FUNC_NAME s_scm_i_product
-{
-  while (!scm_is_null (rest))
-    { x = scm_product (x, y);
-      y = scm_car (rest);
-      rest = scm_cdr (rest);
-    }
-  return scm_product (x, y);
-}
-#undef FUNC_NAME
-                       
-#define s_product s_scm_i_product
-#define g_product g_scm_i_product
-
-SCM
-scm_product (SCM x, SCM y)
+static SCM
+product (SCM x, SCM y)
 {
-  if (SCM_UNLIKELY (SCM_UNBNDP (y)))
-    {
-      if (SCM_UNBNDP (x))
-       return SCM_I_MAKINUM (1L);
-      else if (SCM_NUMBERP (x))
-       return x;
-      else
-       return scm_wta_dispatch_1 (g_product, x, SCM_ARG1, s_product);
-    }
-  
-  if (SCM_LIKELY (SCM_I_INUMP (x)))
+  if (SCM_I_INUMP (x))
     {
-      scm_t_inum xx;
-
-    xinum:
-      xx = SCM_I_INUM (x);
-
-      switch (xx)
-       {
-        case 1:
-         /* exact1 is the universal multiplicative identity */
-         return y;
-         break;
-        case 0:
-         /* exact0 times a fixnum is exact0: optimize this case */
-         if (SCM_LIKELY (SCM_I_INUMP (y)))
-           return SCM_INUM0;
-         /* if the other argument is inexact, the result is inexact,
-            and we must do the multiplication in order to handle
-            infinities and NaNs properly. */
-         else if (SCM_REALP (y))
-           return scm_i_from_double (0.0 * SCM_REAL_VALUE (y));
-         else if (SCM_COMPLEXP (y))
-           return scm_c_make_rectangular (0.0 * SCM_COMPLEX_REAL (y),
-                                          0.0 * SCM_COMPLEX_IMAG (y));
-         /* we've already handled inexact numbers,
-            so y must be exact, and we return exact0 */
-         else if (SCM_NUMP (y))
-           return SCM_INUM0;
-         else
-           return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
-         break;
-       }
-
-      if (SCM_LIKELY (SCM_I_INUMP (y)))
-       {
-         scm_t_inum yy = SCM_I_INUM (y);
-#if SCM_I_FIXNUM_BIT < 32 && SCM_HAVE_T_INT64
-          int64_t kk = xx * (int64_t) yy;
-          if (SCM_FIXABLE (kk))
-            return SCM_I_MAKINUM (kk);
-#else
-          scm_t_inum axx = (xx > 0) ? xx : -xx;
-          scm_t_inum ayy = (yy > 0) ? yy : -yy;
-          if (SCM_MOST_POSITIVE_FIXNUM / axx >= ayy)
-            return SCM_I_MAKINUM (xx * yy);
-#endif
-         else
-           {
-             SCM result = scm_i_inum2big (xx);
-             mpz_mul_si (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), yy);
-             return scm_i_normbig (result);
-           }
-       }
+      if (scm_is_eq (x, SCM_I_MAKINUM (-1)))
+        return negate (y);
+      else if (SCM_I_INUMP (y))
+        return scm_integer_mul_ii (SCM_I_INUM (x), SCM_I_INUM (y));
       else if (SCM_BIGP (y))
-       {
-          /* There is one bignum which, when multiplied by negative one,
-             becomes a non-zero fixnum: (1+ most-positive-fixum).  Since
-             we know the type of X and Y are numbers, delegate this
-             special case to scm_difference.  */
-          if (xx == -1)
-            return scm_difference (y, SCM_UNDEFINED);
-          else
-            {
-              SCM result = scm_i_mkbig ();
-              mpz_mul_si (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (y), xx);
-              scm_remember_upto_here_1 (y);
-              return result;
-            }
-       }
+        return scm_integer_mul_zi (scm_bignum (y), SCM_I_INUM (x));
       else if (SCM_REALP (y))
-       return scm_i_from_double (xx * SCM_REAL_VALUE (y));
+        return scm_i_from_double (SCM_I_INUM (x) * SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
-       return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
-                                xx * SCM_COMPLEX_IMAG (y));
+        return scm_c_make_rectangular (SCM_I_INUM (x) * SCM_COMPLEX_REAL (y),
+                                       SCM_I_INUM (x) * SCM_COMPLEX_IMAG (y));
       else if (SCM_FRACTIONP (y))
-       return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
-                              SCM_FRACTION_DENOMINATOR (y));
-      else
-       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
+        return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
+                                 SCM_FRACTION_DENOMINATOR (y));
+      abort (); /* Unreachable.  */
     }
   else if (SCM_BIGP (x))
     {
-      if (SCM_I_INUMP (y))
-       {
-         SCM_SWAP (x, y);
-         goto xinum;
-       }
-      else if (SCM_BIGP (y))
-       {
-         SCM result = scm_i_mkbig ();
-         mpz_mul (SCM_I_BIG_MPZ (result),
-                  SCM_I_BIG_MPZ (x),
-                  SCM_I_BIG_MPZ (y));
-         scm_remember_upto_here_2 (x, y);
-         return result;
-       }
+      if (SCM_BIGP (y))
+        return scm_integer_mul_zz (scm_bignum (x), scm_bignum (y));
       else if (SCM_REALP (y))
-       {
-         double result = mpz_get_d (SCM_I_BIG_MPZ (x)) * SCM_REAL_VALUE (y);
-         scm_remember_upto_here_1 (x);
-         return scm_i_from_double (result);
-       }
+        return scm_from_double (scm_integer_to_double_z (scm_bignum (x))
+                                * SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
        {
-         double z = mpz_get_d (SCM_I_BIG_MPZ (x));
-         scm_remember_upto_here_1 (x);
+         double z = scm_integer_to_double_z (scm_bignum (x));
          return scm_c_make_rectangular (z * SCM_COMPLEX_REAL (y),
-                                  z * SCM_COMPLEX_IMAG (y));
+                                         z * SCM_COMPLEX_IMAG (y));
        }
       else if (SCM_FRACTIONP (y))
-       return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
-                              SCM_FRACTION_DENOMINATOR (y));
+       return scm_i_make_ratio (product (x, SCM_FRACTION_NUMERATOR (y)),
+                                 SCM_FRACTION_DENOMINATOR (y));
       else
-       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
+        return product (y, x);
     }
   else if (SCM_REALP (x))
     {
-      if (SCM_I_INUMP (y))
-       {
-         SCM_SWAP (x, y);
-         goto xinum;
-       }
-      else if (SCM_BIGP (y))
-       {
-         double result = mpz_get_d (SCM_I_BIG_MPZ (y)) * SCM_REAL_VALUE (x);
-         scm_remember_upto_here_1 (y);
-         return scm_i_from_double (result);
-       }
-      else if (SCM_REALP (y))
+      if (SCM_REALP (y))
        return scm_i_from_double (SCM_REAL_VALUE (x) * SCM_REAL_VALUE (y));
       else if (SCM_COMPLEXP (y))
-       return scm_c_make_rectangular (SCM_REAL_VALUE (x) * SCM_COMPLEX_REAL 
(y),
-                                SCM_REAL_VALUE (x) * SCM_COMPLEX_IMAG (y));
+       return scm_c_make_rectangular
+          (SCM_REAL_VALUE (x) * SCM_COMPLEX_REAL (y),
+           SCM_REAL_VALUE (x) * SCM_COMPLEX_IMAG (y));
       else if (SCM_FRACTIONP (y))
-       return scm_i_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double 
(y));
+       return scm_i_from_double
+          (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
       else
-       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
+        return product (y, x);
     }
   else if (SCM_COMPLEXP (x))
     {
-      if (SCM_I_INUMP (y))
-       {
-         SCM_SWAP (x, y);
-         goto xinum;
-       }
-      else if (SCM_BIGP (y))
-       {
-         double z = mpz_get_d (SCM_I_BIG_MPZ (y));
-         scm_remember_upto_here_1 (y);
-         return scm_c_make_rectangular (z * SCM_COMPLEX_REAL (x),
-                                  z * SCM_COMPLEX_IMAG (x));
-       }
-      else if (SCM_REALP (y))
-       return scm_c_make_rectangular (SCM_REAL_VALUE (y) * SCM_COMPLEX_REAL 
(x),
-                                SCM_REAL_VALUE (y) * SCM_COMPLEX_IMAG (x));
-      else if (SCM_COMPLEXP (y))
+      if (SCM_COMPLEXP (y))
        {
-         return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) * 
SCM_COMPLEX_REAL (y)
-                                  - SCM_COMPLEX_IMAG (x) * SCM_COMPLEX_IMAG 
(y),
-                                  SCM_COMPLEX_REAL (x) * SCM_COMPLEX_IMAG (y)
-                                  + SCM_COMPLEX_IMAG (x) * SCM_COMPLEX_REAL 
(y));
+          double rx = SCM_COMPLEX_REAL (x), ry = SCM_COMPLEX_REAL (y);
+          double ix = SCM_COMPLEX_IMAG (x), iy = SCM_COMPLEX_IMAG (y);
+         return scm_c_make_rectangular (rx * ry - ix * iy, rx * iy + ix * ry);
        }
       else if (SCM_FRACTIONP (y))
        {
          double yy = scm_i_fraction2double (y);
          return scm_c_make_rectangular (yy * SCM_COMPLEX_REAL (x),
-                                  yy * SCM_COMPLEX_IMAG (x));
+                                         yy * SCM_COMPLEX_IMAG (x));
        }
       else
-       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
+        return product (y, x);
     }
   else if (SCM_FRACTIONP (x))
     {
-      if (SCM_I_INUMP (y))
-       return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
-                              SCM_FRACTION_DENOMINATOR (x));
-      else if (SCM_BIGP (y))
-       return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
-                              SCM_FRACTION_DENOMINATOR (x));
-      else if (SCM_REALP (y))
-       return scm_i_from_double (scm_i_fraction2double (x) * SCM_REAL_VALUE 
(y));
-      else if (SCM_COMPLEXP (y))
-       {
-         double xx = scm_i_fraction2double (x);
-         return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
-                                  xx * SCM_COMPLEX_IMAG (y));
-       }
-      else if (SCM_FRACTIONP (y))
+      if (SCM_FRACTIONP (y))
        /* a/b * c/d = ac / bd */
        return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x),
-                                           SCM_FRACTION_NUMERATOR (y)),
-                              scm_product (SCM_FRACTION_DENOMINATOR (x),
-                                           SCM_FRACTION_DENOMINATOR (y)));
+                                              SCM_FRACTION_NUMERATOR (y)),
+                                 scm_product (SCM_FRACTION_DENOMINATOR (x),
+                                              SCM_FRACTION_DENOMINATOR (y)));
       else
-       return scm_wta_dispatch_2 (g_product, x, y, SCM_ARGn, s_product);
+       return product (y, x);
     }
   else
-    return scm_wta_dispatch_2 (g_product, x, y, SCM_ARG1, s_product);
+    abort (); /* Unreachable.  */
+}
+
+SCM_PRIMITIVE_GENERIC (scm_i_product, "*", 0, 2, 1,
+                       (SCM x, SCM y, SCM rest),
+                       "Return the product of all arguments.  If called 
without arguments,\n"
+                       "1 is returned.")
+#define FUNC_NAME s_scm_i_product
+{
+  while (!scm_is_null (rest))
+    { x = scm_product (x, y);
+      y = scm_car (rest);
+      rest = scm_cdr (rest);
+    }
+  return scm_product (x, y);
+}
+#undef FUNC_NAME
+
+SCM
+scm_product (SCM x, SCM y)
+{
+  if (SCM_UNBNDP (y))
+    {
+      if (SCM_UNBNDP (x))
+       return SCM_I_MAKINUM (1L);
+      else if (SCM_NUMBERP (x))
+       return x;
+      else
+       return scm_wta_dispatch_1 (g_scm_i_product, x, SCM_ARG1,
+                                   s_scm_i_product);
+    }
+
+  /* This is pretty gross!  But (* 1 X) is apparently X in Guile, for
+     any type of X, even a pair.  */
+  if (scm_is_eq (x, SCM_INUM1))
+    return y;
+  if (scm_is_eq (y, SCM_INUM1))
+    return x;
+
+  if (!SCM_NUMBERP (x))
+    return scm_wta_dispatch_2 (g_scm_i_product, x, y, SCM_ARG1,
+                               s_scm_i_product);
+  if (!SCM_NUMBERP (y))
+    return scm_wta_dispatch_2 (g_scm_i_product, x, y, SCM_ARG2,
+                               s_scm_i_product);
+
+  return product (x, y);
 }
 
 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \