1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc() */
55 #include "libguile/_scm.h"
56 #include "libguile/feature.h"
57 #include "libguile/ports.h"
58 #include "libguile/root.h"
59 #include "libguile/smob.h"
60 #include "libguile/strings.h"
62 #include "libguile/validate.h"
63 #include "libguile/numbers.h"
64 #include "libguile/deprecation.h"
66 #include "libguile/eq.h"
71 Wonder if this might be faster for some of our code? A switch on
72 the numtag would jump directly to the right case, and the
73 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
75 #define SCM_I_NUMTAG_NOTNUM 0
76 #define SCM_I_NUMTAG_INUM 1
77 #define SCM_I_NUMTAG_BIG scm_tc16_big
78 #define SCM_I_NUMTAG_REAL scm_tc16_real
79 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
80 #define SCM_I_NUMTAG(x) \
81 (SCM_INUMP(x) ? SCM_I_NUMTAG_INUM \
82 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
83 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
84 : SCM_I_NUMTAG_NOTNUM)))
86 /* the macro above will not work as is with fractions */
89 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
91 /* FLOBUFLEN is the maximum number of characters neccessary for the
92 * printed or scm_string representation of an inexact number.
94 #define FLOBUFLEN (10+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
97 #if ! defined (HAVE_ISNAN)
102 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
105 #if ! defined (HAVE_ISINF)
110 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
117 /* mpz_cmp_d only recognises infinities in gmp 4.2 and up.
118 For prior versions use an explicit check here. */
119 #if __GNU_MP_VERSION < 4 \
120 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
121 #define xmpz_cmp_d(z, d) \
122 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
124 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
130 #if defined (HAVE_ISINF)
132 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
133 return (! (finite (x
) || isnan (x
)));
142 #if defined (HAVE_ISNAN)
151 static SCM abs_most_negative_fixnum
;
152 static mpz_t z_negative_one
;
156 SCM_C_INLINE_KEYWORD SCM
159 /* Return a newly created bignum. */
160 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
161 mpz_init (SCM_I_BIG_MPZ (z
));
165 SCM_C_INLINE_KEYWORD
static SCM
166 scm_i_clonebig (SCM src_big
, int same_sign_p
)
168 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
169 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
170 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
172 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
176 SCM_C_INLINE_KEYWORD
int
177 scm_i_bigcmp (SCM x
, SCM y
)
179 /* Return neg if x < y, pos if x > y, and 0 if x == y */
180 /* presume we already know x and y are bignums */
181 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
182 scm_remember_upto_here_2 (x
, y
);
186 SCM_C_INLINE_KEYWORD SCM
187 scm_i_dbl2big (double d
)
189 /* results are only defined if d is an integer */
190 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
191 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
195 /* Convert a integer in double representation to a SCM number. */
197 SCM_C_INLINE_KEYWORD SCM
198 scm_i_dbl2num (double u
)
200 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
201 powers of 2, so there's no rounding when making "double" values
202 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
203 get rounded on a 64-bit machine, hence the "+1".
205 The use of floor() to force to an integer value ensures we get a
206 "numerically closest" value without depending on how a
207 double->long cast or how mpz_set_d will round. For reference,
208 double->long probably follows the hardware rounding mode,
209 mpz_set_d truncates towards zero. */
211 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
212 representable as a double? */
214 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
215 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
216 return SCM_MAKINUM ((long) u
);
218 return scm_i_dbl2big (u
);
221 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
222 with R5RS exact->inexact.
224 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
225 (ie. it truncates towards zero), then adjust to get the closest double by
226 examining the next lower bit and adding 1 if necessary.
228 Note that bignums exactly half way between representable doubles are
229 rounded to the next higher absolute value (ie. away from zero). This
230 seems like an adequate interpretation of R5RS "numerically closest", and
231 it's easier and faster than a full "nearest-even" style.
233 The bit test is done on the absolute value of the mpz_t, which means we
234 must use mpz_getlimbn. mpz_tstbit is not right, it treats negatives as
237 Prior to GMP 4.2, the rounding done by mpz_get_d was unspecified. It
238 happened to follow the hardware rounding mode, but on the absolute value
239 of its operand. This is not what we want, so we put the high
240 DBL_MANT_DIG bits into a temporary. This extra init/clear is a slowdown,
241 but doesn't matter too much since it's only for older GMP. */
244 scm_i_big2dbl (SCM b
)
249 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
251 #if __GNU_MP_VERSION < 4 \
252 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
254 /* GMP prior to 4.2, force truncate towards zero */
256 if (bits
> DBL_MANT_DIG
)
258 size_t shift
= bits
- DBL_MANT_DIG
;
259 mpz_init2 (tmp
, DBL_MANT_DIG
);
260 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
261 result
= ldexp (mpz_get_d (tmp
), shift
);
266 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
271 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
274 if (bits
> DBL_MANT_DIG
)
276 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
277 /* test bit number "pos" in absolute value */
278 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
279 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
281 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
285 scm_remember_upto_here_1 (b
);
289 SCM_C_INLINE_KEYWORD SCM
290 scm_i_normbig (SCM b
)
292 /* convert a big back to a fixnum if it'll fit */
293 /* presume b is a bignum */
294 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
296 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
297 if (SCM_FIXABLE (val
))
298 b
= SCM_MAKINUM (val
);
303 static SCM_C_INLINE_KEYWORD SCM
304 scm_i_mpz2num (mpz_t b
)
306 /* convert a mpz number to a SCM number. */
307 if (mpz_fits_slong_p (b
))
309 long val
= mpz_get_si (b
);
310 if (SCM_FIXABLE (val
))
311 return SCM_MAKINUM (val
);
315 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
316 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
321 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
322 static SCM
scm_divide2real (SCM x
, SCM y
);
325 scm_make_ratio (SCM numerator
, SCM denominator
)
326 #define FUNC_NAME "make-ratio"
328 /* First make sure the arguments are proper.
330 if (SCM_INUMP (denominator
))
332 if (SCM_EQ_P (denominator
, SCM_INUM0
))
333 scm_num_overflow ("make-ratio");
334 if (SCM_EQ_P (denominator
, SCM_MAKINUM(1)))
339 if (!(SCM_BIGP(denominator
)))
340 SCM_WRONG_TYPE_ARG (2, denominator
);
342 if (!SCM_INUMP (numerator
) && !SCM_BIGP (numerator
))
343 SCM_WRONG_TYPE_ARG (1, numerator
);
345 /* Then flip signs so that the denominator is positive.
347 if (SCM_NFALSEP (scm_negative_p (denominator
)))
349 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
350 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
353 /* Now consider for each of the four fixnum/bignum combinations
354 whether the rational number is really an integer.
356 if (SCM_INUMP (numerator
))
358 long x
= SCM_INUM (numerator
);
359 if (SCM_EQ_P (numerator
, SCM_INUM0
))
361 if (SCM_INUMP (denominator
))
364 y
= SCM_INUM (denominator
);
366 return SCM_MAKINUM(1);
368 return SCM_MAKINUM (x
/ y
);
372 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
373 of that value for the denominator, as a bignum. */
374 long abs_x
= (x
>= 0 ? x
: -x
);
375 if (mpz_cmpabs_ui (SCM_I_BIG_MPZ (denominator
), abs_x
) == 0)
376 return SCM_MAKINUM(-1);
379 else if (SCM_BIGP (numerator
))
381 if (SCM_INUMP (denominator
))
383 long yy
= SCM_INUM (denominator
);
384 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
385 return scm_divide (numerator
, denominator
);
389 if (SCM_EQ_P (numerator
, denominator
))
390 return SCM_MAKINUM(1);
391 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
392 SCM_I_BIG_MPZ (denominator
)))
393 return scm_divide(numerator
, denominator
);
397 /* No, it's a proper fraction.
399 return scm_double_cell (scm_tc16_fraction
,
400 SCM_UNPACK (numerator
),
401 SCM_UNPACK (denominator
), 0);
405 static void scm_i_fraction_reduce (SCM z
)
407 if (!(SCM_FRACTION_REDUCED (z
)))
410 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
411 if (!(SCM_EQ_P (divisor
, SCM_MAKINUM(1))))
414 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
415 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
417 SCM_FRACTION_REDUCED_SET (z
);
422 scm_i_fraction2double (SCM z
)
424 return scm_num2dbl (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
425 SCM_FRACTION_DENOMINATOR (z
)),
429 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
431 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
433 #define FUNC_NAME s_scm_exact_p
439 if (SCM_FRACTIONP (x
))
443 SCM_WRONG_TYPE_ARG (1, x
);
448 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
450 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
452 #define FUNC_NAME s_scm_odd_p
456 long val
= SCM_INUM (n
);
457 return SCM_BOOL ((val
& 1L) != 0);
459 else if (SCM_BIGP (n
))
461 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
462 scm_remember_upto_here_1 (n
);
463 return SCM_BOOL (odd_p
);
465 else if (!SCM_FALSEP (scm_inf_p (n
)))
467 else if (SCM_REALP (n
))
469 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
475 SCM_WRONG_TYPE_ARG (1, n
);
478 SCM_WRONG_TYPE_ARG (1, n
);
483 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
485 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
487 #define FUNC_NAME s_scm_even_p
491 long val
= SCM_INUM (n
);
492 return SCM_BOOL ((val
& 1L) == 0);
494 else if (SCM_BIGP (n
))
496 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
497 scm_remember_upto_here_1 (n
);
498 return SCM_BOOL (even_p
);
500 else if (!SCM_FALSEP (scm_inf_p (n
)))
502 else if (SCM_REALP (n
))
504 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
510 SCM_WRONG_TYPE_ARG (1, n
);
513 SCM_WRONG_TYPE_ARG (1, n
);
517 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
519 "Return @code{#t} if @var{n} is infinite, @code{#f}\n"
521 #define FUNC_NAME s_scm_inf_p
524 return SCM_BOOL (xisinf (SCM_REAL_VALUE (n
)));
525 else if (SCM_COMPLEXP (n
))
526 return SCM_BOOL (xisinf (SCM_COMPLEX_REAL (n
))
527 || xisinf (SCM_COMPLEX_IMAG (n
)));
533 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
535 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
537 #define FUNC_NAME s_scm_nan_p
540 return SCM_BOOL (xisnan (SCM_REAL_VALUE (n
)));
541 else if (SCM_COMPLEXP (n
))
542 return SCM_BOOL (xisnan (SCM_COMPLEX_REAL (n
))
543 || xisnan (SCM_COMPLEX_IMAG (n
)));
549 /* Guile's idea of infinity. */
550 static double guile_Inf
;
552 /* Guile's idea of not a number. */
553 static double guile_NaN
;
556 guile_ieee_init (void)
558 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
560 /* Some version of gcc on some old version of Linux used to crash when
561 trying to make Inf and NaN. */
564 /* C99 INFINITY, when available.
565 FIXME: The standard allows for INFINITY to be something that overflows
566 at compile time. We ought to have a configure test to check for that
567 before trying to use it. (But in practice we believe this is not a
568 problem on any system guile is likely to target.) */
569 guile_Inf
= INFINITY
;
572 extern unsigned int DINFINITY
[2];
573 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
580 if (guile_Inf
== tmp
)
588 #if defined (HAVE_ISNAN)
591 /* C99 NAN, when available */
595 extern unsigned int DQNAN
[2];
596 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
598 guile_NaN
= guile_Inf
/ guile_Inf
;
604 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
607 #define FUNC_NAME s_scm_inf
609 static int initialized
= 0;
615 return scm_make_real (guile_Inf
);
619 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
622 #define FUNC_NAME s_scm_nan
624 static int initialized
= 0;
630 return scm_make_real (guile_NaN
);
635 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
637 "Return the absolute value of @var{x}.")
642 long int xx
= SCM_INUM (x
);
645 else if (SCM_POSFIXABLE (-xx
))
646 return SCM_MAKINUM (-xx
);
648 return scm_i_long2big (-xx
);
650 else if (SCM_BIGP (x
))
652 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
654 return scm_i_clonebig (x
, 0);
658 else if (SCM_REALP (x
))
660 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
661 double xx
= SCM_REAL_VALUE (x
);
663 return scm_make_real (-xx
);
667 else if (SCM_FRACTIONP (x
))
669 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
671 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
672 SCM_FRACTION_DENOMINATOR (x
));
675 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
680 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
681 /* "Return the quotient of the numbers @var{x} and @var{y}."
684 scm_quotient (SCM x
, SCM y
)
688 long xx
= SCM_INUM (x
);
691 long yy
= SCM_INUM (y
);
693 scm_num_overflow (s_quotient
);
698 return SCM_MAKINUM (z
);
700 return scm_i_long2big (z
);
703 else if (SCM_BIGP (y
))
705 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
706 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
707 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
708 return SCM_MAKINUM (-1);
710 return SCM_MAKINUM (0);
713 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
715 else if (SCM_BIGP (x
))
719 long yy
= SCM_INUM (y
);
721 scm_num_overflow (s_quotient
);
726 SCM result
= scm_i_mkbig ();
729 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
732 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
735 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
736 scm_remember_upto_here_1 (x
);
737 return scm_i_normbig (result
);
740 else if (SCM_BIGP (y
))
742 SCM result
= scm_i_mkbig ();
743 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
746 scm_remember_upto_here_2 (x
, y
);
747 return scm_i_normbig (result
);
750 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
753 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
756 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
757 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
759 * "(remainder 13 4) @result{} 1\n"
760 * "(remainder -13 4) @result{} -1\n"
764 scm_remainder (SCM x
, SCM y
)
770 long yy
= SCM_INUM (y
);
772 scm_num_overflow (s_remainder
);
775 long z
= SCM_INUM (x
) % yy
;
776 return SCM_MAKINUM (z
);
779 else if (SCM_BIGP (y
))
781 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
782 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
783 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
784 return SCM_MAKINUM (0);
789 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
791 else if (SCM_BIGP (x
))
795 long yy
= SCM_INUM (y
);
797 scm_num_overflow (s_remainder
);
800 SCM result
= scm_i_mkbig ();
803 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
804 scm_remember_upto_here_1 (x
);
805 return scm_i_normbig (result
);
808 else if (SCM_BIGP (y
))
810 SCM result
= scm_i_mkbig ();
811 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
814 scm_remember_upto_here_2 (x
, y
);
815 return scm_i_normbig (result
);
818 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
821 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
825 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
826 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
828 * "(modulo 13 4) @result{} 1\n"
829 * "(modulo -13 4) @result{} 3\n"
833 scm_modulo (SCM x
, SCM y
)
837 long xx
= SCM_INUM (x
);
840 long yy
= SCM_INUM (y
);
842 scm_num_overflow (s_modulo
);
845 /* FIXME: I think this may be a bug on some arches -- results
846 of % with negative second arg are undefined... */
864 return SCM_MAKINUM (result
);
867 else if (SCM_BIGP (y
))
869 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
872 scm_num_overflow (s_modulo
);
880 SCM pos_y
= scm_i_clonebig (y
, 0);
881 /* do this after the last scm_op */
882 mpz_init_set_si (z_x
, xx
);
883 result
= pos_y
; /* re-use this bignum */
884 mpz_mod (SCM_I_BIG_MPZ (result
),
886 SCM_I_BIG_MPZ (pos_y
));
887 scm_remember_upto_here_1 (pos_y
);
891 result
= scm_i_mkbig ();
892 /* do this after the last scm_op */
893 mpz_init_set_si (z_x
, xx
);
894 mpz_mod (SCM_I_BIG_MPZ (result
),
897 scm_remember_upto_here_1 (y
);
900 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
901 mpz_add (SCM_I_BIG_MPZ (result
),
903 SCM_I_BIG_MPZ (result
));
904 scm_remember_upto_here_1 (y
);
905 /* and do this before the next one */
907 return scm_i_normbig (result
);
911 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
913 else if (SCM_BIGP (x
))
917 long yy
= SCM_INUM (y
);
919 scm_num_overflow (s_modulo
);
922 SCM result
= scm_i_mkbig ();
923 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
925 (yy
< 0) ? - yy
: yy
);
926 scm_remember_upto_here_1 (x
);
927 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
928 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
929 SCM_I_BIG_MPZ (result
),
931 return scm_i_normbig (result
);
934 else if (SCM_BIGP (y
))
936 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
938 scm_num_overflow (s_modulo
);
941 SCM result
= scm_i_mkbig ();
942 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
943 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
944 mpz_mod (SCM_I_BIG_MPZ (result
),
946 SCM_I_BIG_MPZ (pos_y
));
948 scm_remember_upto_here_1 (x
);
949 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
950 mpz_add (SCM_I_BIG_MPZ (result
),
952 SCM_I_BIG_MPZ (result
));
953 scm_remember_upto_here_2 (y
, pos_y
);
954 return scm_i_normbig (result
);
958 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
961 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
964 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
965 /* "Return the greatest common divisor of all arguments.\n"
966 * "If called without arguments, 0 is returned."
969 scm_gcd (SCM x
, SCM y
)
972 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
978 long xx
= SCM_INUM (x
);
979 long yy
= SCM_INUM (y
);
980 long u
= xx
< 0 ? -xx
: xx
;
981 long v
= yy
< 0 ? -yy
: yy
;
991 /* Determine a common factor 2^k */
992 while (!(1 & (u
| v
)))
998 /* Now, any factor 2^n can be eliminated */
1018 return (SCM_POSFIXABLE (result
)
1019 ? SCM_MAKINUM (result
)
1020 : scm_i_long2big (result
));
1022 else if (SCM_BIGP (y
))
1024 SCM result
= scm_i_mkbig ();
1025 SCM mx
= scm_i_mkbig ();
1026 mpz_set_si (SCM_I_BIG_MPZ (mx
), SCM_INUM (x
));
1027 scm_remember_upto_here_1 (x
);
1028 mpz_gcd (SCM_I_BIG_MPZ (result
),
1031 scm_remember_upto_here_2 (mx
, y
);
1032 return scm_i_normbig (result
);
1035 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1037 else if (SCM_BIGP (x
))
1041 unsigned long result
;
1042 long yy
= SCM_INUM (y
);
1047 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1048 scm_remember_upto_here_1 (x
);
1049 return (SCM_POSFIXABLE (result
)
1050 ? SCM_MAKINUM (result
)
1051 : scm_ulong2num (result
));
1053 else if (SCM_BIGP (y
))
1055 SCM result
= scm_i_mkbig ();
1056 mpz_gcd (SCM_I_BIG_MPZ (result
),
1059 scm_remember_upto_here_2 (x
, y
);
1060 return scm_i_normbig (result
);
1063 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1066 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1069 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1070 /* "Return the least common multiple of the arguments.\n"
1071 * "If called without arguments, 1 is returned."
1074 scm_lcm (SCM n1
, SCM n2
)
1076 if (SCM_UNBNDP (n2
))
1078 if (SCM_UNBNDP (n1
))
1079 return SCM_MAKINUM (1L);
1080 n2
= SCM_MAKINUM (1L);
1083 SCM_GASSERT2 (SCM_INUMP (n1
) || SCM_BIGP (n1
),
1084 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1085 SCM_GASSERT2 (SCM_INUMP (n2
) || SCM_BIGP (n2
),
1086 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1092 SCM d
= scm_gcd (n1
, n2
);
1093 if (SCM_EQ_P (d
, SCM_INUM0
))
1096 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1100 /* inum n1, big n2 */
1103 SCM result
= scm_i_mkbig ();
1104 long nn1
= SCM_INUM (n1
);
1105 if (nn1
== 0) return SCM_INUM0
;
1106 if (nn1
< 0) nn1
= - nn1
;
1107 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1108 scm_remember_upto_here_1 (n2
);
1123 SCM result
= scm_i_mkbig ();
1124 mpz_lcm(SCM_I_BIG_MPZ (result
),
1126 SCM_I_BIG_MPZ (n2
));
1127 scm_remember_upto_here_2(n1
, n2
);
1128 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1134 #ifndef scm_long2num
1135 #define SCM_LOGOP_RETURN(x) scm_ulong2num(x)
1137 #define SCM_LOGOP_RETURN(x) SCM_MAKINUM(x)
1140 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1145 + + + x (map digit:logand X Y)
1146 + - + x (map digit:logand X (lognot (+ -1 Y)))
1147 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1148 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1153 + + + (map digit:logior X Y)
1154 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1155 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1156 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1161 + + + (map digit:logxor X Y)
1162 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1163 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1164 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1169 + + (any digit:logand X Y)
1170 + - (any digit:logand X (lognot (+ -1 Y)))
1171 - + (any digit:logand (lognot (+ -1 X)) Y)
1176 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1178 "Return the bitwise AND of the integer arguments.\n\n"
1180 "(logand) @result{} -1\n"
1181 "(logand 7) @result{} 7\n"
1182 "(logand #b111 #b011 #b001) @result{} 1\n"
1184 #define FUNC_NAME s_scm_logand
1188 if (SCM_UNBNDP (n2
))
1190 if (SCM_UNBNDP (n1
))
1191 return SCM_MAKINUM (-1);
1192 else if (!SCM_NUMBERP (n1
))
1193 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1194 else if (SCM_NUMBERP (n1
))
1197 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1202 nn1
= SCM_INUM (n1
);
1205 long nn2
= SCM_INUM (n2
);
1206 return SCM_MAKINUM (nn1
& nn2
);
1208 else if SCM_BIGP (n2
)
1214 SCM result_z
= scm_i_mkbig ();
1216 mpz_init_set_si (nn1_z
, nn1
);
1217 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1218 scm_remember_upto_here_1 (n2
);
1220 return scm_i_normbig (result_z
);
1224 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1226 else if (SCM_BIGP (n1
))
1231 nn1
= SCM_INUM (n1
);
1234 else if (SCM_BIGP (n2
))
1236 SCM result_z
= scm_i_mkbig ();
1237 mpz_and (SCM_I_BIG_MPZ (result_z
),
1239 SCM_I_BIG_MPZ (n2
));
1240 scm_remember_upto_here_2 (n1
, n2
);
1241 return scm_i_normbig (result_z
);
1244 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1247 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1252 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1254 "Return the bitwise OR of the integer arguments.\n\n"
1256 "(logior) @result{} 0\n"
1257 "(logior 7) @result{} 7\n"
1258 "(logior #b000 #b001 #b011) @result{} 3\n"
1260 #define FUNC_NAME s_scm_logior
1264 if (SCM_UNBNDP (n2
))
1266 if (SCM_UNBNDP (n1
))
1268 else if (SCM_NUMBERP (n1
))
1271 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1276 nn1
= SCM_INUM (n1
);
1279 long nn2
= SCM_INUM (n2
);
1280 return SCM_MAKINUM (nn1
| nn2
);
1282 else if (SCM_BIGP (n2
))
1288 SCM result_z
= scm_i_mkbig ();
1290 mpz_init_set_si (nn1_z
, nn1
);
1291 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1292 scm_remember_upto_here_1 (n2
);
1298 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1300 else if (SCM_BIGP (n1
))
1305 nn1
= SCM_INUM (n1
);
1308 else if (SCM_BIGP (n2
))
1310 SCM result_z
= scm_i_mkbig ();
1311 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1313 SCM_I_BIG_MPZ (n2
));
1314 scm_remember_upto_here_2 (n1
, n2
);
1318 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1321 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1326 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1328 "Return the bitwise XOR of the integer arguments. A bit is\n"
1329 "set in the result if it is set in an odd number of arguments.\n"
1331 "(logxor) @result{} 0\n"
1332 "(logxor 7) @result{} 7\n"
1333 "(logxor #b000 #b001 #b011) @result{} 2\n"
1334 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1336 #define FUNC_NAME s_scm_logxor
1340 if (SCM_UNBNDP (n2
))
1342 if (SCM_UNBNDP (n1
))
1344 else if (SCM_NUMBERP (n1
))
1347 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1352 nn1
= SCM_INUM (n1
);
1355 long nn2
= SCM_INUM (n2
);
1356 return SCM_MAKINUM (nn1
^ nn2
);
1358 else if (SCM_BIGP (n2
))
1362 SCM result_z
= scm_i_mkbig ();
1364 mpz_init_set_si (nn1_z
, nn1
);
1365 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1366 scm_remember_upto_here_1 (n2
);
1368 return scm_i_normbig (result_z
);
1372 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1374 else if (SCM_BIGP (n1
))
1379 nn1
= SCM_INUM (n1
);
1382 else if (SCM_BIGP (n2
))
1384 SCM result_z
= scm_i_mkbig ();
1385 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1387 SCM_I_BIG_MPZ (n2
));
1388 scm_remember_upto_here_2 (n1
, n2
);
1389 return scm_i_normbig (result_z
);
1392 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1395 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1400 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1403 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1404 "(logtest #b0100 #b1011) @result{} #f\n"
1405 "(logtest #b0100 #b0111) @result{} #t\n"
1407 #define FUNC_NAME s_scm_logtest
1416 long nk
= SCM_INUM (k
);
1417 return SCM_BOOL (nj
& nk
);
1419 else if (SCM_BIGP (k
))
1427 mpz_init_set_si (nj_z
, nj
);
1428 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1429 scm_remember_upto_here_1 (k
);
1430 result
= SCM_BOOL (mpz_sgn (nj_z
) != 0);
1436 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1438 else if (SCM_BIGP (j
))
1446 else if (SCM_BIGP (k
))
1450 mpz_init (result_z
);
1454 scm_remember_upto_here_2 (j
, k
);
1455 result
= SCM_BOOL (mpz_sgn (result_z
) != 0);
1456 mpz_clear (result_z
);
1460 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1463 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1468 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1471 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1472 "(logbit? 0 #b1101) @result{} #t\n"
1473 "(logbit? 1 #b1101) @result{} #f\n"
1474 "(logbit? 2 #b1101) @result{} #t\n"
1475 "(logbit? 3 #b1101) @result{} #t\n"
1476 "(logbit? 4 #b1101) @result{} #f\n"
1478 #define FUNC_NAME s_scm_logbit_p
1480 unsigned long int iindex
;
1482 SCM_VALIDATE_INUM_MIN (SCM_ARG1
, index
, 0);
1483 iindex
= (unsigned long int) SCM_INUM (index
);
1486 return SCM_BOOL ((1L << iindex
) & SCM_INUM (j
));
1487 else if (SCM_BIGP (j
))
1489 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1490 scm_remember_upto_here_1 (j
);
1491 return SCM_BOOL (val
);
1494 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1499 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1501 "Return the integer which is the ones-complement of the integer\n"
1505 "(number->string (lognot #b10000000) 2)\n"
1506 " @result{} \"-10000001\"\n"
1507 "(number->string (lognot #b0) 2)\n"
1508 " @result{} \"-1\"\n"
1510 #define FUNC_NAME s_scm_lognot
1512 if (SCM_INUMP (n
)) {
1513 /* No overflow here, just need to toggle all the bits making up the inum.
1514 Enhancement: No need to strip the tag and add it back, could just xor
1515 a block of 1 bits, if that worked with the various debug versions of
1517 return SCM_MAKINUM (~ SCM_INUM (n
));
1519 } else if (SCM_BIGP (n
)) {
1520 SCM result
= scm_i_mkbig ();
1521 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1522 scm_remember_upto_here_1 (n
);
1526 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1531 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1533 "Return @var{n} raised to the non-negative integer exponent\n"
1537 "(integer-expt 2 5)\n"
1539 "(integer-expt -3 3)\n"
1542 #define FUNC_NAME s_scm_integer_expt
1545 SCM z_i2
= SCM_BOOL_F
;
1547 SCM acc
= SCM_MAKINUM (1L);
1549 /* 0^0 == 1 according to R5RS */
1550 if (SCM_EQ_P (n
, SCM_INUM0
) || SCM_EQ_P (n
, acc
))
1551 return SCM_FALSEP (scm_zero_p(k
)) ? n
: acc
;
1552 else if (SCM_EQ_P (n
, SCM_MAKINUM (-1L)))
1553 return SCM_FALSEP (scm_even_p (k
)) ? n
: acc
;
1557 else if (SCM_BIGP (k
))
1559 z_i2
= scm_i_clonebig (k
, 1);
1560 scm_remember_upto_here_1 (k
);
1563 else if (SCM_REALP (k
))
1565 double r
= SCM_REAL_VALUE (k
);
1567 SCM_WRONG_TYPE_ARG (2, k
);
1568 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1570 z_i2
= scm_i_mkbig ();
1571 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1580 SCM_WRONG_TYPE_ARG (2, k
);
1584 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1586 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1587 n
= scm_divide (n
, SCM_UNDEFINED
);
1591 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1595 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1597 return scm_product (acc
, n
);
1599 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1600 acc
= scm_product (acc
, n
);
1601 n
= scm_product (n
, n
);
1602 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1610 n
= scm_divide (n
, SCM_UNDEFINED
);
1617 return scm_product (acc
, n
);
1619 acc
= scm_product (acc
, n
);
1620 n
= scm_product (n
, n
);
1627 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1629 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1630 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1632 "This is effectively a multiplication by 2^@var{cnt}}, and when\n"
1633 "@var{cnt} is negative it's a division, rounded towards negative\n"
1634 "infinity. (Note that this is not the same rounding as\n"
1635 "@code{quotient} does.)\n"
1637 "With @var{n} viewed as an infinite precision twos complement,\n"
1638 "@code{ash} means a left shift introducing zero bits, or a right\n"
1639 "shift dropping bits.\n"
1642 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1643 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1645 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1646 "(ash -23 -2) @result{} -6\n"
1648 #define FUNC_NAME s_scm_ash
1652 SCM_VALIDATE_INUM (2, cnt
);
1654 bits_to_shift
= SCM_INUM (cnt
);
1656 if (bits_to_shift
< 0)
1658 /* Shift right by abs(cnt) bits. This is realized as a division
1659 by div:=2^abs(cnt). However, to guarantee the floor
1660 rounding, negative values require some special treatment.
1662 SCM div
= scm_integer_expt (SCM_MAKINUM (2),
1663 SCM_MAKINUM (-bits_to_shift
));
1665 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1666 if (SCM_FALSEP (scm_negative_p (n
)))
1667 return scm_quotient (n
, div
);
1669 return scm_sum (SCM_MAKINUM (-1L),
1670 scm_quotient (scm_sum (SCM_MAKINUM (1L), n
), div
));
1673 /* Shift left is done by multiplication with 2^CNT */
1674 return scm_product (n
, scm_integer_expt (SCM_MAKINUM (2), cnt
));
1679 #define MIN(x,y) ((x) < (y) ? (x) : (y))
1681 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1682 (SCM n
, SCM start
, SCM end
),
1683 "Return the integer composed of the @var{start} (inclusive)\n"
1684 "through @var{end} (exclusive) bits of @var{n}. The\n"
1685 "@var{start}th bit becomes the 0-th bit in the result.\n"
1688 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1689 " @result{} \"1010\"\n"
1690 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1691 " @result{} \"10110\"\n"
1693 #define FUNC_NAME s_scm_bit_extract
1695 unsigned long int istart
, iend
, bits
;
1696 SCM_VALIDATE_INUM_MIN_COPY (2, start
,0, istart
);
1697 SCM_VALIDATE_INUM_MIN_COPY (3, end
, 0, iend
);
1698 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1700 /* how many bits to keep */
1701 bits
= iend
- istart
;
1705 long int in
= SCM_INUM (n
);
1707 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1708 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in".
1709 FIXME: This shift relies on signed right shifts being arithmetic,
1710 which is not guaranteed by C99. */
1711 in
>>= MIN (istart
, SCM_I_FIXNUM_BIT
-1);
1713 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1715 /* Since we emulate two's complement encoded numbers, this
1716 * special case requires us to produce a result that has
1717 * more bits than can be stored in a fixnum.
1719 SCM result
= scm_i_long2big (in
);
1720 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1725 /* mask down to requisite bits */
1726 bits
= MIN (bits
, SCM_I_FIXNUM_BIT
);
1727 return SCM_MAKINUM (in
& ((1L << bits
) - 1));
1729 else if (SCM_BIGP (n
))
1734 result
= SCM_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1738 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1739 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1740 such bits into a ulong. */
1741 result
= scm_i_mkbig ();
1742 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1743 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1744 result
= scm_i_normbig (result
);
1746 scm_remember_upto_here_1 (n
);
1750 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1755 static const char scm_logtab
[] = {
1756 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1759 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1761 "Return the number of bits in integer @var{n}. If integer is\n"
1762 "positive, the 1-bits in its binary representation are counted.\n"
1763 "If negative, the 0-bits in its two's-complement binary\n"
1764 "representation are counted. If 0, 0 is returned.\n"
1767 "(logcount #b10101010)\n"
1774 #define FUNC_NAME s_scm_logcount
1778 unsigned long int c
= 0;
1779 long int nn
= SCM_INUM (n
);
1784 c
+= scm_logtab
[15 & nn
];
1787 return SCM_MAKINUM (c
);
1789 else if (SCM_BIGP (n
))
1791 unsigned long count
;
1792 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1793 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1795 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1796 scm_remember_upto_here_1 (n
);
1797 return SCM_MAKINUM (count
);
1800 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1805 static const char scm_ilentab
[] = {
1806 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1810 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1812 "Return the number of bits necessary to represent @var{n}.\n"
1815 "(integer-length #b10101010)\n"
1817 "(integer-length 0)\n"
1819 "(integer-length #b1111)\n"
1822 #define FUNC_NAME s_scm_integer_length
1826 unsigned long int c
= 0;
1828 long int nn
= SCM_INUM (n
);
1834 l
= scm_ilentab
[15 & nn
];
1837 return SCM_MAKINUM (c
- 4 + l
);
1839 else if (SCM_BIGP (n
))
1841 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1842 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1843 1 too big, so check for that and adjust. */
1844 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1845 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1846 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1847 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1849 scm_remember_upto_here_1 (n
);
1850 return SCM_MAKINUM (size
);
1853 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1857 /*** NUMBERS -> STRINGS ***/
1859 static const double fx
[] =
1860 { 0.0, 5e-1, 5e-2, 5e-3, 5e-4, 5e-5,
1861 5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
1862 5e-11, 5e-12, 5e-13, 5e-14, 5e-15,
1863 5e-16, 5e-17, 5e-18, 5e-19, 5e-20};
1866 idbl2str (double f
, char *a
)
1868 int efmt
, dpt
, d
, i
, wp
= scm_dblprec
;
1874 #ifdef HAVE_COPYSIGN
1875 double sgn
= copysign (1.0, f
);
1881 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
1887 strcpy (a
, "-inf.0");
1889 strcpy (a
, "+inf.0");
1892 else if (xisnan (f
))
1894 strcpy (a
, "+nan.0");
1904 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
1905 make-uniform-vector, from causing infinite loops. */
1909 if (exp
-- < DBL_MIN_10_EXP
)
1920 if (exp
++ > DBL_MAX_10_EXP
)
1940 if (f
+ fx
[wp
] >= 10.0)
1947 dpt
= (exp
+ 9999) % 3;
1951 efmt
= (exp
< -3) || (exp
> wp
+ 2);
1976 if (f
+ fx
[wp
] >= 1.0)
1990 if ((dpt
> 4) && (exp
> 6))
1992 d
= (a
[0] == '-' ? 2 : 1);
1993 for (i
= ch
++; i
> d
; i
--)
2006 if (a
[ch
- 1] == '.')
2007 a
[ch
++] = '0'; /* trailing zero */
2016 for (i
= 10; i
<= exp
; i
*= 10);
2017 for (i
/= 10; i
; i
/= 10)
2019 a
[ch
++] = exp
/ i
+ '0';
2028 iflo2str (SCM flt
, char *str
)
2031 if (SCM_REALP (flt
))
2032 i
= idbl2str (SCM_REAL_VALUE (flt
), str
);
2035 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
);
2036 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2038 double imag
= SCM_COMPLEX_IMAG (flt
);
2039 /* Don't output a '+' for negative numbers or for Inf and
2040 NaN. They will provide their own sign. */
2041 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2043 i
+= idbl2str (imag
, &str
[i
]);
2050 /* convert a long to a string (unterminated). returns the number of
2051 characters in the result.
2053 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2055 scm_iint2str (long num
, int rad
, char *p
)
2059 unsigned long n
= (num
< 0) ? -num
: num
;
2061 for (n
/= rad
; n
> 0; n
/= rad
)
2078 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2083 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2085 "Return a string holding the external representation of the\n"
2086 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2087 "inexact, a radix of 10 will be used.")
2088 #define FUNC_NAME s_scm_number_to_string
2092 if (SCM_UNBNDP (radix
))
2096 SCM_VALIDATE_INUM (2, radix
);
2097 base
= SCM_INUM (radix
);
2098 /* FIXME: ask if range limit was OK, and if so, document */
2099 SCM_ASSERT_RANGE (2, radix
, (base
>= 2) && (base
<= 36));
2104 char num_buf
[SCM_INTBUFLEN
];
2105 size_t length
= scm_iint2str (SCM_INUM (n
), base
, num_buf
);
2106 return scm_mem2string (num_buf
, length
);
2108 else if (SCM_BIGP (n
))
2110 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2111 scm_remember_upto_here_1 (n
);
2112 return scm_take0str (str
);
2114 else if (SCM_FRACTIONP (n
))
2116 scm_i_fraction_reduce (n
);
2117 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2118 scm_mem2string ("/", 1),
2119 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2121 else if (SCM_INEXACTP (n
))
2123 char num_buf
[FLOBUFLEN
];
2124 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
));
2127 SCM_WRONG_TYPE_ARG (1, n
);
2132 /* These print routines used to be stubbed here so that scm_repl.c
2133 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2136 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2138 char num_buf
[FLOBUFLEN
];
2139 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2144 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2147 char num_buf
[FLOBUFLEN
];
2148 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2153 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2156 scm_i_fraction_reduce (sexp
);
2157 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2158 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2159 scm_remember_upto_here_1 (str
);
2164 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2166 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2167 scm_remember_upto_here_1 (exp
);
2168 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2172 /*** END nums->strs ***/
2175 /*** STRINGS -> NUMBERS ***/
2177 /* The following functions implement the conversion from strings to numbers.
2178 * The implementation somehow follows the grammar for numbers as it is given
2179 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2180 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2181 * points should be noted about the implementation:
2182 * * Each function keeps a local index variable 'idx' that points at the
2183 * current position within the parsed string. The global index is only
2184 * updated if the function could parse the corresponding syntactic unit
2186 * * Similarly, the functions keep track of indicators of inexactness ('#',
2187 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2188 * global exactness information is only updated after each part has been
2189 * successfully parsed.
2190 * * Sequences of digits are parsed into temporary variables holding fixnums.
2191 * Only if these fixnums would overflow, the result variables are updated
2192 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2193 * the temporary variables holding the fixnums are cleared, and the process
2194 * starts over again. If for example fixnums were able to store five decimal
2195 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2196 * and the result was computed as 12345 * 100000 + 67890. In other words,
2197 * only every five digits two bignum operations were performed.
2200 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2202 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2204 /* In non ASCII-style encodings the following macro might not work. */
2205 #define XDIGIT2UINT(d) (isdigit (d) ? (d) - '0' : tolower (d) - 'a' + 10)
2208 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2209 unsigned int radix
, enum t_exactness
*p_exactness
)
2211 unsigned int idx
= *p_idx
;
2212 unsigned int hash_seen
= 0;
2213 scm_t_bits shift
= 1;
2215 unsigned int digit_value
;
2225 digit_value
= XDIGIT2UINT (c
);
2226 if (digit_value
>= radix
)
2230 result
= SCM_MAKINUM (digit_value
);
2238 digit_value
= XDIGIT2UINT (c
);
2239 if (digit_value
>= radix
)
2251 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2253 result
= scm_product (result
, SCM_MAKINUM (shift
));
2255 result
= scm_sum (result
, SCM_MAKINUM (add
));
2262 shift
= shift
* radix
;
2263 add
= add
* radix
+ digit_value
;
2268 result
= scm_product (result
, SCM_MAKINUM (shift
));
2270 result
= scm_sum (result
, SCM_MAKINUM (add
));
2274 *p_exactness
= INEXACT
;
2280 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2281 * covers the parts of the rules that start at a potential point. The value
2282 * of the digits up to the point have been parsed by the caller and are given
2283 * in variable result. The content of *p_exactness indicates, whether a hash
2284 * has already been seen in the digits before the point.
2287 /* In non ASCII-style encodings the following macro might not work. */
2288 #define DIGIT2UINT(d) ((d) - '0')
2291 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2292 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2294 unsigned int idx
= *p_idx
;
2295 enum t_exactness x
= *p_exactness
;
2300 if (mem
[idx
] == '.')
2302 scm_t_bits shift
= 1;
2304 unsigned int digit_value
;
2305 SCM big_shift
= SCM_MAKINUM (1);
2316 digit_value
= DIGIT2UINT (c
);
2327 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2329 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2330 result
= scm_product (result
, SCM_MAKINUM (shift
));
2332 result
= scm_sum (result
, SCM_MAKINUM (add
));
2340 add
= add
* 10 + digit_value
;
2346 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2347 result
= scm_product (result
, SCM_MAKINUM (shift
));
2348 result
= scm_sum (result
, SCM_MAKINUM (add
));
2351 result
= scm_divide (result
, big_shift
);
2353 /* We've seen a decimal point, thus the value is implicitly inexact. */
2365 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2396 exponent
= DIGIT2UINT (c
);
2403 if (exponent
<= SCM_MAXEXP
)
2404 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2410 if (exponent
> SCM_MAXEXP
)
2412 size_t exp_len
= idx
- start
;
2413 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2414 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2415 scm_out_of_range ("string->number", exp_num
);
2418 e
= scm_integer_expt (SCM_MAKINUM (10), SCM_MAKINUM (exponent
));
2420 result
= scm_product (result
, e
);
2422 result
= scm_divide2real (result
, e
);
2424 /* We've seen an exponent, thus the value is implicitly inexact. */
2442 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2445 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2446 unsigned int radix
, enum t_exactness
*p_exactness
)
2448 unsigned int idx
= *p_idx
;
2454 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2460 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2462 enum t_exactness x
= EXACT
;
2464 /* Cobble up the fractional part. We might want to set the
2465 NaN's mantissa from it. */
2467 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2472 if (mem
[idx
] == '.')
2476 else if (idx
+ 1 == len
)
2478 else if (!isdigit (mem
[idx
+ 1]))
2481 result
= mem2decimal_from_point (SCM_MAKINUM (0), mem
, len
,
2482 p_idx
, p_exactness
);
2486 enum t_exactness x
= EXACT
;
2489 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2490 if (SCM_FALSEP (uinteger
))
2495 else if (mem
[idx
] == '/')
2501 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2502 if (SCM_FALSEP (divisor
))
2505 /* both are int/big here, I assume */
2506 result
= scm_make_ratio (uinteger
, divisor
);
2508 else if (radix
== 10)
2510 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2511 if (SCM_FALSEP (result
))
2522 /* When returning an inexact zero, make sure it is represented as a
2523 floating point value so that we can change its sign.
2525 if (SCM_EQ_P (result
, SCM_MAKINUM(0)) && *p_exactness
== INEXACT
)
2526 result
= scm_make_real (0.0);
2532 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2535 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2536 unsigned int radix
, enum t_exactness
*p_exactness
)
2560 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2561 if (SCM_FALSEP (ureal
))
2563 /* input must be either +i or -i */
2568 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2574 return scm_make_rectangular (SCM_MAKINUM (0), SCM_MAKINUM (sign
));
2581 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2582 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2591 /* either +<ureal>i or -<ureal>i */
2598 return scm_make_rectangular (SCM_MAKINUM (0), ureal
);
2601 /* polar input: <real>@<real>. */
2626 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2627 if (SCM_FALSEP (angle
))
2632 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2633 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2635 result
= scm_make_polar (ureal
, angle
);
2640 /* expecting input matching <real>[+-]<ureal>?i */
2647 int sign
= (c
== '+') ? 1 : -1;
2648 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2650 if (SCM_FALSEP (imag
))
2651 imag
= SCM_MAKINUM (sign
);
2652 else if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2653 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2657 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2664 return scm_make_rectangular (ureal
, imag
);
2673 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2675 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2678 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2680 unsigned int idx
= 0;
2681 unsigned int radix
= NO_RADIX
;
2682 enum t_exactness forced_x
= NO_EXACTNESS
;
2683 enum t_exactness implicit_x
= EXACT
;
2686 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2687 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2689 switch (mem
[idx
+ 1])
2692 if (radix
!= NO_RADIX
)
2697 if (radix
!= NO_RADIX
)
2702 if (forced_x
!= NO_EXACTNESS
)
2707 if (forced_x
!= NO_EXACTNESS
)
2712 if (radix
!= NO_RADIX
)
2717 if (radix
!= NO_RADIX
)
2727 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2728 if (radix
== NO_RADIX
)
2729 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2731 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2733 if (SCM_FALSEP (result
))
2739 if (SCM_INEXACTP (result
))
2740 return scm_inexact_to_exact (result
);
2744 if (SCM_INEXACTP (result
))
2747 return scm_exact_to_inexact (result
);
2750 if (implicit_x
== INEXACT
)
2752 if (SCM_INEXACTP (result
))
2755 return scm_exact_to_inexact (result
);
2763 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2764 (SCM string
, SCM radix
),
2765 "Return a number of the maximally precise representation\n"
2766 "expressed by the given @var{string}. @var{radix} must be an\n"
2767 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2768 "is a default radix that may be overridden by an explicit radix\n"
2769 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2770 "supplied, then the default radix is 10. If string is not a\n"
2771 "syntactically valid notation for a number, then\n"
2772 "@code{string->number} returns @code{#f}.")
2773 #define FUNC_NAME s_scm_string_to_number
2777 SCM_VALIDATE_STRING (1, string
);
2778 SCM_VALIDATE_INUM_MIN_DEF_COPY (2, radix
,2,10, base
);
2779 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2780 SCM_STRING_LENGTH (string
),
2782 return scm_return_first (answer
, string
);
2787 /*** END strs->nums ***/
2791 scm_make_real (double x
)
2793 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2795 SCM_REAL_VALUE (z
) = x
;
2801 scm_make_complex (double x
, double y
)
2804 return scm_make_real (x
);
2808 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2810 SCM_COMPLEX_REAL (z
) = x
;
2811 SCM_COMPLEX_IMAG (z
) = y
;
2818 scm_bigequal (SCM x
, SCM y
)
2820 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2821 scm_remember_upto_here_2 (x
, y
);
2822 return SCM_BOOL (0 == result
);
2826 scm_real_equalp (SCM x
, SCM y
)
2828 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2832 scm_complex_equalp (SCM x
, SCM y
)
2834 return SCM_BOOL (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2835 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
2839 scm_i_fraction_equalp (SCM x
, SCM y
)
2841 scm_i_fraction_reduce (x
);
2842 scm_i_fraction_reduce (y
);
2843 if (SCM_FALSEP (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
2844 SCM_FRACTION_NUMERATOR (y
)))
2845 || SCM_FALSEP (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
2846 SCM_FRACTION_DENOMINATOR (y
))))
2853 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
2854 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
2855 * "else. Note that the sets of complex, real, rational and\n"
2856 * "integer values form subsets of the set of numbers, i. e. the\n"
2857 * "predicate will be fulfilled for any number."
2859 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
2861 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
2862 "otherwise. Note that the sets of real, rational and integer\n"
2863 "values form subsets of the set of complex numbers, i. e. the\n"
2864 "predicate will also be fulfilled if @var{x} is a real,\n"
2865 "rational or integer number.")
2866 #define FUNC_NAME s_scm_number_p
2868 return SCM_BOOL (SCM_NUMBERP (x
));
2873 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
2875 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
2876 "otherwise. Note that the set of integer values forms a subset of\n"
2877 "the set of real numbers, i. e. the predicate will also be\n"
2878 "fulfilled if @var{x} is an integer number.")
2879 #define FUNC_NAME s_scm_real_p
2881 /* we can't represent irrational numbers. */
2882 return scm_rational_p (x
);
2886 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
2888 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
2889 "otherwise. Note that the set of integer values forms a subset of\n"
2890 "the set of rational numbers, i. e. the predicate will also be\n"
2891 "fulfilled if @var{x} is an integer number.")
2892 #define FUNC_NAME s_scm_rational_p
2896 else if (SCM_IMP (x
))
2898 else if (SCM_BIGP (x
))
2900 else if (SCM_FRACTIONP (x
))
2902 else if (SCM_REALP (x
))
2903 /* due to their limited precision, all floating point numbers are
2904 rational as well. */
2912 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
2914 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
2916 #define FUNC_NAME s_scm_integer_p
2925 if (!SCM_INEXACTP (x
))
2927 if (SCM_COMPLEXP (x
))
2929 r
= SCM_REAL_VALUE (x
);
2937 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
2939 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
2941 #define FUNC_NAME s_scm_inexact_p
2943 if (SCM_INEXACTP (x
))
2945 if (SCM_NUMBERP (x
))
2947 SCM_WRONG_TYPE_ARG (1, x
);
2952 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
2953 /* "Return @code{#t} if all parameters are numerically equal." */
2955 scm_num_eq_p (SCM x
, SCM y
)
2960 long xx
= SCM_INUM (x
);
2963 long yy
= SCM_INUM (y
);
2964 return SCM_BOOL (xx
== yy
);
2966 else if (SCM_BIGP (y
))
2968 else if (SCM_REALP (y
))
2969 return SCM_BOOL ((double) xx
== SCM_REAL_VALUE (y
));
2970 else if (SCM_COMPLEXP (y
))
2971 return SCM_BOOL (((double) xx
== SCM_COMPLEX_REAL (y
))
2972 && (0.0 == SCM_COMPLEX_IMAG (y
)));
2973 else if (SCM_FRACTIONP (y
))
2976 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2978 else if (SCM_BIGP (x
))
2982 else if (SCM_BIGP (y
))
2984 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2985 scm_remember_upto_here_2 (x
, y
);
2986 return SCM_BOOL (0 == cmp
);
2988 else if (SCM_REALP (y
))
2991 if (xisnan (SCM_REAL_VALUE (y
)))
2993 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
2994 scm_remember_upto_here_1 (x
);
2995 return SCM_BOOL (0 == cmp
);
2997 else if (SCM_COMPLEXP (y
))
3000 if (0.0 != SCM_COMPLEX_IMAG (y
))
3002 if (xisnan (SCM_COMPLEX_REAL (y
)))
3004 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3005 scm_remember_upto_here_1 (x
);
3006 return SCM_BOOL (0 == cmp
);
3008 else if (SCM_FRACTIONP (y
))
3011 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3013 else if (SCM_REALP (x
))
3016 return SCM_BOOL (SCM_REAL_VALUE (x
) == (double) SCM_INUM (y
));
3017 else if (SCM_BIGP (y
))
3020 if (xisnan (SCM_REAL_VALUE (x
)))
3022 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3023 scm_remember_upto_here_1 (y
);
3024 return SCM_BOOL (0 == cmp
);
3026 else if (SCM_REALP (y
))
3027 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3028 else if (SCM_COMPLEXP (y
))
3029 return SCM_BOOL ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3030 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3031 else if (SCM_FRACTIONP (y
))
3033 double xx
= SCM_REAL_VALUE (x
);
3037 return SCM_BOOL (xx
< 0.0);
3038 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3042 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3044 else if (SCM_COMPLEXP (x
))
3047 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == (double) SCM_INUM (y
))
3048 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3049 else if (SCM_BIGP (y
))
3052 if (0.0 != SCM_COMPLEX_IMAG (x
))
3054 if (xisnan (SCM_COMPLEX_REAL (x
)))
3056 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3057 scm_remember_upto_here_1 (y
);
3058 return SCM_BOOL (0 == cmp
);
3060 else if (SCM_REALP (y
))
3061 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3062 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3063 else if (SCM_COMPLEXP (y
))
3064 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3065 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3066 else if (SCM_FRACTIONP (y
))
3069 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3071 xx
= SCM_COMPLEX_REAL (x
);
3075 return SCM_BOOL (xx
< 0.0);
3076 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3080 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3082 else if (SCM_FRACTIONP (x
))
3086 else if (SCM_BIGP (y
))
3088 else if (SCM_REALP (y
))
3090 double yy
= SCM_REAL_VALUE (y
);
3094 return SCM_BOOL (0.0 < yy
);
3095 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3098 else if (SCM_COMPLEXP (y
))
3101 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3103 yy
= SCM_COMPLEX_REAL (y
);
3107 return SCM_BOOL (0.0 < yy
);
3108 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3111 else if (SCM_FRACTIONP (y
))
3112 return scm_i_fraction_equalp (x
, y
);
3114 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3117 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3121 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3122 done are good for inums, but for bignums an answer can almost always be
3123 had by just examining a few high bits of the operands, as done by GMP in
3124 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3125 of the float exponent to take into account. */
3127 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3128 /* "Return @code{#t} if the list of parameters is monotonically\n"
3132 scm_less_p (SCM x
, SCM y
)
3137 long xx
= SCM_INUM (x
);
3140 long yy
= SCM_INUM (y
);
3141 return SCM_BOOL (xx
< yy
);
3143 else if (SCM_BIGP (y
))
3145 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3146 scm_remember_upto_here_1 (y
);
3147 return SCM_BOOL (sgn
> 0);
3149 else if (SCM_REALP (y
))
3150 return SCM_BOOL ((double) xx
< SCM_REAL_VALUE (y
));
3151 else if (SCM_FRACTIONP (y
))
3153 /* "x < a/b" becomes "x*b < a" */
3155 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3156 y
= SCM_FRACTION_NUMERATOR (y
);
3160 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3162 else if (SCM_BIGP (x
))
3166 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3167 scm_remember_upto_here_1 (x
);
3168 return SCM_BOOL (sgn
< 0);
3170 else if (SCM_BIGP (y
))
3172 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3173 scm_remember_upto_here_2 (x
, y
);
3174 return SCM_BOOL (cmp
< 0);
3176 else if (SCM_REALP (y
))
3179 if (xisnan (SCM_REAL_VALUE (y
)))
3181 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3182 scm_remember_upto_here_1 (x
);
3183 return SCM_BOOL (cmp
< 0);
3185 else if (SCM_FRACTIONP (y
))
3188 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3190 else if (SCM_REALP (x
))
3193 return SCM_BOOL (SCM_REAL_VALUE (x
) < (double) SCM_INUM (y
));
3194 else if (SCM_BIGP (y
))
3197 if (xisnan (SCM_REAL_VALUE (x
)))
3199 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3200 scm_remember_upto_here_1 (y
);
3201 return SCM_BOOL (cmp
> 0);
3203 else if (SCM_REALP (y
))
3204 return SCM_BOOL (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3205 else if (SCM_FRACTIONP (y
))
3207 double xx
= SCM_REAL_VALUE (x
);
3211 return SCM_BOOL (xx
< 0.0);
3212 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3216 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3218 else if (SCM_FRACTIONP (x
))
3220 if (SCM_INUMP (y
) || SCM_BIGP (y
))
3222 /* "a/b < y" becomes "a < y*b" */
3223 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3224 x
= SCM_FRACTION_NUMERATOR (x
);
3227 else if (SCM_REALP (y
))
3229 double yy
= SCM_REAL_VALUE (y
);
3233 return SCM_BOOL (0.0 < yy
);
3234 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3237 else if (SCM_FRACTIONP (y
))
3239 /* "a/b < c/d" becomes "a*d < c*b" */
3240 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3241 SCM_FRACTION_DENOMINATOR (y
));
3242 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3243 SCM_FRACTION_DENOMINATOR (x
));
3249 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3252 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3256 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3257 /* "Return @code{#t} if the list of parameters is monotonically\n"
3260 #define FUNC_NAME s_scm_gr_p
3262 scm_gr_p (SCM x
, SCM y
)
3264 if (!SCM_NUMBERP (x
))
3265 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3266 else if (!SCM_NUMBERP (y
))
3267 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3269 return scm_less_p (y
, x
);
3274 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3275 /* "Return @code{#t} if the list of parameters is monotonically\n"
3278 #define FUNC_NAME s_scm_leq_p
3280 scm_leq_p (SCM x
, SCM y
)
3282 if (!SCM_NUMBERP (x
))
3283 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3284 else if (!SCM_NUMBERP (y
))
3285 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3286 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3289 return SCM_BOOL_NOT (scm_less_p (y
, x
));
3294 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3295 /* "Return @code{#t} if the list of parameters is monotonically\n"
3298 #define FUNC_NAME s_scm_geq_p
3300 scm_geq_p (SCM x
, SCM y
)
3302 if (!SCM_NUMBERP (x
))
3303 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3304 else if (!SCM_NUMBERP (y
))
3305 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3306 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3309 return SCM_BOOL_NOT (scm_less_p (x
, y
));
3314 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3315 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3322 return SCM_BOOL (SCM_EQ_P (z
, SCM_INUM0
));
3323 else if (SCM_BIGP (z
))
3325 else if (SCM_REALP (z
))
3326 return SCM_BOOL (SCM_REAL_VALUE (z
) == 0.0);
3327 else if (SCM_COMPLEXP (z
))
3328 return SCM_BOOL (SCM_COMPLEX_REAL (z
) == 0.0
3329 && SCM_COMPLEX_IMAG (z
) == 0.0);
3330 else if (SCM_FRACTIONP (z
))
3333 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3337 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3338 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3342 scm_positive_p (SCM x
)
3345 return SCM_BOOL (SCM_INUM (x
) > 0);
3346 else if (SCM_BIGP (x
))
3348 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3349 scm_remember_upto_here_1 (x
);
3350 return SCM_BOOL (sgn
> 0);
3352 else if (SCM_REALP (x
))
3353 return SCM_BOOL(SCM_REAL_VALUE (x
) > 0.0);
3354 else if (SCM_FRACTIONP (x
))
3355 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3357 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3361 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3362 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3366 scm_negative_p (SCM x
)
3369 return SCM_BOOL (SCM_INUM (x
) < 0);
3370 else if (SCM_BIGP (x
))
3372 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3373 scm_remember_upto_here_1 (x
);
3374 return SCM_BOOL (sgn
< 0);
3376 else if (SCM_REALP (x
))
3377 return SCM_BOOL(SCM_REAL_VALUE (x
) < 0.0);
3378 else if (SCM_FRACTIONP (x
))
3379 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3381 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3385 /* scm_min and scm_max return an inexact when either argument is inexact, as
3386 required by r5rs. On that basis, for exact/inexact combinations the
3387 exact is converted to inexact to compare and possibly return. This is
3388 unlike scm_less_p above which takes some trouble to preserve all bits in
3389 its test, such trouble is not required for min and max. */
3391 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3392 /* "Return the maximum of all parameter values."
3395 scm_max (SCM x
, SCM y
)
3400 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3401 else if (SCM_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3404 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3409 long xx
= SCM_INUM (x
);
3412 long yy
= SCM_INUM (y
);
3413 return (xx
< yy
) ? y
: x
;
3415 else if (SCM_BIGP (y
))
3417 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3418 scm_remember_upto_here_1 (y
);
3419 return (sgn
< 0) ? x
: y
;
3421 else if (SCM_REALP (y
))
3424 /* if y==NaN then ">" is false and we return NaN */
3425 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3427 else if (SCM_FRACTIONP (y
))
3430 return (z
> scm_i_fraction2double (y
)) ? x
: y
;
3433 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3435 else if (SCM_BIGP (x
))
3439 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3440 scm_remember_upto_here_1 (x
);
3441 return (sgn
< 0) ? y
: x
;
3443 else if (SCM_BIGP (y
))
3445 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3446 scm_remember_upto_here_2 (x
, y
);
3447 return (cmp
> 0) ? x
: y
;
3449 else if (SCM_REALP (y
))
3451 /* if y==NaN then xx>yy is false, so we return the NaN y */
3454 xx
= scm_i_big2dbl (x
);
3455 yy
= SCM_REAL_VALUE (y
);
3456 return (xx
> yy
? scm_make_real (xx
) : y
);
3458 else if (SCM_FRACTIONP (y
))
3460 double yy
= scm_i_fraction2double (y
);
3462 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3463 scm_remember_upto_here_1 (x
);
3464 return (cmp
> 0) ? x
: y
;
3467 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3469 else if (SCM_REALP (x
))
3473 double z
= SCM_INUM (y
);
3474 /* if x==NaN then "<" is false and we return NaN */
3475 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3477 else if (SCM_BIGP (y
))
3479 SCM t
= x
; x
= y
; y
= t
;
3482 else if (SCM_REALP (y
))
3484 /* if x==NaN then our explicit check means we return NaN
3485 if y==NaN then ">" is false and we return NaN
3486 calling isnan is unavoidable, since it's the only way to know
3487 which of x or y causes any compares to be false */
3488 double xx
= SCM_REAL_VALUE (x
);
3489 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3491 else if (SCM_FRACTIONP (y
))
3493 double yy
= scm_i_fraction2double (y
);
3494 double xx
= SCM_REAL_VALUE (x
);
3495 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3498 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3500 else if (SCM_FRACTIONP (x
))
3504 double z
= SCM_INUM (y
);
3505 return (scm_i_fraction2double (x
) < z
) ? y
: x
;
3507 else if (SCM_BIGP (y
))
3509 double xx
= scm_i_fraction2double (x
);
3511 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3512 scm_remember_upto_here_1 (y
);
3513 return (cmp
< 0) ? x
: y
;
3515 else if (SCM_REALP (y
))
3517 double xx
= scm_i_fraction2double (x
);
3518 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3520 else if (SCM_FRACTIONP (y
))
3522 double yy
= scm_i_fraction2double (y
);
3523 double xx
= scm_i_fraction2double (x
);
3524 return (xx
< yy
) ? y
: x
;
3527 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3530 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3534 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3535 /* "Return the minium of all parameter values."
3538 scm_min (SCM x
, SCM y
)
3543 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3544 else if (SCM_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3547 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3552 long xx
= SCM_INUM (x
);
3555 long yy
= SCM_INUM (y
);
3556 return (xx
< yy
) ? x
: y
;
3558 else if (SCM_BIGP (y
))
3560 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3561 scm_remember_upto_here_1 (y
);
3562 return (sgn
< 0) ? y
: x
;
3564 else if (SCM_REALP (y
))
3567 /* if y==NaN then "<" is false and we return NaN */
3568 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3570 else if (SCM_FRACTIONP (y
))
3573 return (z
< scm_i_fraction2double (y
)) ? x
: y
;
3576 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3578 else if (SCM_BIGP (x
))
3582 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3583 scm_remember_upto_here_1 (x
);
3584 return (sgn
< 0) ? x
: y
;
3586 else if (SCM_BIGP (y
))
3588 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3589 scm_remember_upto_here_2 (x
, y
);
3590 return (cmp
> 0) ? y
: x
;
3592 else if (SCM_REALP (y
))
3594 /* if y==NaN then xx<yy is false, so we return the NaN y */
3597 xx
= scm_i_big2dbl (x
);
3598 yy
= SCM_REAL_VALUE (y
);
3599 return (xx
< yy
? scm_make_real (xx
) : y
);
3601 else if (SCM_FRACTIONP (y
))
3603 double yy
= scm_i_fraction2double (y
);
3605 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3606 scm_remember_upto_here_1 (x
);
3607 return (cmp
> 0) ? y
: x
;
3610 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3612 else if (SCM_REALP (x
))
3616 double z
= SCM_INUM (y
);
3617 /* if x==NaN then "<" is false and we return NaN */
3618 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3620 else if (SCM_BIGP (y
))
3622 SCM t
= x
; x
= y
; y
= t
;
3625 else if (SCM_REALP (y
))
3627 /* if x==NaN then our explicit check means we return NaN
3628 if y==NaN then "<" is false and we return NaN
3629 calling isnan is unavoidable, since it's the only way to know
3630 which of x or y causes any compares to be false */
3631 double xx
= SCM_REAL_VALUE (x
);
3632 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3634 else if (SCM_FRACTIONP (y
))
3636 double yy
= scm_i_fraction2double (y
);
3637 double xx
= SCM_REAL_VALUE (x
);
3638 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3641 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3643 else if (SCM_FRACTIONP (x
))
3647 double z
= SCM_INUM (y
);
3648 return (scm_i_fraction2double (x
) < z
) ? x
: y
;
3650 else if (SCM_BIGP (y
))
3652 double xx
= scm_i_fraction2double (x
);
3654 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3655 scm_remember_upto_here_1 (y
);
3656 return (cmp
< 0) ? y
: x
;
3658 else if (SCM_REALP (y
))
3660 double xx
= scm_i_fraction2double (x
);
3661 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3663 else if (SCM_FRACTIONP (y
))
3665 double yy
= scm_i_fraction2double (y
);
3666 double xx
= scm_i_fraction2double (x
);
3667 return (xx
< yy
) ? x
: y
;
3670 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3673 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3677 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3678 /* "Return the sum of all parameter values. Return 0 if called without\n"
3682 scm_sum (SCM x
, SCM y
)
3686 if (SCM_NUMBERP (x
)) return x
;
3687 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3688 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3695 long xx
= SCM_INUM (x
);
3696 long yy
= SCM_INUM (y
);
3697 long int z
= xx
+ yy
;
3698 return SCM_FIXABLE (z
) ? SCM_MAKINUM (z
) : scm_i_long2big (z
);
3700 else if (SCM_BIGP (y
))
3705 else if (SCM_REALP (y
))
3707 long int xx
= SCM_INUM (x
);
3708 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3710 else if (SCM_COMPLEXP (y
))
3712 long int xx
= SCM_INUM (x
);
3713 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3714 SCM_COMPLEX_IMAG (y
));
3716 else if (SCM_FRACTIONP (y
))
3717 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3718 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3719 SCM_FRACTION_DENOMINATOR (y
));
3721 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3722 } else if (SCM_BIGP (x
))
3729 inum
= SCM_INUM (y
);
3732 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3735 SCM result
= scm_i_mkbig ();
3736 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3737 scm_remember_upto_here_1 (x
);
3738 /* we know the result will have to be a bignum */
3741 return scm_i_normbig (result
);
3745 SCM result
= scm_i_mkbig ();
3746 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3747 scm_remember_upto_here_1 (x
);
3748 /* we know the result will have to be a bignum */
3751 return scm_i_normbig (result
);
3754 else if (SCM_BIGP (y
))
3756 SCM result
= scm_i_mkbig ();
3757 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3758 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3759 mpz_add (SCM_I_BIG_MPZ (result
),
3762 scm_remember_upto_here_2 (x
, y
);
3763 /* we know the result will have to be a bignum */
3766 return scm_i_normbig (result
);
3768 else if (SCM_REALP (y
))
3770 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3771 scm_remember_upto_here_1 (x
);
3772 return scm_make_real (result
);
3774 else if (SCM_COMPLEXP (y
))
3776 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3777 + SCM_COMPLEX_REAL (y
));
3778 scm_remember_upto_here_1 (x
);
3779 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3781 else if (SCM_FRACTIONP (y
))
3782 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3783 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3784 SCM_FRACTION_DENOMINATOR (y
));
3786 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3788 else if (SCM_REALP (x
))
3791 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_INUM (y
));
3792 else if (SCM_BIGP (y
))
3794 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3795 scm_remember_upto_here_1 (y
);
3796 return scm_make_real (result
);
3798 else if (SCM_REALP (y
))
3799 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3800 else if (SCM_COMPLEXP (y
))
3801 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3802 SCM_COMPLEX_IMAG (y
));
3803 else if (SCM_FRACTIONP (y
))
3804 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3806 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3808 else if (SCM_COMPLEXP (x
))
3811 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_INUM (y
),
3812 SCM_COMPLEX_IMAG (x
));
3813 else if (SCM_BIGP (y
))
3815 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3816 + SCM_COMPLEX_REAL (x
));
3817 scm_remember_upto_here_1 (y
);
3818 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3820 else if (SCM_REALP (y
))
3821 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3822 SCM_COMPLEX_IMAG (x
));
3823 else if (SCM_COMPLEXP (y
))
3824 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3825 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3826 else if (SCM_FRACTIONP (y
))
3827 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3828 SCM_COMPLEX_IMAG (x
));
3830 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3832 else if (SCM_FRACTIONP (x
))
3835 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3836 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3837 SCM_FRACTION_DENOMINATOR (x
));
3838 else if (SCM_BIGP (y
))
3839 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3840 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3841 SCM_FRACTION_DENOMINATOR (x
));
3842 else if (SCM_REALP (y
))
3843 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3844 else if (SCM_COMPLEXP (y
))
3845 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3846 SCM_COMPLEX_IMAG (y
));
3847 else if (SCM_FRACTIONP (y
))
3848 /* a/b + c/d = (ad + bc) / bd */
3849 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3850 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3851 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3853 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3856 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
3860 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
3861 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
3862 * the sum of all but the first argument are subtracted from the first
3864 #define FUNC_NAME s_difference
3866 scm_difference (SCM x
, SCM y
)
3871 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
3875 long xx
= -SCM_INUM (x
);
3876 if (SCM_FIXABLE (xx
))
3877 return SCM_MAKINUM (xx
);
3879 return scm_i_long2big (xx
);
3881 else if (SCM_BIGP (x
))
3882 /* FIXME: do we really need to normalize here? */
3883 return scm_i_normbig (scm_i_clonebig (x
, 0));
3884 else if (SCM_REALP (x
))
3885 return scm_make_real (-SCM_REAL_VALUE (x
));
3886 else if (SCM_COMPLEXP (x
))
3887 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
3888 -SCM_COMPLEX_IMAG (x
));
3889 else if (SCM_FRACTIONP (x
))
3890 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
3891 SCM_FRACTION_DENOMINATOR (x
));
3893 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
3900 long int xx
= SCM_INUM (x
);
3901 long int yy
= SCM_INUM (y
);
3902 long int z
= xx
- yy
;
3903 if (SCM_FIXABLE (z
))
3904 return SCM_MAKINUM (z
);
3906 return scm_i_long2big (z
);
3908 else if (SCM_BIGP (y
))
3910 /* inum-x - big-y */
3911 long xx
= SCM_INUM (x
);
3914 return scm_i_clonebig (y
, 0);
3917 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3918 SCM result
= scm_i_mkbig ();
3921 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
3924 /* x - y == -(y + -x) */
3925 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
3926 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
3928 scm_remember_upto_here_1 (y
);
3930 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
3931 /* we know the result will have to be a bignum */
3934 return scm_i_normbig (result
);
3937 else if (SCM_REALP (y
))
3939 long int xx
= SCM_INUM (x
);
3940 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
3942 else if (SCM_COMPLEXP (y
))
3944 long int xx
= SCM_INUM (x
);
3945 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
3946 - SCM_COMPLEX_IMAG (y
));
3948 else if (SCM_FRACTIONP (y
))
3949 /* a - b/c = (ac - b) / c */
3950 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3951 SCM_FRACTION_NUMERATOR (y
)),
3952 SCM_FRACTION_DENOMINATOR (y
));
3954 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3956 else if (SCM_BIGP (x
))
3960 /* big-x - inum-y */
3961 long yy
= SCM_INUM (y
);
3962 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3964 scm_remember_upto_here_1 (x
);
3966 return SCM_FIXABLE (-yy
) ? SCM_MAKINUM (-yy
) : scm_long2num (-yy
);
3969 SCM result
= scm_i_mkbig ();
3972 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
3974 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
3975 scm_remember_upto_here_1 (x
);
3977 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
3978 /* we know the result will have to be a bignum */
3981 return scm_i_normbig (result
);
3984 else if (SCM_BIGP (y
))
3986 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3987 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3988 SCM result
= scm_i_mkbig ();
3989 mpz_sub (SCM_I_BIG_MPZ (result
),
3992 scm_remember_upto_here_2 (x
, y
);
3993 /* we know the result will have to be a bignum */
3994 if ((sgn_x
== 1) && (sgn_y
== -1))
3996 if ((sgn_x
== -1) && (sgn_y
== 1))
3998 return scm_i_normbig (result
);
4000 else if (SCM_REALP (y
))
4002 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4003 scm_remember_upto_here_1 (x
);
4004 return scm_make_real (result
);
4006 else if (SCM_COMPLEXP (y
))
4008 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4009 - SCM_COMPLEX_REAL (y
));
4010 scm_remember_upto_here_1 (x
);
4011 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
4013 else if (SCM_FRACTIONP (y
))
4014 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4015 SCM_FRACTION_NUMERATOR (y
)),
4016 SCM_FRACTION_DENOMINATOR (y
));
4017 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4019 else if (SCM_REALP (x
))
4022 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_INUM (y
));
4023 else if (SCM_BIGP (y
))
4025 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4026 scm_remember_upto_here_1 (x
);
4027 return scm_make_real (result
);
4029 else if (SCM_REALP (y
))
4030 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4031 else if (SCM_COMPLEXP (y
))
4032 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4033 -SCM_COMPLEX_IMAG (y
));
4034 else if (SCM_FRACTIONP (y
))
4035 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4037 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4039 else if (SCM_COMPLEXP (x
))
4042 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_INUM (y
),
4043 SCM_COMPLEX_IMAG (x
));
4044 else if (SCM_BIGP (y
))
4046 double real_part
= (SCM_COMPLEX_REAL (x
)
4047 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4048 scm_remember_upto_here_1 (x
);
4049 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
4051 else if (SCM_REALP (y
))
4052 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4053 SCM_COMPLEX_IMAG (x
));
4054 else if (SCM_COMPLEXP (y
))
4055 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4056 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4057 else if (SCM_FRACTIONP (y
))
4058 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4059 SCM_COMPLEX_IMAG (x
));
4061 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4063 else if (SCM_FRACTIONP (x
))
4066 /* a/b - c = (a - cb) / b */
4067 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4068 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4069 SCM_FRACTION_DENOMINATOR (x
));
4070 else if (SCM_BIGP (y
))
4071 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4072 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4073 SCM_FRACTION_DENOMINATOR (x
));
4074 else if (SCM_REALP (y
))
4075 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4076 else if (SCM_COMPLEXP (y
))
4077 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4078 -SCM_COMPLEX_IMAG (y
));
4079 else if (SCM_FRACTIONP (y
))
4080 /* a/b - c/d = (ad - bc) / bd */
4081 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4082 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4083 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4085 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4088 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4093 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4094 /* "Return the product of all arguments. If called without arguments,\n"
4098 scm_product (SCM x
, SCM y
)
4103 return SCM_MAKINUM (1L);
4104 else if (SCM_NUMBERP (x
))
4107 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4119 case 0: return x
; break;
4120 case 1: return y
; break;
4125 long yy
= SCM_INUM (y
);
4127 SCM k
= SCM_MAKINUM (kk
);
4128 if ((kk
== SCM_INUM (k
)) && (kk
/ xx
== yy
))
4132 SCM result
= scm_i_long2big (xx
);
4133 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4134 return scm_i_normbig (result
);
4137 else if (SCM_BIGP (y
))
4139 SCM result
= scm_i_mkbig ();
4140 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4141 scm_remember_upto_here_1 (y
);
4144 else if (SCM_REALP (y
))
4145 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4146 else if (SCM_COMPLEXP (y
))
4147 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4148 xx
* SCM_COMPLEX_IMAG (y
));
4149 else if (SCM_FRACTIONP (y
))
4150 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4151 SCM_FRACTION_DENOMINATOR (y
));
4153 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4155 else if (SCM_BIGP (x
))
4162 else if (SCM_BIGP (y
))
4164 SCM result
= scm_i_mkbig ();
4165 mpz_mul (SCM_I_BIG_MPZ (result
),
4168 scm_remember_upto_here_2 (x
, y
);
4171 else if (SCM_REALP (y
))
4173 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4174 scm_remember_upto_here_1 (x
);
4175 return scm_make_real (result
);
4177 else if (SCM_COMPLEXP (y
))
4179 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4180 scm_remember_upto_here_1 (x
);
4181 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4182 z
* SCM_COMPLEX_IMAG (y
));
4184 else if (SCM_FRACTIONP (y
))
4185 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4186 SCM_FRACTION_DENOMINATOR (y
));
4188 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4190 else if (SCM_REALP (x
))
4193 return scm_make_real (SCM_INUM (y
) * SCM_REAL_VALUE (x
));
4194 else if (SCM_BIGP (y
))
4196 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4197 scm_remember_upto_here_1 (y
);
4198 return scm_make_real (result
);
4200 else if (SCM_REALP (y
))
4201 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4202 else if (SCM_COMPLEXP (y
))
4203 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4204 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4205 else if (SCM_FRACTIONP (y
))
4206 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4208 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4210 else if (SCM_COMPLEXP (x
))
4213 return scm_make_complex (SCM_INUM (y
) * SCM_COMPLEX_REAL (x
),
4214 SCM_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4215 else if (SCM_BIGP (y
))
4217 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4218 scm_remember_upto_here_1 (y
);
4219 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4220 z
* SCM_COMPLEX_IMAG (x
));
4222 else if (SCM_REALP (y
))
4223 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4224 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4225 else if (SCM_COMPLEXP (y
))
4227 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4228 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4229 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4230 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4232 else if (SCM_FRACTIONP (y
))
4234 double yy
= scm_i_fraction2double (y
);
4235 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4236 yy
* SCM_COMPLEX_IMAG (x
));
4239 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4241 else if (SCM_FRACTIONP (x
))
4244 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4245 SCM_FRACTION_DENOMINATOR (x
));
4246 else if (SCM_BIGP (y
))
4247 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4248 SCM_FRACTION_DENOMINATOR (x
));
4249 else if (SCM_REALP (y
))
4250 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4251 else if (SCM_COMPLEXP (y
))
4253 double xx
= scm_i_fraction2double (x
);
4254 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4255 xx
* SCM_COMPLEX_IMAG (y
));
4257 else if (SCM_FRACTIONP (y
))
4258 /* a/b * c/d = ac / bd */
4259 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4260 SCM_FRACTION_NUMERATOR (y
)),
4261 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4262 SCM_FRACTION_DENOMINATOR (y
)));
4264 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4267 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4271 scm_num2dbl (SCM a
, const char *why
)
4272 #define FUNC_NAME why
4275 return (double) SCM_INUM (a
);
4276 else if (SCM_BIGP (a
))
4278 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4279 scm_remember_upto_here_1 (a
);
4282 else if (SCM_REALP (a
))
4283 return (SCM_REAL_VALUE (a
));
4284 else if (SCM_FRACTIONP (a
))
4285 return scm_i_fraction2double (a
);
4287 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4291 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4292 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4293 #define ALLOW_DIVIDE_BY_ZERO
4294 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4297 /* The code below for complex division is adapted from the GNU
4298 libstdc++, which adapted it from f2c's libF77, and is subject to
4301 /****************************************************************
4302 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4304 Permission to use, copy, modify, and distribute this software
4305 and its documentation for any purpose and without fee is hereby
4306 granted, provided that the above copyright notice appear in all
4307 copies and that both that the copyright notice and this
4308 permission notice and warranty disclaimer appear in supporting
4309 documentation, and that the names of AT&T Bell Laboratories or
4310 Bellcore or any of their entities not be used in advertising or
4311 publicity pertaining to distribution of the software without
4312 specific, written prior permission.
4314 AT&T and Bellcore disclaim all warranties with regard to this
4315 software, including all implied warranties of merchantability
4316 and fitness. In no event shall AT&T or Bellcore be liable for
4317 any special, indirect or consequential damages or any damages
4318 whatsoever resulting from loss of use, data or profits, whether
4319 in an action of contract, negligence or other tortious action,
4320 arising out of or in connection with the use or performance of
4322 ****************************************************************/
4324 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4325 /* Divide the first argument by the product of the remaining
4326 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4328 #define FUNC_NAME s_divide
4330 scm_i_divide (SCM x
, SCM y
, int inexact
)
4337 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4338 else if (SCM_INUMP (x
))
4340 long xx
= SCM_INUM (x
);
4341 if (xx
== 1 || xx
== -1)
4343 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4345 scm_num_overflow (s_divide
);
4350 return scm_make_real (1.0 / (double) xx
);
4351 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4354 else if (SCM_BIGP (x
))
4357 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4358 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4360 else if (SCM_REALP (x
))
4362 double xx
= SCM_REAL_VALUE (x
);
4363 #ifndef ALLOW_DIVIDE_BY_ZERO
4365 scm_num_overflow (s_divide
);
4368 return scm_make_real (1.0 / xx
);
4370 else if (SCM_COMPLEXP (x
))
4372 double r
= SCM_COMPLEX_REAL (x
);
4373 double i
= SCM_COMPLEX_IMAG (x
);
4377 double d
= i
* (1.0 + t
* t
);
4378 return scm_make_complex (t
/ d
, -1.0 / d
);
4383 double d
= r
* (1.0 + t
* t
);
4384 return scm_make_complex (1.0 / d
, -t
/ d
);
4387 else if (SCM_FRACTIONP (x
))
4388 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4389 SCM_FRACTION_NUMERATOR (x
));
4391 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4396 long xx
= SCM_INUM (x
);
4399 long yy
= SCM_INUM (y
);
4402 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4403 scm_num_overflow (s_divide
);
4405 return scm_make_real ((double) xx
/ (double) yy
);
4408 else if (xx
% yy
!= 0)
4411 return scm_make_real ((double) xx
/ (double) yy
);
4412 else return scm_make_ratio (x
, y
);
4417 if (SCM_FIXABLE (z
))
4418 return SCM_MAKINUM (z
);
4420 return scm_i_long2big (z
);
4423 else if (SCM_BIGP (y
))
4426 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4427 else return scm_make_ratio (x
, y
);
4429 else if (SCM_REALP (y
))
4431 double yy
= SCM_REAL_VALUE (y
);
4432 #ifndef ALLOW_DIVIDE_BY_ZERO
4434 scm_num_overflow (s_divide
);
4437 return scm_make_real ((double) xx
/ yy
);
4439 else if (SCM_COMPLEXP (y
))
4442 complex_div
: /* y _must_ be a complex number */
4444 double r
= SCM_COMPLEX_REAL (y
);
4445 double i
= SCM_COMPLEX_IMAG (y
);
4449 double d
= i
* (1.0 + t
* t
);
4450 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4455 double d
= r
* (1.0 + t
* t
);
4456 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4460 else if (SCM_FRACTIONP (y
))
4461 /* a / b/c = ac / b */
4462 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4463 SCM_FRACTION_NUMERATOR (y
));
4465 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4467 else if (SCM_BIGP (x
))
4471 long int yy
= SCM_INUM (y
);
4474 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4475 scm_num_overflow (s_divide
);
4477 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4478 scm_remember_upto_here_1 (x
);
4479 return (sgn
== 0) ? scm_nan () : scm_inf ();
4486 /* FIXME: HMM, what are the relative performance issues here?
4487 We need to test. Is it faster on average to test
4488 divisible_p, then perform whichever operation, or is it
4489 faster to perform the integer div opportunistically and
4490 switch to real if there's a remainder? For now we take the
4491 middle ground: test, then if divisible, use the faster div
4494 long abs_yy
= yy
< 0 ? -yy
: yy
;
4495 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4499 SCM result
= scm_i_mkbig ();
4500 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4501 scm_remember_upto_here_1 (x
);
4503 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4504 return scm_i_normbig (result
);
4509 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4510 else return scm_make_ratio (x
, y
);
4514 else if (SCM_BIGP (y
))
4516 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4519 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4520 scm_num_overflow (s_divide
);
4522 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4523 scm_remember_upto_here_1 (x
);
4524 return (sgn
== 0) ? scm_nan () : scm_inf ();
4530 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4534 SCM result
= scm_i_mkbig ();
4535 mpz_divexact (SCM_I_BIG_MPZ (result
),
4538 scm_remember_upto_here_2 (x
, y
);
4539 return scm_i_normbig (result
);
4545 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4546 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4547 scm_remember_upto_here_2 (x
, y
);
4548 return scm_make_real (dbx
/ dby
);
4550 else return scm_make_ratio (x
, y
);
4554 else if (SCM_REALP (y
))
4556 double yy
= SCM_REAL_VALUE (y
);
4557 #ifndef ALLOW_DIVIDE_BY_ZERO
4559 scm_num_overflow (s_divide
);
4562 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4564 else if (SCM_COMPLEXP (y
))
4566 a
= scm_i_big2dbl (x
);
4569 else if (SCM_FRACTIONP (y
))
4570 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4571 SCM_FRACTION_NUMERATOR (y
));
4573 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4575 else if (SCM_REALP (x
))
4577 double rx
= SCM_REAL_VALUE (x
);
4580 long int yy
= SCM_INUM (y
);
4581 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4583 scm_num_overflow (s_divide
);
4586 return scm_make_real (rx
/ (double) yy
);
4588 else if (SCM_BIGP (y
))
4590 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4591 scm_remember_upto_here_1 (y
);
4592 return scm_make_real (rx
/ dby
);
4594 else if (SCM_REALP (y
))
4596 double yy
= SCM_REAL_VALUE (y
);
4597 #ifndef ALLOW_DIVIDE_BY_ZERO
4599 scm_num_overflow (s_divide
);
4602 return scm_make_real (rx
/ yy
);
4604 else if (SCM_COMPLEXP (y
))
4609 else if (SCM_FRACTIONP (y
))
4610 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4612 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4614 else if (SCM_COMPLEXP (x
))
4616 double rx
= SCM_COMPLEX_REAL (x
);
4617 double ix
= SCM_COMPLEX_IMAG (x
);
4620 long int yy
= SCM_INUM (y
);
4621 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4623 scm_num_overflow (s_divide
);
4628 return scm_make_complex (rx
/ d
, ix
/ d
);
4631 else if (SCM_BIGP (y
))
4633 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4634 scm_remember_upto_here_1 (y
);
4635 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4637 else if (SCM_REALP (y
))
4639 double yy
= SCM_REAL_VALUE (y
);
4640 #ifndef ALLOW_DIVIDE_BY_ZERO
4642 scm_num_overflow (s_divide
);
4645 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4647 else if (SCM_COMPLEXP (y
))
4649 double ry
= SCM_COMPLEX_REAL (y
);
4650 double iy
= SCM_COMPLEX_IMAG (y
);
4654 double d
= iy
* (1.0 + t
* t
);
4655 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4660 double d
= ry
* (1.0 + t
* t
);
4661 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4664 else if (SCM_FRACTIONP (y
))
4666 double yy
= scm_i_fraction2double (y
);
4667 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4670 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4672 else if (SCM_FRACTIONP (x
))
4676 long int yy
= SCM_INUM (y
);
4677 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4679 scm_num_overflow (s_divide
);
4682 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4683 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4685 else if (SCM_BIGP (y
))
4687 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4688 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4690 else if (SCM_REALP (y
))
4692 double yy
= SCM_REAL_VALUE (y
);
4693 #ifndef ALLOW_DIVIDE_BY_ZERO
4695 scm_num_overflow (s_divide
);
4698 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4700 else if (SCM_COMPLEXP (y
))
4702 a
= scm_i_fraction2double (x
);
4705 else if (SCM_FRACTIONP (y
))
4706 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4707 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4709 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4712 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4716 scm_divide (SCM x
, SCM y
)
4718 return scm_i_divide (x
, y
, 0);
4721 static SCM
scm_divide2real (SCM x
, SCM y
)
4723 return scm_i_divide (x
, y
, 1);
4729 scm_asinh (double x
)
4734 #define asinh scm_asinh
4735 return log (x
+ sqrt (x
* x
+ 1));
4738 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4739 /* "Return the inverse hyperbolic sine of @var{x}."
4744 scm_acosh (double x
)
4749 #define acosh scm_acosh
4750 return log (x
+ sqrt (x
* x
- 1));
4753 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4754 /* "Return the inverse hyperbolic cosine of @var{x}."
4759 scm_atanh (double x
)
4764 #define atanh scm_atanh
4765 return 0.5 * log ((1 + x
) / (1 - x
));
4768 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4769 /* "Return the inverse hyperbolic tangent of @var{x}."
4773 /* XXX - eventually, we should remove this definition of scm_round and
4774 rename scm_round_number to scm_round. Likewise for scm_truncate
4775 and scm_truncate_number.
4779 scm_truncate (double x
)
4784 #define trunc scm_truncate
4792 scm_round (double x
)
4794 double plus_half
= x
+ 0.5;
4795 double result
= floor (plus_half
);
4796 /* Adjust so that the scm_round is towards even. */
4797 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4802 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4804 "Round the number @var{x} towards zero.")
4805 #define FUNC_NAME s_scm_truncate_number
4807 if (SCM_FALSEP (scm_negative_p (x
)))
4808 return scm_floor (x
);
4810 return scm_ceiling (x
);
4814 static SCM exactly_one_half
;
4816 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4818 "Round the number @var{x} towards the nearest integer. "
4819 "When it is exactly halfway between two integers, "
4820 "round towards the even one.")
4821 #define FUNC_NAME s_scm_round_number
4823 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4824 SCM result
= scm_floor (plus_half
);
4825 /* Adjust so that the scm_round is towards even. */
4826 if (!SCM_FALSEP (scm_num_eq_p (plus_half
, result
))
4827 && !SCM_FALSEP (scm_odd_p (result
)))
4828 return scm_difference (result
, SCM_MAKINUM (1));
4834 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4836 "Round the number @var{x} towards minus infinity.")
4837 #define FUNC_NAME s_scm_floor
4839 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4841 else if (SCM_REALP (x
))
4842 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
4843 else if (SCM_FRACTIONP (x
))
4845 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4846 SCM_FRACTION_DENOMINATOR (x
));
4847 if (SCM_FALSEP (scm_negative_p (x
)))
4849 /* For positive x, rounding towards zero is correct. */
4854 /* For negative x, we need to return q-1 unless x is an
4855 integer. But fractions are never integer, per our
4857 return scm_difference (q
, SCM_MAKINUM (1));
4861 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
4865 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
4867 "Round the number @var{x} towards infinity.")
4868 #define FUNC_NAME s_scm_ceiling
4870 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4872 else if (SCM_REALP (x
))
4873 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
4874 else if (SCM_FRACTIONP (x
))
4876 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4877 SCM_FRACTION_DENOMINATOR (x
));
4878 if (SCM_FALSEP (scm_positive_p (x
)))
4880 /* For negative x, rounding towards zero is correct. */
4885 /* For positive x, we need to return q+1 unless x is an
4886 integer. But fractions are never integer, per our
4888 return scm_sum (q
, SCM_MAKINUM (1));
4892 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
4896 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
4897 /* "Return the square root of the real number @var{x}."
4899 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
4900 /* "Return the absolute value of the real number @var{x}."
4902 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
4903 /* "Return the @var{x}th power of e."
4905 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
4906 /* "Return the natural logarithm of the real number @var{x}."
4908 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
4909 /* "Return the sine of the real number @var{x}."
4911 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
4912 /* "Return the cosine of the real number @var{x}."
4914 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
4915 /* "Return the tangent of the real number @var{x}."
4917 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
4918 /* "Return the arc sine of the real number @var{x}."
4920 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
4921 /* "Return the arc cosine of the real number @var{x}."
4923 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
4924 /* "Return the arc tangent of the real number @var{x}."
4926 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
4927 /* "Return the hyperbolic sine of the real number @var{x}."
4929 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
4930 /* "Return the hyperbolic cosine of the real number @var{x}."
4932 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
4933 /* "Return the hyperbolic tangent of the real number @var{x}."
4941 static void scm_two_doubles (SCM x
,
4943 const char *sstring
,
4947 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
4950 xy
->x
= SCM_INUM (x
);
4951 else if (SCM_BIGP (x
))
4952 xy
->x
= scm_i_big2dbl (x
);
4953 else if (SCM_REALP (x
))
4954 xy
->x
= SCM_REAL_VALUE (x
);
4955 else if (SCM_FRACTIONP (x
))
4956 xy
->x
= scm_i_fraction2double (x
);
4958 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
4961 xy
->y
= SCM_INUM (y
);
4962 else if (SCM_BIGP (y
))
4963 xy
->y
= scm_i_big2dbl (y
);
4964 else if (SCM_REALP (y
))
4965 xy
->y
= SCM_REAL_VALUE (y
);
4966 else if (SCM_FRACTIONP (y
))
4967 xy
->y
= scm_i_fraction2double (y
);
4969 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
4973 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
4975 "Return @var{x} raised to the power of @var{y}. This\n"
4976 "procedure does not accept complex arguments.")
4977 #define FUNC_NAME s_scm_sys_expt
4980 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4981 return scm_make_real (pow (xy
.x
, xy
.y
));
4986 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
4988 "Return the arc tangent of the two arguments @var{x} and\n"
4989 "@var{y}. This is similar to calculating the arc tangent of\n"
4990 "@var{x} / @var{y}, except that the signs of both arguments\n"
4991 "are used to determine the quadrant of the result. This\n"
4992 "procedure does not accept complex arguments.")
4993 #define FUNC_NAME s_scm_sys_atan2
4996 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4997 return scm_make_real (atan2 (xy
.x
, xy
.y
));
5002 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5003 (SCM real
, SCM imaginary
),
5004 "Return a complex number constructed of the given @var{real} and\n"
5005 "@var{imaginary} parts.")
5006 #define FUNC_NAME s_scm_make_rectangular
5009 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5010 return scm_make_complex (xy
.x
, xy
.y
);
5016 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5018 "Return the complex number @var{x} * e^(i * @var{y}).")
5019 #define FUNC_NAME s_scm_make_polar
5023 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5025 sincos (xy
.y
, &s
, &c
);
5030 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
5035 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5036 /* "Return the real part of the number @var{z}."
5039 scm_real_part (SCM z
)
5043 else if (SCM_BIGP (z
))
5045 else if (SCM_REALP (z
))
5047 else if (SCM_COMPLEXP (z
))
5048 return scm_make_real (SCM_COMPLEX_REAL (z
));
5049 else if (SCM_FRACTIONP (z
))
5052 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5056 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5057 /* "Return the imaginary part of the number @var{z}."
5060 scm_imag_part (SCM z
)
5064 else if (SCM_BIGP (z
))
5066 else if (SCM_REALP (z
))
5068 else if (SCM_COMPLEXP (z
))
5069 return scm_make_real (SCM_COMPLEX_IMAG (z
));
5070 else if (SCM_FRACTIONP (z
))
5073 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5076 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5077 /* "Return the numerator of the number @var{z}."
5080 scm_numerator (SCM z
)
5084 else if (SCM_BIGP (z
))
5086 else if (SCM_FRACTIONP (z
))
5088 scm_i_fraction_reduce (z
);
5089 return SCM_FRACTION_NUMERATOR (z
);
5091 else if (SCM_REALP (z
))
5092 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5094 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5098 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5099 /* "Return the denominator of the number @var{z}."
5102 scm_denominator (SCM z
)
5105 return SCM_MAKINUM (1);
5106 else if (SCM_BIGP (z
))
5107 return SCM_MAKINUM (1);
5108 else if (SCM_FRACTIONP (z
))
5110 scm_i_fraction_reduce (z
);
5111 return SCM_FRACTION_DENOMINATOR (z
);
5113 else if (SCM_REALP (z
))
5114 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5116 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5119 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5120 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5121 * "@code{abs} for real arguments, but also allows complex numbers."
5124 scm_magnitude (SCM z
)
5128 long int zz
= SCM_INUM (z
);
5131 else if (SCM_POSFIXABLE (-zz
))
5132 return SCM_MAKINUM (-zz
);
5134 return scm_i_long2big (-zz
);
5136 else if (SCM_BIGP (z
))
5138 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5139 scm_remember_upto_here_1 (z
);
5141 return scm_i_clonebig (z
, 0);
5145 else if (SCM_REALP (z
))
5146 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5147 else if (SCM_COMPLEXP (z
))
5148 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5149 else if (SCM_FRACTIONP (z
))
5151 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5153 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5154 SCM_FRACTION_DENOMINATOR (z
));
5157 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5161 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5162 /* "Return the angle of the complex number @var{z}."
5167 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5168 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5169 But if atan2 follows the floating point rounding mode, then the value
5170 is not a constant. Maybe it'd be close enough though. */
5173 if (SCM_INUM (z
) >= 0)
5176 return scm_make_real (atan2 (0.0, -1.0));
5178 else if (SCM_BIGP (z
))
5180 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5181 scm_remember_upto_here_1 (z
);
5183 return scm_make_real (atan2 (0.0, -1.0));
5187 else if (SCM_REALP (z
))
5189 if (SCM_REAL_VALUE (z
) >= 0)
5192 return scm_make_real (atan2 (0.0, -1.0));
5194 else if (SCM_COMPLEXP (z
))
5195 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5196 else if (SCM_FRACTIONP (z
))
5198 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5200 else return scm_make_real (atan2 (0.0, -1.0));
5203 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5207 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5208 /* Convert the number @var{x} to its inexact representation.\n"
5211 scm_exact_to_inexact (SCM z
)
5214 return scm_make_real ((double) SCM_INUM (z
));
5215 else if (SCM_BIGP (z
))
5216 return scm_make_real (scm_i_big2dbl (z
));
5217 else if (SCM_FRACTIONP (z
))
5218 return scm_make_real (scm_i_fraction2double (z
));
5219 else if (SCM_INEXACTP (z
))
5222 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5226 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5228 "Return an exact number that is numerically closest to @var{z}.")
5229 #define FUNC_NAME s_scm_inexact_to_exact
5233 else if (SCM_BIGP (z
))
5235 else if (SCM_REALP (z
))
5237 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5238 SCM_OUT_OF_RANGE (1, z
);
5245 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5246 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5247 scm_i_mpz2num (mpq_denref (frac
)));
5249 /* When scm_make_ratio throws, we leak the memory allocated
5256 else if (SCM_FRACTIONP (z
))
5259 SCM_WRONG_TYPE_ARG (1, z
);
5263 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5265 "Return an exact number that is within @var{err} of @var{x}.")
5266 #define FUNC_NAME s_scm_rationalize
5270 else if (SCM_BIGP (x
))
5272 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5274 /* Use continued fractions to find closest ratio. All
5275 arithmetic is done with exact numbers.
5278 SCM ex
= scm_inexact_to_exact (x
);
5279 SCM int_part
= scm_floor (ex
);
5280 SCM tt
= SCM_MAKINUM (1);
5281 SCM a1
= SCM_MAKINUM (0), a2
= SCM_MAKINUM (1), a
= SCM_MAKINUM (0);
5282 SCM b1
= SCM_MAKINUM (1), b2
= SCM_MAKINUM (0), b
= SCM_MAKINUM (0);
5286 if (!SCM_FALSEP (scm_num_eq_p (ex
, int_part
)))
5289 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5290 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5292 /* We stop after a million iterations just to be absolutely sure
5293 that we don't go into an infinite loop. The process normally
5294 converges after less than a dozen iterations.
5297 err
= scm_abs (err
);
5298 while (++i
< 1000000)
5300 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5301 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5302 if (SCM_FALSEP (scm_zero_p (b
)) && /* b != 0 */
5304 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5305 err
))) /* abs(x-a/b) <= err */
5307 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5308 if (SCM_FALSEP (scm_exact_p (x
))
5309 || SCM_FALSEP (scm_exact_p (err
)))
5310 return scm_exact_to_inexact (res
);
5314 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5316 tt
= scm_floor (rx
); /* tt = floor (rx) */
5322 scm_num_overflow (s_scm_rationalize
);
5325 SCM_WRONG_TYPE_ARG (1, x
);
5329 /* if you need to change this, change test-num2integral.c as well */
5330 #if SCM_SIZEOF_LONG_LONG != 0
5332 # define ULLONG_MAX ((unsigned long long) (-1))
5333 # define LLONG_MAX ((long long) (ULLONG_MAX >> 1))
5334 # define LLONG_MIN (~LLONG_MAX)
5338 /* Parameters for creating integer conversion routines.
5340 Define the following preprocessor macros before including
5341 "libguile/num2integral.i.c":
5343 NUM2INTEGRAL - the name of the function for converting from a
5344 Scheme object to the integral type. This function will be
5345 defined when including "num2integral.i.c".
5347 INTEGRAL2NUM - the name of the function for converting from the
5348 integral type to a Scheme object. This function will be defined.
5350 INTEGRAL2BIG - the name of an internal function that createas a
5351 bignum from the integral type. This function will be defined.
5352 The name should start with "scm_i_".
5354 ITYPE - the name of the integral type.
5356 UNSIGNED - Define this to 1 when ITYPE is an unsigned type. Define
5359 UNSIGNED_ITYPE - the name of the the unsigned variant of the
5360 integral type. If you don't define this, it defaults to
5361 "unsigned ITYPE" for signed types and simply "ITYPE" for unsigned
5364 SIZEOF_ITYPE - an expression giving the size of the integral type
5365 in bytes. This expression must be computable by the
5366 preprocessor. (SIZEOF_FOO values are calculated by configure.in
5371 #define NUM2INTEGRAL scm_num2short
5372 #define INTEGRAL2NUM scm_short2num
5373 #define INTEGRAL2BIG scm_i_short2big
5376 #define SIZEOF_ITYPE SIZEOF_SHORT
5377 #include "libguile/num2integral.i.c"
5379 #define NUM2INTEGRAL scm_num2ushort
5380 #define INTEGRAL2NUM scm_ushort2num
5381 #define INTEGRAL2BIG scm_i_ushort2big
5383 #define ITYPE unsigned short
5384 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_SHORT
5385 #include "libguile/num2integral.i.c"
5387 #define NUM2INTEGRAL scm_num2int
5388 #define INTEGRAL2NUM scm_int2num
5389 #define INTEGRAL2BIG scm_i_int2big
5392 #define SIZEOF_ITYPE SIZEOF_INT
5393 #include "libguile/num2integral.i.c"
5395 #define NUM2INTEGRAL scm_num2uint
5396 #define INTEGRAL2NUM scm_uint2num
5397 #define INTEGRAL2BIG scm_i_uint2big
5399 #define ITYPE unsigned int
5400 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_INT
5401 #include "libguile/num2integral.i.c"
5403 #define NUM2INTEGRAL scm_num2long
5404 #define INTEGRAL2NUM scm_long2num
5405 #define INTEGRAL2BIG scm_i_long2big
5408 #define SIZEOF_ITYPE SIZEOF_LONG
5409 #include "libguile/num2integral.i.c"
5411 #define NUM2INTEGRAL scm_num2ulong
5412 #define INTEGRAL2NUM scm_ulong2num
5413 #define INTEGRAL2BIG scm_i_ulong2big
5415 #define ITYPE unsigned long
5416 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG
5417 #include "libguile/num2integral.i.c"
5419 #define NUM2INTEGRAL scm_num2ptrdiff
5420 #define INTEGRAL2NUM scm_ptrdiff2num
5421 #define INTEGRAL2BIG scm_i_ptrdiff2big
5423 #define ITYPE scm_t_ptrdiff
5424 #define UNSIGNED_ITYPE size_t
5425 #define SIZEOF_ITYPE SCM_SIZEOF_SCM_T_PTRDIFF
5426 #include "libguile/num2integral.i.c"
5428 #define NUM2INTEGRAL scm_num2size
5429 #define INTEGRAL2NUM scm_size2num
5430 #define INTEGRAL2BIG scm_i_size2big
5432 #define ITYPE size_t
5433 #define SIZEOF_ITYPE SIZEOF_SIZE_T
5434 #include "libguile/num2integral.i.c"
5436 #if SCM_SIZEOF_LONG_LONG != 0
5438 #ifndef ULONG_LONG_MAX
5439 #define ULONG_LONG_MAX (~0ULL)
5442 #define NUM2INTEGRAL scm_num2long_long
5443 #define INTEGRAL2NUM scm_long_long2num
5444 #define INTEGRAL2BIG scm_i_long_long2big
5446 #define ITYPE long long
5447 #define SIZEOF_ITYPE SIZEOF_LONG_LONG
5448 #include "libguile/num2integral.i.c"
5450 #define NUM2INTEGRAL scm_num2ulong_long
5451 #define INTEGRAL2NUM scm_ulong_long2num
5452 #define INTEGRAL2BIG scm_i_ulong_long2big
5454 #define ITYPE unsigned long long
5455 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG_LONG
5456 #include "libguile/num2integral.i.c"
5458 #endif /* SCM_SIZEOF_LONG_LONG != 0 */
5460 #define NUM2FLOAT scm_num2float
5461 #define FLOAT2NUM scm_float2num
5463 #include "libguile/num2float.i.c"
5465 #define NUM2FLOAT scm_num2double
5466 #define FLOAT2NUM scm_double2num
5467 #define FTYPE double
5468 #include "libguile/num2float.i.c"
5473 #define SIZE_MAX ((size_t) (-1))
5476 #define PTRDIFF_MIN \
5477 ((scm_t_ptrdiff) ((scm_t_ptrdiff) 1 \
5478 << ((sizeof (scm_t_ptrdiff) * SCM_CHAR_BIT) - 1)))
5481 #define PTRDIFF_MAX (~ PTRDIFF_MIN)
5484 #define CHECK(type, v) \
5487 if ((v) != scm_num2##type (scm_##type##2num (v), 1, "check_sanity")) \
5507 CHECK (ptrdiff
, -1);
5509 CHECK (short, SHRT_MAX
);
5510 CHECK (short, SHRT_MIN
);
5511 CHECK (ushort
, USHRT_MAX
);
5512 CHECK (int, INT_MAX
);
5513 CHECK (int, INT_MIN
);
5514 CHECK (uint
, UINT_MAX
);
5515 CHECK (long, LONG_MAX
);
5516 CHECK (long, LONG_MIN
);
5517 CHECK (ulong
, ULONG_MAX
);
5518 CHECK (size
, SIZE_MAX
);
5519 CHECK (ptrdiff
, PTRDIFF_MAX
);
5520 CHECK (ptrdiff
, PTRDIFF_MIN
);
5522 #if SCM_SIZEOF_LONG_LONG != 0
5523 CHECK (long_long
, 0LL);
5524 CHECK (ulong_long
, 0ULL);
5525 CHECK (long_long
, -1LL);
5526 CHECK (long_long
, LLONG_MAX
);
5527 CHECK (long_long
, LLONG_MIN
);
5528 CHECK (ulong_long
, ULLONG_MAX
);
5535 scm_internal_catch (SCM_BOOL_T, check_body, &data, check_handler, &data); \
5536 if (!SCM_FALSEP (data)) abort();
5539 check_body (void *data
)
5541 SCM num
= *(SCM
*) data
;
5542 scm_num2ulong (num
, 1, NULL
);
5544 return SCM_UNSPECIFIED
;
5548 check_handler (void *data
, SCM tag
, SCM throw_args
)
5550 SCM
*num
= (SCM
*) data
;
5553 return SCM_UNSPECIFIED
;
5556 SCM_DEFINE (scm_sys_check_number_conversions
, "%check-number-conversions", 0, 0, 0,
5558 "Number conversion sanity checking.")
5559 #define FUNC_NAME s_scm_sys_check_number_conversions
5561 SCM data
= SCM_MAKINUM (-1);
5563 data
= scm_int2num (INT_MIN
);
5565 data
= scm_ulong2num (ULONG_MAX
);
5566 data
= scm_difference (SCM_INUM0
, data
);
5568 data
= scm_ulong2num (ULONG_MAX
);
5569 data
= scm_sum (SCM_MAKINUM (1), data
); data
= scm_difference (SCM_INUM0
, data
);
5571 data
= scm_int2num (-10000); data
= scm_product (data
, data
); data
= scm_product (data
, data
);
5574 return SCM_UNSPECIFIED
;
5583 abs_most_negative_fixnum
= scm_i_long2big (- SCM_MOST_NEGATIVE_FIXNUM
);
5584 scm_permanent_object (abs_most_negative_fixnum
);
5586 mpz_init_set_si (z_negative_one
, -1);
5588 /* It may be possible to tune the performance of some algorithms by using
5589 * the following constants to avoid the creation of bignums. Please, before
5590 * using these values, remember the two rules of program optimization:
5591 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5592 scm_c_define ("most-positive-fixnum",
5593 SCM_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5594 scm_c_define ("most-negative-fixnum",
5595 SCM_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5597 scm_add_feature ("complex");
5598 scm_add_feature ("inexact");
5599 scm_flo0
= scm_make_real (0.0);
5601 scm_dblprec
= (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5603 { /* determine floating point precision */
5605 double fsum
= 1.0 + f
;
5608 if (++scm_dblprec
> 20)
5616 scm_dblprec
= scm_dblprec
- 1;
5618 #endif /* DBL_DIG */
5624 exactly_one_half
= scm_permanent_object (scm_divide (SCM_MAKINUM (1),
5626 #include "libguile/numbers.x"