-/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003 Free Software Foundation, Inc.
+/* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 Free Software Foundation, Inc.
*
* Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
* and Bellcore. See scm_divide.
\f
-static SCM abs_most_negative_fixnum;
static mpz_t z_negative_one;
\f
-static const char s_bignum[] = "bignum";
-
SCM_C_INLINE_KEYWORD SCM
scm_i_mkbig ()
{
else
{
/* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
- of that value for the denominator, as a bignum. */
- long abs_x = (x >= 0 ? x : -x);
- if (mpz_cmpabs_ui (SCM_I_BIG_MPZ (denominator), abs_x) == 0)
+ of that value for the denominator, as a bignum. Apart from
+ that case, abs(bignum) > abs(inum) so inum/bignum is not an
+ integer. */
+ if (x == SCM_MOST_NEGATIVE_FIXNUM
+ && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator),
+ - SCM_MOST_NEGATIVE_FIXNUM) == 0)
return SCM_MAKINUM(-1);
}
}
/* Some version of gcc on some old version of Linux used to crash when
trying to make Inf and NaN. */
-#if defined (SCO)
- double tmp = 1.0;
- guile_Inf = 1.0 / (tmp - tmp);
-#elif defined (__alpha__) && ! defined (linux)
+#ifdef INFINITY
+ /* C99 INFINITY, when available.
+ FIXME: The standard allows for INFINITY to be something that overflows
+ at compile time. We ought to have a configure test to check for that
+ before trying to use it. (But in practice we believe this is not a
+ problem on any system guile is likely to target.) */
+ guile_Inf = INFINITY;
+#elif HAVE_DINFINITY
+ /* OSF */
extern unsigned int DINFINITY[2];
guile_Inf = (*(X_CAST(double *, DINFINITY)));
#else
#if defined (HAVE_ISNAN)
-#if defined (__alpha__) && ! defined (linux)
+#ifdef NAN
+ /* C99 NAN, when available */
+ guile_NaN = NAN;
+#elif HAVE_DQNAN
+ /* OSF */
extern unsigned int DQNAN[2];
guile_NaN = (*(X_CAST(double *, DQNAN)));
#else
else if (SCM_BIGP (y))
{
if ((SCM_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
- && (scm_i_bigcmp (abs_most_negative_fixnum, y) == 0))
- /* Special case: x == fixnum-min && y == abs (fixnum-min) */
- return SCM_MAKINUM (-1);
+ && (mpz_cmp_ui (SCM_I_BIG_MPZ (y),
+ - SCM_MOST_NEGATIVE_FIXNUM) == 0))
+ {
+ /* Special case: x == fixnum-min && y == abs (fixnum-min) */
+ scm_remember_upto_here_1 (y);
+ return SCM_MAKINUM (-1);
+ }
else
return SCM_MAKINUM (0);
}
else if (SCM_BIGP (y))
{
if ((SCM_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
- && (scm_i_bigcmp (abs_most_negative_fixnum, y) == 0))
- /* Special case: x == fixnum-min && y == abs (fixnum-min) */
- return SCM_MAKINUM (0);
+ && (mpz_cmp_ui (SCM_I_BIG_MPZ (y),
+ - SCM_MOST_NEGATIVE_FIXNUM) == 0))
+ {
+ /* Special case: x == fixnum-min && y == abs (fixnum-min) */
+ scm_remember_upto_here_1 (y);
+ return SCM_MAKINUM (0);
+ }
else
return x;
}
else if (SCM_BIGP (y))
{
int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
-
- if (sgn_y == 0)
- scm_num_overflow (s_modulo);
- else
{
mpz_t z_x;
SCM result;
}
else if (SCM_BIGP (y))
{
- int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
- if (sgn_y == 0)
- scm_num_overflow (s_modulo);
- else
{
SCM result = scm_i_mkbig ();
int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
}
#undef FUNC_NAME
+/* returns 0 if IN is not an integer. OUT must already be
+ initialized. */
+static int
+coerce_to_big (SCM in, mpz_t out)
+{
+ if (SCM_BIGP (in))
+ mpz_set (out, SCM_I_BIG_MPZ (in));
+ else if (SCM_INUMP (in))
+ mpz_set_si (out, SCM_INUM (in));
+ else
+ return 0;
+
+ return 1;
+}
+
+SCM_DEFINE (scm_modulo_expt, "modulo-expt", 3, 0, 0,
+ (SCM n, SCM k, SCM m),
+ "Return @var{n} raised to the integer exponent\n"
+ "@var{k}, modulo @var{m}.\n"
+ "\n"
+ "@lisp\n"
+ "(modulo-expt 2 3 5)\n"
+ " @result{} 3\n"
+ "@end lisp")
+#define FUNC_NAME s_scm_modulo_expt
+{
+ mpz_t n_tmp;
+ mpz_t k_tmp;
+ mpz_t m_tmp;
+
+ /* There are two classes of error we might encounter --
+ 1) Math errors, which we'll report by calling scm_num_overflow,
+ and
+ 2) wrong-type errors, which of course we'll report by calling
+ SCM_WRONG_TYPE_ARG.
+ We don't report those errors immediately, however; instead we do
+ some cleanup first. These variables tell us which error (if
+ any) we should report after cleaning up.
+ */
+ int report_overflow = 0;
+
+ int position_of_wrong_type = 0;
+ SCM value_of_wrong_type = SCM_INUM0;
+
+ SCM result = SCM_UNDEFINED;
+
+ mpz_init (n_tmp);
+ mpz_init (k_tmp);
+ mpz_init (m_tmp);
+
+ if (SCM_EQ_P (m, SCM_INUM0))
+ {
+ report_overflow = 1;
+ goto cleanup;
+ }
+
+ if (!coerce_to_big (n, n_tmp))
+ {
+ value_of_wrong_type = n;
+ position_of_wrong_type = 1;
+ goto cleanup;
+ }
+
+ if (!coerce_to_big (k, k_tmp))
+ {
+ value_of_wrong_type = k;
+ position_of_wrong_type = 2;
+ goto cleanup;
+ }
+
+ if (!coerce_to_big (m, m_tmp))
+ {
+ value_of_wrong_type = m;
+ position_of_wrong_type = 3;
+ goto cleanup;
+ }
+
+ /* if the exponent K is negative, and we simply call mpz_powm, we
+ will get a divide-by-zero exception when an inverse 1/n mod m
+ doesn't exist (or is not unique). Since exceptions are hard to
+ handle, we'll attempt the inversion "by hand" -- that way, we get
+ a simple failure code, which is easy to handle. */
+
+ if (-1 == mpz_sgn (k_tmp))
+ {
+ if (!mpz_invert (n_tmp, n_tmp, m_tmp))
+ {
+ report_overflow = 1;
+ goto cleanup;
+ }
+ mpz_neg (k_tmp, k_tmp);
+ }
+
+ result = scm_i_mkbig ();
+ mpz_powm (SCM_I_BIG_MPZ (result),
+ n_tmp,
+ k_tmp,
+ m_tmp);
+
+ if (mpz_sgn (m_tmp) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)
+ mpz_add (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), m_tmp);
+
+ cleanup:
+ mpz_clear (m_tmp);
+ mpz_clear (k_tmp);
+ mpz_clear (n_tmp);
+
+ if (report_overflow)
+ scm_num_overflow (FUNC_NAME);
+
+ if (position_of_wrong_type)
+ SCM_WRONG_TYPE_ARG (position_of_wrong_type,
+ value_of_wrong_type);
+
+ return scm_i_normbig (result);
+}
+#undef FUNC_NAME
+
SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
(SCM n, SCM k),
"Return @var{n} raised to the non-negative integer exponent\n"
SCM
scm_num_eq_p (SCM x, SCM y)
{
+ again:
if (SCM_INUMP (x))
{
long xx = SCM_INUM (x);
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);
}
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);
}
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
}
+/* 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 by GMP in
+ 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."
SCM
scm_less_p (SCM x, SCM y)
{
+ again:
if (SCM_INUMP (x))
{
long xx = SCM_INUM (x);
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);
}
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);
}
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);
}
}
+/* scm_min and scm_max return an inexact when either argument is inexact, as
+ required by r5rs. On that basis, for exact/inexact combinations the
+ exact is converted to inexact to compare and possibly return. This is
+ unlike scm_less_p above which takes some trouble to preserve all bits in
+ its test, such trouble is not required for min and max. */
+
SCM_GPROC1 (s_max, "max", scm_tc7_asubr, scm_max, g_max);
/* "Return the maximum of all parameter values."
*/
{
if (SCM_UNBNDP (x))
SCM_WTA_DISPATCH_0 (g_max, s_max);
- else if (SCM_NUMBERP (x))
+ else if (SCM_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
return x;
else
SCM_WTA_DISPATCH_1 (g_max, x, SCM_ARG1, s_max);
}
else if (SCM_REALP (y))
{
- double yy = SCM_REAL_VALUE (y);
- int cmp;
- if (xisnan (yy))
- return y;
- cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), yy);
- scm_remember_upto_here_1 (x);
- return (cmp > 0) ? x : y;
+ /* if y==NaN then xx>yy is false, so we return the NaN y */
+ double xx, yy;
+ big_real:
+ xx = scm_i_big2dbl (x);
+ yy = SCM_REAL_VALUE (y);
+ return (xx > yy ? scm_make_real (xx) : y);
}
else if (SCM_FRACTIONP (y))
{
}
else if (SCM_BIGP (y))
{
- double xx = SCM_REAL_VALUE (x);
- int cmp;
- if (xisnan (xx))
- return x;
- cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), xx);
- scm_remember_upto_here_1 (y);
- return (cmp < 0) ? x : y;
+ SCM t = x; x = y; y = t;
+ goto big_real;
}
else if (SCM_REALP (y))
{
{
if (SCM_UNBNDP (x))
SCM_WTA_DISPATCH_0 (g_min, s_min);
- else if (SCM_NUMBERP (x))
+ else if (SCM_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
return x;
else
SCM_WTA_DISPATCH_1 (g_min, x, SCM_ARG1, s_min);
}
else if (SCM_REALP (y))
{
- double yy = SCM_REAL_VALUE (y);
- int cmp;
- if (xisnan (yy))
- return y;
- cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), yy);
- scm_remember_upto_here_1 (x);
- return (cmp > 0) ? y : x;
+ /* if y==NaN then xx<yy is false, so we return the NaN y */
+ double xx, yy;
+ big_real:
+ xx = scm_i_big2dbl (x);
+ yy = SCM_REAL_VALUE (y);
+ return (xx < yy ? scm_make_real (xx) : y);
}
else if (SCM_FRACTIONP (y))
{
}
else if (SCM_BIGP (y))
{
- double xx = SCM_REAL_VALUE (x);
- int cmp;
- if (xisnan (xx))
- return x;
- cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), xx);
- scm_remember_upto_here_1 (y);
- return (cmp < 0) ? y : x;
+ SCM t = x; x = y; y = t;
+ goto big_real;
}
else if (SCM_REALP (y))
{
void
scm_init_numbers ()
{
- abs_most_negative_fixnum = scm_i_long2big (- SCM_MOST_NEGATIVE_FIXNUM);
- scm_permanent_object (abs_most_negative_fixnum);
-
mpz_init_set_si (z_negative_one, -1);
/* It may be possible to tune the performance of some algorithms by using
HCoop / bpt / guile.git / blobdiff