merge from 1.8 branch
[bpt/guile.git] / libguile / numbers.c
index d33edde..3b6d781 100644 (file)
@@ -1,4 +1,4 @@
-/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 Free Software Foundation, Inc.
+/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006 Free Software Foundation, Inc.
  *
  * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
  * and Bellcore.  See scm_divide.
@@ -16,7 +16,7 @@
  *
  * You should have received a copy of the GNU Lesser General Public
  * License along with this library; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  */
 
 \f
@@ -183,7 +183,7 @@ scm_i_ulong2big (unsigned long x)
   return z;
 }
 
-SCM_C_INLINE_KEYWORD static SCM
+SCM_C_INLINE_KEYWORD SCM
 scm_i_clonebig (SCM src_big, int same_sign_p)
 {
   /* Copy src_big's value, negate it if same_sign_p is false, and return. */
@@ -1309,7 +1309,7 @@ SCM_DEFINE1 (scm_logior, "logior", scm_tc7_asubr,
            mpz_ior (SCM_I_BIG_MPZ (result_z), nn1_z, SCM_I_BIG_MPZ (n2));
            scm_remember_upto_here_1 (n2);
            mpz_clear (nn1_z);
-           return result_z;
+           return scm_i_normbig (result_z);
          }
        }
       else
@@ -1330,7 +1330,7 @@ SCM_DEFINE1 (scm_logior, "logior", scm_tc7_asubr,
                   SCM_I_BIG_MPZ (n1),
                   SCM_I_BIG_MPZ (n2));
          scm_remember_upto_here_2 (n1, n2);
-         return result_z;
+         return scm_i_normbig (result_z);
        }
       else
        SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
@@ -2937,7 +2937,8 @@ mem2complex (const char* mem, size_t len, unsigned int idx,
 enum t_radix {NO_RADIX=0, DUAL=2, OCT=8, DEC=10, HEX=16};
 
 SCM
-scm_i_mem2number (const char* mem, size_t len, unsigned int default_radix)
+scm_c_locale_stringn_to_number (const char* mem, size_t len,
+                               unsigned int default_radix)
 {
   unsigned int idx = 0;
   unsigned int radix = NO_RADIX;
@@ -3043,9 +3044,9 @@ SCM_DEFINE (scm_string_to_number, "string->number", 1, 1, 0,
   else
     base = scm_to_unsigned_integer (radix, 2, INT_MAX);
 
-  answer = scm_i_mem2number (scm_i_string_chars (string),
-                            scm_i_string_length (string),
-                            base);
+  answer = scm_c_locale_stringn_to_number (scm_i_string_chars (string),
+                                          scm_i_string_length (string),
+                                          base);
   scm_remember_upto_here_1 (string);
   return answer;
 }
@@ -3211,7 +3212,27 @@ scm_num_eq_p (SCM x, SCM y)
       else if (SCM_BIGP (y))
        return SCM_BOOL_F;
       else if (SCM_REALP (y))
-       return scm_from_bool ((double) xx == SCM_REAL_VALUE (y));
+        {
+          /* On a 32-bit system an inum fits a double, we can cast the inum
+             to a double and compare.
+
+             But on a 64-bit system an inum is bigger than a double and
+             casting it to a double (call that dxx) will round.  dxx is at
+             worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
+             an integer and fits a long.  So we cast yy to a long and
+             compare with plain xx.
+
+             An alternative (for any size system actually) would be to check
+             yy is an integer (with floor) and is in range of an inum
+             (compare against appropriate powers of 2) then test
+             xx==(long)yy.  It's just a matter of which casts/comparisons
+             might be fastest or easiest for the cpu.  */
+
+          double yy = SCM_REAL_VALUE (y);
+          return scm_from_bool ((double) xx == yy
+                               && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
+                                   || xx == (long) yy));
+        }
       else if (SCM_COMPLEXP (y))
        return scm_from_bool (((double) xx == SCM_COMPLEX_REAL (y))
                         && (0.0 == SCM_COMPLEX_IMAG (y)));
@@ -3257,8 +3278,15 @@ scm_num_eq_p (SCM x, SCM y)
     }
   else if (SCM_REALP (x))
     {
+      double xx = SCM_REAL_VALUE (x);
       if (SCM_I_INUMP (y))
-       return scm_from_bool (SCM_REAL_VALUE (x) == (double) SCM_I_INUM (y));
+        {
+          /* see comments with inum/real above */
+          long yy = SCM_I_INUM (y);
+          return scm_from_bool (xx == (double) yy
+                               && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
+                                   || (long) xx == yy));
+        }
       else if (SCM_BIGP (y))
        {
          int cmp;
@@ -4080,6 +4108,16 @@ scm_sum (SCM x, SCM y)
 }
 
 
+SCM_DEFINE (scm_oneplus, "1+", 1, 0, 0, 
+            (SCM x),
+           "Return @math{@var{x}+1}.")
+#define FUNC_NAME s_scm_oneplus
+{
+  return scm_sum (x, SCM_I_MAKINUM (1));
+}
+#undef FUNC_NAME
+
+
 SCM_GPROC1 (s_difference, "-", scm_tc7_asubr, scm_difference, g_difference);
 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
  * the sum of all but the first argument are subtracted from the first
@@ -4315,6 +4353,16 @@ scm_difference (SCM x, SCM y)
 #undef FUNC_NAME
 
 
+SCM_DEFINE (scm_oneminus, "1-", 1, 0, 0, 
+            (SCM x),
+           "Return @math{@var{x}-1}.")
+#define FUNC_NAME s_scm_oneminus
+{
+  return scm_difference (x, SCM_I_MAKINUM (1));
+}
+#undef FUNC_NAME
+
+
 SCM_GPROC1 (s_product, "*", scm_tc7_asubr, scm_product, g_product);
 /* "Return the product of all arguments.  If called without arguments,\n"
  * "1 is returned."
@@ -4575,7 +4623,7 @@ scm_i_divide (SCM x, SCM y, int inexact)
        {
          double r = SCM_COMPLEX_REAL (x);
          double i = SCM_COMPLEX_IMAG (x);
-         if (r <= i)
+         if (fabs(r) <= fabs(i))
            {
              double t = r / i;
              double d = i * (1.0 + t * t);
@@ -4647,7 +4695,7 @@ scm_i_divide (SCM x, SCM y, int inexact)
          {
            double r = SCM_COMPLEX_REAL (y);
            double i = SCM_COMPLEX_IMAG (y);
-           if (r <= i)
+           if (fabs(r) <= fabs(i))
              {
                double t = r / i;
                double d = i * (1.0 + t * t);
@@ -4731,28 +4779,33 @@ scm_i_divide (SCM x, SCM y, int inexact)
          else
            {
              /* big_x / big_y */
-             int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
-                                                SCM_I_BIG_MPZ (y));
-             if (divisible_p)
-               {
-                 SCM result = scm_i_mkbig ();
-                 mpz_divexact (SCM_I_BIG_MPZ (result),
-                               SCM_I_BIG_MPZ (x),
-                               SCM_I_BIG_MPZ (y));
-                 scm_remember_upto_here_2 (x, y);
-                 return scm_i_normbig (result);
-               }
-             else
-               {
-                 if (inexact)
-                   {
-                     double dbx = mpz_get_d (SCM_I_BIG_MPZ (x));
-                     double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
-                     scm_remember_upto_here_2 (x, y);
-                     return scm_from_double (dbx / dby);
-                   }
-                 else return scm_i_make_ratio (x, y);
-               }
+              if (inexact)
+                {
+                  /* It's easily possible for the ratio x/y to fit a double
+                     but one or both x and y be too big to fit a double,
+                     hence the use of mpq_get_d rather than converting and
+                     dividing.  */
+                  mpq_t q;
+                  *mpq_numref(q) = *SCM_I_BIG_MPZ (x);
+                  *mpq_denref(q) = *SCM_I_BIG_MPZ (y);
+                  return scm_from_double (mpq_get_d (q));
+                }
+              else
+                {
+                  int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
+                                                     SCM_I_BIG_MPZ (y));
+                  if (divisible_p)
+                    {
+                      SCM result = scm_i_mkbig ();
+                      mpz_divexact (SCM_I_BIG_MPZ (result),
+                                    SCM_I_BIG_MPZ (x),
+                                    SCM_I_BIG_MPZ (y));
+                      scm_remember_upto_here_2 (x, y);
+                      return scm_i_normbig (result);
+                    }
+                  else
+                    return scm_i_make_ratio (x, y);
+                }
            }
        }
       else if (SCM_REALP (y))
@@ -4852,7 +4905,7 @@ scm_i_divide (SCM x, SCM y, int inexact)
        {
          double ry = SCM_COMPLEX_REAL (y);
          double iy = SCM_COMPLEX_IMAG (y);
-         if (ry <= iy)
+         if (fabs(ry) <= fabs(iy))
            {
              double t = ry / iy;
              double d = iy * (1.0 + t * t);