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. */
565 guile_Inf
= 1.0 / (tmp
- tmp
);
566 #elif defined (__alpha__) && ! defined (linux)
567 extern unsigned int DINFINITY
[2];
568 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
575 if (guile_Inf
== tmp
)
583 #if defined (HAVE_ISNAN)
585 #if defined (__alpha__) && ! defined (linux)
586 extern unsigned int DQNAN
[2];
587 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
589 guile_NaN
= guile_Inf
/ guile_Inf
;
595 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
598 #define FUNC_NAME s_scm_inf
600 static int initialized
= 0;
606 return scm_make_real (guile_Inf
);
610 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
613 #define FUNC_NAME s_scm_nan
615 static int initialized
= 0;
621 return scm_make_real (guile_NaN
);
626 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
628 "Return the absolute value of @var{x}.")
633 long int xx
= SCM_INUM (x
);
636 else if (SCM_POSFIXABLE (-xx
))
637 return SCM_MAKINUM (-xx
);
639 return scm_i_long2big (-xx
);
641 else if (SCM_BIGP (x
))
643 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
645 return scm_i_clonebig (x
, 0);
649 else if (SCM_REALP (x
))
651 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
652 double xx
= SCM_REAL_VALUE (x
);
654 return scm_make_real (-xx
);
658 else if (SCM_FRACTIONP (x
))
660 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
662 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
663 SCM_FRACTION_DENOMINATOR (x
));
666 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
671 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
672 /* "Return the quotient of the numbers @var{x} and @var{y}."
675 scm_quotient (SCM x
, SCM y
)
679 long xx
= SCM_INUM (x
);
682 long yy
= SCM_INUM (y
);
684 scm_num_overflow (s_quotient
);
689 return SCM_MAKINUM (z
);
691 return scm_i_long2big (z
);
694 else if (SCM_BIGP (y
))
696 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
697 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
698 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
699 return SCM_MAKINUM (-1);
701 return SCM_MAKINUM (0);
704 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
706 else if (SCM_BIGP (x
))
710 long yy
= SCM_INUM (y
);
712 scm_num_overflow (s_quotient
);
717 SCM result
= scm_i_mkbig ();
720 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
723 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
726 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
727 scm_remember_upto_here_1 (x
);
728 return scm_i_normbig (result
);
731 else if (SCM_BIGP (y
))
733 SCM result
= scm_i_mkbig ();
734 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
737 scm_remember_upto_here_2 (x
, y
);
738 return scm_i_normbig (result
);
741 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
744 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
747 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
748 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
750 * "(remainder 13 4) @result{} 1\n"
751 * "(remainder -13 4) @result{} -1\n"
755 scm_remainder (SCM x
, SCM y
)
761 long yy
= SCM_INUM (y
);
763 scm_num_overflow (s_remainder
);
766 long z
= SCM_INUM (x
) % yy
;
767 return SCM_MAKINUM (z
);
770 else if (SCM_BIGP (y
))
772 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
773 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
774 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
775 return SCM_MAKINUM (0);
780 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
782 else if (SCM_BIGP (x
))
786 long yy
= SCM_INUM (y
);
788 scm_num_overflow (s_remainder
);
791 SCM result
= scm_i_mkbig ();
794 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
795 scm_remember_upto_here_1 (x
);
796 return scm_i_normbig (result
);
799 else if (SCM_BIGP (y
))
801 SCM result
= scm_i_mkbig ();
802 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
805 scm_remember_upto_here_2 (x
, y
);
806 return scm_i_normbig (result
);
809 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
812 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
816 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
817 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
819 * "(modulo 13 4) @result{} 1\n"
820 * "(modulo -13 4) @result{} 3\n"
824 scm_modulo (SCM x
, SCM y
)
828 long xx
= SCM_INUM (x
);
831 long yy
= SCM_INUM (y
);
833 scm_num_overflow (s_modulo
);
836 /* FIXME: I think this may be a bug on some arches -- results
837 of % with negative second arg are undefined... */
855 return SCM_MAKINUM (result
);
858 else if (SCM_BIGP (y
))
860 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
863 scm_num_overflow (s_modulo
);
871 SCM pos_y
= scm_i_clonebig (y
, 0);
872 /* do this after the last scm_op */
873 mpz_init_set_si (z_x
, xx
);
874 result
= pos_y
; /* re-use this bignum */
875 mpz_mod (SCM_I_BIG_MPZ (result
),
877 SCM_I_BIG_MPZ (pos_y
));
878 scm_remember_upto_here_1 (pos_y
);
882 result
= scm_i_mkbig ();
883 /* do this after the last scm_op */
884 mpz_init_set_si (z_x
, xx
);
885 mpz_mod (SCM_I_BIG_MPZ (result
),
888 scm_remember_upto_here_1 (y
);
891 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
892 mpz_add (SCM_I_BIG_MPZ (result
),
894 SCM_I_BIG_MPZ (result
));
895 scm_remember_upto_here_1 (y
);
896 /* and do this before the next one */
898 return scm_i_normbig (result
);
902 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
904 else if (SCM_BIGP (x
))
908 long yy
= SCM_INUM (y
);
910 scm_num_overflow (s_modulo
);
913 SCM result
= scm_i_mkbig ();
914 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
916 (yy
< 0) ? - yy
: yy
);
917 scm_remember_upto_here_1 (x
);
918 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
919 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
920 SCM_I_BIG_MPZ (result
),
922 return scm_i_normbig (result
);
925 else if (SCM_BIGP (y
))
927 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
929 scm_num_overflow (s_modulo
);
932 SCM result
= scm_i_mkbig ();
933 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
934 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
935 mpz_mod (SCM_I_BIG_MPZ (result
),
937 SCM_I_BIG_MPZ (pos_y
));
939 scm_remember_upto_here_1 (x
);
940 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
941 mpz_add (SCM_I_BIG_MPZ (result
),
943 SCM_I_BIG_MPZ (result
));
944 scm_remember_upto_here_2 (y
, pos_y
);
945 return scm_i_normbig (result
);
949 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
952 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
955 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
956 /* "Return the greatest common divisor of all arguments.\n"
957 * "If called without arguments, 0 is returned."
960 scm_gcd (SCM x
, SCM y
)
963 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
969 long xx
= SCM_INUM (x
);
970 long yy
= SCM_INUM (y
);
971 long u
= xx
< 0 ? -xx
: xx
;
972 long v
= yy
< 0 ? -yy
: yy
;
982 /* Determine a common factor 2^k */
983 while (!(1 & (u
| v
)))
989 /* Now, any factor 2^n can be eliminated */
1009 return (SCM_POSFIXABLE (result
)
1010 ? SCM_MAKINUM (result
)
1011 : scm_i_long2big (result
));
1013 else if (SCM_BIGP (y
))
1015 SCM result
= scm_i_mkbig ();
1016 SCM mx
= scm_i_mkbig ();
1017 mpz_set_si (SCM_I_BIG_MPZ (mx
), SCM_INUM (x
));
1018 scm_remember_upto_here_1 (x
);
1019 mpz_gcd (SCM_I_BIG_MPZ (result
),
1022 scm_remember_upto_here_2 (mx
, y
);
1023 return scm_i_normbig (result
);
1026 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1028 else if (SCM_BIGP (x
))
1032 unsigned long result
;
1033 long yy
= SCM_INUM (y
);
1038 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1039 scm_remember_upto_here_1 (x
);
1040 return (SCM_POSFIXABLE (result
)
1041 ? SCM_MAKINUM (result
)
1042 : scm_ulong2num (result
));
1044 else if (SCM_BIGP (y
))
1046 SCM result
= scm_i_mkbig ();
1047 mpz_gcd (SCM_I_BIG_MPZ (result
),
1050 scm_remember_upto_here_2 (x
, y
);
1051 return scm_i_normbig (result
);
1054 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1057 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1060 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1061 /* "Return the least common multiple of the arguments.\n"
1062 * "If called without arguments, 1 is returned."
1065 scm_lcm (SCM n1
, SCM n2
)
1067 if (SCM_UNBNDP (n2
))
1069 if (SCM_UNBNDP (n1
))
1070 return SCM_MAKINUM (1L);
1071 n2
= SCM_MAKINUM (1L);
1074 SCM_GASSERT2 (SCM_INUMP (n1
) || SCM_BIGP (n1
),
1075 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1076 SCM_GASSERT2 (SCM_INUMP (n2
) || SCM_BIGP (n2
),
1077 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1083 SCM d
= scm_gcd (n1
, n2
);
1084 if (SCM_EQ_P (d
, SCM_INUM0
))
1087 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1091 /* inum n1, big n2 */
1094 SCM result
= scm_i_mkbig ();
1095 long nn1
= SCM_INUM (n1
);
1096 if (nn1
== 0) return SCM_INUM0
;
1097 if (nn1
< 0) nn1
= - nn1
;
1098 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1099 scm_remember_upto_here_1 (n2
);
1114 SCM result
= scm_i_mkbig ();
1115 mpz_lcm(SCM_I_BIG_MPZ (result
),
1117 SCM_I_BIG_MPZ (n2
));
1118 scm_remember_upto_here_2(n1
, n2
);
1119 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1125 #ifndef scm_long2num
1126 #define SCM_LOGOP_RETURN(x) scm_ulong2num(x)
1128 #define SCM_LOGOP_RETURN(x) SCM_MAKINUM(x)
1131 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1136 + + + x (map digit:logand X Y)
1137 + - + x (map digit:logand X (lognot (+ -1 Y)))
1138 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1139 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1144 + + + (map digit:logior X Y)
1145 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1146 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1147 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1152 + + + (map digit:logxor X Y)
1153 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1154 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1155 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1160 + + (any digit:logand X Y)
1161 + - (any digit:logand X (lognot (+ -1 Y)))
1162 - + (any digit:logand (lognot (+ -1 X)) Y)
1167 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1169 "Return the bitwise AND of the integer arguments.\n\n"
1171 "(logand) @result{} -1\n"
1172 "(logand 7) @result{} 7\n"
1173 "(logand #b111 #b011 #b001) @result{} 1\n"
1175 #define FUNC_NAME s_scm_logand
1179 if (SCM_UNBNDP (n2
))
1181 if (SCM_UNBNDP (n1
))
1182 return SCM_MAKINUM (-1);
1183 else if (!SCM_NUMBERP (n1
))
1184 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1185 else if (SCM_NUMBERP (n1
))
1188 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1193 nn1
= SCM_INUM (n1
);
1196 long nn2
= SCM_INUM (n2
);
1197 return SCM_MAKINUM (nn1
& nn2
);
1199 else if SCM_BIGP (n2
)
1205 SCM result_z
= scm_i_mkbig ();
1207 mpz_init_set_si (nn1_z
, nn1
);
1208 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1209 scm_remember_upto_here_1 (n2
);
1211 return scm_i_normbig (result_z
);
1215 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1217 else if (SCM_BIGP (n1
))
1222 nn1
= SCM_INUM (n1
);
1225 else if (SCM_BIGP (n2
))
1227 SCM result_z
= scm_i_mkbig ();
1228 mpz_and (SCM_I_BIG_MPZ (result_z
),
1230 SCM_I_BIG_MPZ (n2
));
1231 scm_remember_upto_here_2 (n1
, n2
);
1232 return scm_i_normbig (result_z
);
1235 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1238 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1243 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1245 "Return the bitwise OR of the integer arguments.\n\n"
1247 "(logior) @result{} 0\n"
1248 "(logior 7) @result{} 7\n"
1249 "(logior #b000 #b001 #b011) @result{} 3\n"
1251 #define FUNC_NAME s_scm_logior
1255 if (SCM_UNBNDP (n2
))
1257 if (SCM_UNBNDP (n1
))
1259 else if (SCM_NUMBERP (n1
))
1262 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1267 nn1
= SCM_INUM (n1
);
1270 long nn2
= SCM_INUM (n2
);
1271 return SCM_MAKINUM (nn1
| nn2
);
1273 else if (SCM_BIGP (n2
))
1279 SCM result_z
= scm_i_mkbig ();
1281 mpz_init_set_si (nn1_z
, nn1
);
1282 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1283 scm_remember_upto_here_1 (n2
);
1289 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1291 else if (SCM_BIGP (n1
))
1296 nn1
= SCM_INUM (n1
);
1299 else if (SCM_BIGP (n2
))
1301 SCM result_z
= scm_i_mkbig ();
1302 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1304 SCM_I_BIG_MPZ (n2
));
1305 scm_remember_upto_here_2 (n1
, n2
);
1309 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1312 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1317 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1319 "Return the bitwise XOR of the integer arguments. A bit is\n"
1320 "set in the result if it is set in an odd number of arguments.\n"
1322 "(logxor) @result{} 0\n"
1323 "(logxor 7) @result{} 7\n"
1324 "(logxor #b000 #b001 #b011) @result{} 2\n"
1325 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1327 #define FUNC_NAME s_scm_logxor
1331 if (SCM_UNBNDP (n2
))
1333 if (SCM_UNBNDP (n1
))
1335 else if (SCM_NUMBERP (n1
))
1338 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1343 nn1
= SCM_INUM (n1
);
1346 long nn2
= SCM_INUM (n2
);
1347 return SCM_MAKINUM (nn1
^ nn2
);
1349 else if (SCM_BIGP (n2
))
1353 SCM result_z
= scm_i_mkbig ();
1355 mpz_init_set_si (nn1_z
, nn1
);
1356 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1357 scm_remember_upto_here_1 (n2
);
1359 return scm_i_normbig (result_z
);
1363 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1365 else if (SCM_BIGP (n1
))
1370 nn1
= SCM_INUM (n1
);
1373 else if (SCM_BIGP (n2
))
1375 SCM result_z
= scm_i_mkbig ();
1376 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1378 SCM_I_BIG_MPZ (n2
));
1379 scm_remember_upto_here_2 (n1
, n2
);
1380 return scm_i_normbig (result_z
);
1383 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1386 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1391 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1394 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1395 "(logtest #b0100 #b1011) @result{} #f\n"
1396 "(logtest #b0100 #b0111) @result{} #t\n"
1398 #define FUNC_NAME s_scm_logtest
1407 long nk
= SCM_INUM (k
);
1408 return SCM_BOOL (nj
& nk
);
1410 else if (SCM_BIGP (k
))
1418 mpz_init_set_si (nj_z
, nj
);
1419 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1420 scm_remember_upto_here_1 (k
);
1421 result
= SCM_BOOL (mpz_sgn (nj_z
) != 0);
1427 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1429 else if (SCM_BIGP (j
))
1437 else if (SCM_BIGP (k
))
1441 mpz_init (result_z
);
1445 scm_remember_upto_here_2 (j
, k
);
1446 result
= SCM_BOOL (mpz_sgn (result_z
) != 0);
1447 mpz_clear (result_z
);
1451 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1454 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1459 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1462 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1463 "(logbit? 0 #b1101) @result{} #t\n"
1464 "(logbit? 1 #b1101) @result{} #f\n"
1465 "(logbit? 2 #b1101) @result{} #t\n"
1466 "(logbit? 3 #b1101) @result{} #t\n"
1467 "(logbit? 4 #b1101) @result{} #f\n"
1469 #define FUNC_NAME s_scm_logbit_p
1471 unsigned long int iindex
;
1473 SCM_VALIDATE_INUM_MIN (SCM_ARG1
, index
, 0);
1474 iindex
= (unsigned long int) SCM_INUM (index
);
1477 return SCM_BOOL ((1L << iindex
) & SCM_INUM (j
));
1478 else if (SCM_BIGP (j
))
1480 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1481 scm_remember_upto_here_1 (j
);
1482 return SCM_BOOL (val
);
1485 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1490 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1492 "Return the integer which is the ones-complement of the integer\n"
1496 "(number->string (lognot #b10000000) 2)\n"
1497 " @result{} \"-10000001\"\n"
1498 "(number->string (lognot #b0) 2)\n"
1499 " @result{} \"-1\"\n"
1501 #define FUNC_NAME s_scm_lognot
1503 if (SCM_INUMP (n
)) {
1504 /* No overflow here, just need to toggle all the bits making up the inum.
1505 Enhancement: No need to strip the tag and add it back, could just xor
1506 a block of 1 bits, if that worked with the various debug versions of
1508 return SCM_MAKINUM (~ SCM_INUM (n
));
1510 } else if (SCM_BIGP (n
)) {
1511 SCM result
= scm_i_mkbig ();
1512 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1513 scm_remember_upto_here_1 (n
);
1517 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1522 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1524 "Return @var{n} raised to the non-negative integer exponent\n"
1528 "(integer-expt 2 5)\n"
1530 "(integer-expt -3 3)\n"
1533 #define FUNC_NAME s_scm_integer_expt
1536 SCM z_i2
= SCM_BOOL_F
;
1538 SCM acc
= SCM_MAKINUM (1L);
1540 /* 0^0 == 1 according to R5RS */
1541 if (SCM_EQ_P (n
, SCM_INUM0
) || SCM_EQ_P (n
, acc
))
1542 return SCM_FALSEP (scm_zero_p(k
)) ? n
: acc
;
1543 else if (SCM_EQ_P (n
, SCM_MAKINUM (-1L)))
1544 return SCM_FALSEP (scm_even_p (k
)) ? n
: acc
;
1548 else if (SCM_BIGP (k
))
1550 z_i2
= scm_i_clonebig (k
, 1);
1551 scm_remember_upto_here_1 (k
);
1554 else if (SCM_REALP (k
))
1556 double r
= SCM_REAL_VALUE (k
);
1558 SCM_WRONG_TYPE_ARG (2, k
);
1559 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1561 z_i2
= scm_i_mkbig ();
1562 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1571 SCM_WRONG_TYPE_ARG (2, k
);
1575 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1577 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1578 n
= scm_divide (n
, SCM_UNDEFINED
);
1582 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1586 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1588 return scm_product (acc
, n
);
1590 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1591 acc
= scm_product (acc
, n
);
1592 n
= scm_product (n
, n
);
1593 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1601 n
= scm_divide (n
, SCM_UNDEFINED
);
1608 return scm_product (acc
, n
);
1610 acc
= scm_product (acc
, n
);
1611 n
= scm_product (n
, n
);
1618 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1620 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1621 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1623 "This is effectively a multiplication by 2^@var{cnt}}, and when\n"
1624 "@var{cnt} is negative it's a division, rounded towards negative\n"
1625 "infinity. (Note that this is not the same rounding as\n"
1626 "@code{quotient} does.)\n"
1628 "With @var{n} viewed as an infinite precision twos complement,\n"
1629 "@code{ash} means a left shift introducing zero bits, or a right\n"
1630 "shift dropping bits.\n"
1633 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1634 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1636 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1637 "(ash -23 -2) @result{} -6\n"
1639 #define FUNC_NAME s_scm_ash
1643 SCM_VALIDATE_INUM (2, cnt
);
1645 bits_to_shift
= SCM_INUM (cnt
);
1647 if (bits_to_shift
< 0)
1649 /* Shift right by abs(cnt) bits. This is realized as a division
1650 by div:=2^abs(cnt). However, to guarantee the floor
1651 rounding, negative values require some special treatment.
1653 SCM div
= scm_integer_expt (SCM_MAKINUM (2),
1654 SCM_MAKINUM (-bits_to_shift
));
1656 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1657 if (SCM_FALSEP (scm_negative_p (n
)))
1658 return scm_quotient (n
, div
);
1660 return scm_sum (SCM_MAKINUM (-1L),
1661 scm_quotient (scm_sum (SCM_MAKINUM (1L), n
), div
));
1664 /* Shift left is done by multiplication with 2^CNT */
1665 return scm_product (n
, scm_integer_expt (SCM_MAKINUM (2), cnt
));
1670 #define MIN(x,y) ((x) < (y) ? (x) : (y))
1672 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1673 (SCM n
, SCM start
, SCM end
),
1674 "Return the integer composed of the @var{start} (inclusive)\n"
1675 "through @var{end} (exclusive) bits of @var{n}. The\n"
1676 "@var{start}th bit becomes the 0-th bit in the result.\n"
1679 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1680 " @result{} \"1010\"\n"
1681 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1682 " @result{} \"10110\"\n"
1684 #define FUNC_NAME s_scm_bit_extract
1686 unsigned long int istart
, iend
, bits
;
1687 SCM_VALIDATE_INUM_MIN_COPY (2, start
,0, istart
);
1688 SCM_VALIDATE_INUM_MIN_COPY (3, end
, 0, iend
);
1689 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1691 /* how many bits to keep */
1692 bits
= iend
- istart
;
1696 long int in
= SCM_INUM (n
);
1698 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1699 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in".
1700 FIXME: This shift relies on signed right shifts being arithmetic,
1701 which is not guaranteed by C99. */
1702 in
>>= MIN (istart
, SCM_I_FIXNUM_BIT
-1);
1704 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1706 /* Since we emulate two's complement encoded numbers, this
1707 * special case requires us to produce a result that has
1708 * more bits than can be stored in a fixnum.
1710 SCM result
= scm_i_long2big (in
);
1711 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1716 /* mask down to requisite bits */
1717 bits
= MIN (bits
, SCM_I_FIXNUM_BIT
);
1718 return SCM_MAKINUM (in
& ((1L << bits
) - 1));
1720 else if (SCM_BIGP (n
))
1725 result
= SCM_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1729 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1730 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1731 such bits into a ulong. */
1732 result
= scm_i_mkbig ();
1733 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1734 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1735 result
= scm_i_normbig (result
);
1737 scm_remember_upto_here_1 (n
);
1741 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1746 static const char scm_logtab
[] = {
1747 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1750 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1752 "Return the number of bits in integer @var{n}. If integer is\n"
1753 "positive, the 1-bits in its binary representation are counted.\n"
1754 "If negative, the 0-bits in its two's-complement binary\n"
1755 "representation are counted. If 0, 0 is returned.\n"
1758 "(logcount #b10101010)\n"
1765 #define FUNC_NAME s_scm_logcount
1769 unsigned long int c
= 0;
1770 long int nn
= SCM_INUM (n
);
1775 c
+= scm_logtab
[15 & nn
];
1778 return SCM_MAKINUM (c
);
1780 else if (SCM_BIGP (n
))
1782 unsigned long count
;
1783 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1784 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1786 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1787 scm_remember_upto_here_1 (n
);
1788 return SCM_MAKINUM (count
);
1791 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1796 static const char scm_ilentab
[] = {
1797 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1801 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1803 "Return the number of bits necessary to represent @var{n}.\n"
1806 "(integer-length #b10101010)\n"
1808 "(integer-length 0)\n"
1810 "(integer-length #b1111)\n"
1813 #define FUNC_NAME s_scm_integer_length
1817 unsigned long int c
= 0;
1819 long int nn
= SCM_INUM (n
);
1825 l
= scm_ilentab
[15 & nn
];
1828 return SCM_MAKINUM (c
- 4 + l
);
1830 else if (SCM_BIGP (n
))
1832 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1833 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1834 1 too big, so check for that and adjust. */
1835 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1836 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1837 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1838 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1840 scm_remember_upto_here_1 (n
);
1841 return SCM_MAKINUM (size
);
1844 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1848 /*** NUMBERS -> STRINGS ***/
1850 static const double fx
[] =
1851 { 0.0, 5e-1, 5e-2, 5e-3, 5e-4, 5e-5,
1852 5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
1853 5e-11, 5e-12, 5e-13, 5e-14, 5e-15,
1854 5e-16, 5e-17, 5e-18, 5e-19, 5e-20};
1857 idbl2str (double f
, char *a
)
1859 int efmt
, dpt
, d
, i
, wp
= scm_dblprec
;
1865 #ifdef HAVE_COPYSIGN
1866 double sgn
= copysign (1.0, f
);
1872 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
1878 strcpy (a
, "-inf.0");
1880 strcpy (a
, "+inf.0");
1883 else if (xisnan (f
))
1885 strcpy (a
, "+nan.0");
1895 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
1896 make-uniform-vector, from causing infinite loops. */
1900 if (exp
-- < DBL_MIN_10_EXP
)
1911 if (exp
++ > DBL_MAX_10_EXP
)
1931 if (f
+ fx
[wp
] >= 10.0)
1938 dpt
= (exp
+ 9999) % 3;
1942 efmt
= (exp
< -3) || (exp
> wp
+ 2);
1967 if (f
+ fx
[wp
] >= 1.0)
1981 if ((dpt
> 4) && (exp
> 6))
1983 d
= (a
[0] == '-' ? 2 : 1);
1984 for (i
= ch
++; i
> d
; i
--)
1997 if (a
[ch
- 1] == '.')
1998 a
[ch
++] = '0'; /* trailing zero */
2007 for (i
= 10; i
<= exp
; i
*= 10);
2008 for (i
/= 10; i
; i
/= 10)
2010 a
[ch
++] = exp
/ i
+ '0';
2019 iflo2str (SCM flt
, char *str
)
2022 if (SCM_REALP (flt
))
2023 i
= idbl2str (SCM_REAL_VALUE (flt
), str
);
2026 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
);
2027 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2029 double imag
= SCM_COMPLEX_IMAG (flt
);
2030 /* Don't output a '+' for negative numbers or for Inf and
2031 NaN. They will provide their own sign. */
2032 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2034 i
+= idbl2str (imag
, &str
[i
]);
2041 /* convert a long to a string (unterminated). returns the number of
2042 characters in the result.
2044 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2046 scm_iint2str (long num
, int rad
, char *p
)
2050 unsigned long n
= (num
< 0) ? -num
: num
;
2052 for (n
/= rad
; n
> 0; n
/= rad
)
2069 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2074 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2076 "Return a string holding the external representation of the\n"
2077 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2078 "inexact, a radix of 10 will be used.")
2079 #define FUNC_NAME s_scm_number_to_string
2083 if (SCM_UNBNDP (radix
))
2087 SCM_VALIDATE_INUM (2, radix
);
2088 base
= SCM_INUM (radix
);
2089 /* FIXME: ask if range limit was OK, and if so, document */
2090 SCM_ASSERT_RANGE (2, radix
, (base
>= 2) && (base
<= 36));
2095 char num_buf
[SCM_INTBUFLEN
];
2096 size_t length
= scm_iint2str (SCM_INUM (n
), base
, num_buf
);
2097 return scm_mem2string (num_buf
, length
);
2099 else if (SCM_BIGP (n
))
2101 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2102 scm_remember_upto_here_1 (n
);
2103 return scm_take0str (str
);
2105 else if (SCM_FRACTIONP (n
))
2107 scm_i_fraction_reduce (n
);
2108 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2109 scm_mem2string ("/", 1),
2110 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2112 else if (SCM_INEXACTP (n
))
2114 char num_buf
[FLOBUFLEN
];
2115 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
));
2118 SCM_WRONG_TYPE_ARG (1, n
);
2123 /* These print routines used to be stubbed here so that scm_repl.c
2124 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2127 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2129 char num_buf
[FLOBUFLEN
];
2130 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2135 scm_print_complex (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_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2147 scm_i_fraction_reduce (sexp
);
2148 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2149 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2150 scm_remember_upto_here_1 (str
);
2155 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2157 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2158 scm_remember_upto_here_1 (exp
);
2159 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2163 /*** END nums->strs ***/
2166 /*** STRINGS -> NUMBERS ***/
2168 /* The following functions implement the conversion from strings to numbers.
2169 * The implementation somehow follows the grammar for numbers as it is given
2170 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2171 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2172 * points should be noted about the implementation:
2173 * * Each function keeps a local index variable 'idx' that points at the
2174 * current position within the parsed string. The global index is only
2175 * updated if the function could parse the corresponding syntactic unit
2177 * * Similarly, the functions keep track of indicators of inexactness ('#',
2178 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2179 * global exactness information is only updated after each part has been
2180 * successfully parsed.
2181 * * Sequences of digits are parsed into temporary variables holding fixnums.
2182 * Only if these fixnums would overflow, the result variables are updated
2183 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2184 * the temporary variables holding the fixnums are cleared, and the process
2185 * starts over again. If for example fixnums were able to store five decimal
2186 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2187 * and the result was computed as 12345 * 100000 + 67890. In other words,
2188 * only every five digits two bignum operations were performed.
2191 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2193 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2195 /* In non ASCII-style encodings the following macro might not work. */
2196 #define XDIGIT2UINT(d) (isdigit (d) ? (d) - '0' : tolower (d) - 'a' + 10)
2199 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2200 unsigned int radix
, enum t_exactness
*p_exactness
)
2202 unsigned int idx
= *p_idx
;
2203 unsigned int hash_seen
= 0;
2204 scm_t_bits shift
= 1;
2206 unsigned int digit_value
;
2216 digit_value
= XDIGIT2UINT (c
);
2217 if (digit_value
>= radix
)
2221 result
= SCM_MAKINUM (digit_value
);
2229 digit_value
= XDIGIT2UINT (c
);
2230 if (digit_value
>= radix
)
2242 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2244 result
= scm_product (result
, SCM_MAKINUM (shift
));
2246 result
= scm_sum (result
, SCM_MAKINUM (add
));
2253 shift
= shift
* radix
;
2254 add
= add
* radix
+ digit_value
;
2259 result
= scm_product (result
, SCM_MAKINUM (shift
));
2261 result
= scm_sum (result
, SCM_MAKINUM (add
));
2265 *p_exactness
= INEXACT
;
2271 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2272 * covers the parts of the rules that start at a potential point. The value
2273 * of the digits up to the point have been parsed by the caller and are given
2274 * in variable result. The content of *p_exactness indicates, whether a hash
2275 * has already been seen in the digits before the point.
2278 /* In non ASCII-style encodings the following macro might not work. */
2279 #define DIGIT2UINT(d) ((d) - '0')
2282 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2283 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2285 unsigned int idx
= *p_idx
;
2286 enum t_exactness x
= *p_exactness
;
2291 if (mem
[idx
] == '.')
2293 scm_t_bits shift
= 1;
2295 unsigned int digit_value
;
2296 SCM big_shift
= SCM_MAKINUM (1);
2307 digit_value
= DIGIT2UINT (c
);
2318 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2320 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2321 result
= scm_product (result
, SCM_MAKINUM (shift
));
2323 result
= scm_sum (result
, SCM_MAKINUM (add
));
2331 add
= add
* 10 + digit_value
;
2337 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2338 result
= scm_product (result
, SCM_MAKINUM (shift
));
2339 result
= scm_sum (result
, SCM_MAKINUM (add
));
2342 result
= scm_divide (result
, big_shift
);
2344 /* We've seen a decimal point, thus the value is implicitly inexact. */
2356 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2387 exponent
= DIGIT2UINT (c
);
2394 if (exponent
<= SCM_MAXEXP
)
2395 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2401 if (exponent
> SCM_MAXEXP
)
2403 size_t exp_len
= idx
- start
;
2404 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2405 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2406 scm_out_of_range ("string->number", exp_num
);
2409 e
= scm_integer_expt (SCM_MAKINUM (10), SCM_MAKINUM (exponent
));
2411 result
= scm_product (result
, e
);
2413 result
= scm_divide2real (result
, e
);
2415 /* We've seen an exponent, thus the value is implicitly inexact. */
2433 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2436 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2437 unsigned int radix
, enum t_exactness
*p_exactness
)
2439 unsigned int idx
= *p_idx
;
2445 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2451 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2453 enum t_exactness x
= EXACT
;
2455 /* Cobble up the fractional part. We might want to set the
2456 NaN's mantissa from it. */
2458 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2463 if (mem
[idx
] == '.')
2467 else if (idx
+ 1 == len
)
2469 else if (!isdigit (mem
[idx
+ 1]))
2472 result
= mem2decimal_from_point (SCM_MAKINUM (0), mem
, len
,
2473 p_idx
, p_exactness
);
2477 enum t_exactness x
= EXACT
;
2480 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2481 if (SCM_FALSEP (uinteger
))
2486 else if (mem
[idx
] == '/')
2492 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2493 if (SCM_FALSEP (divisor
))
2496 /* both are int/big here, I assume */
2497 result
= scm_make_ratio (uinteger
, divisor
);
2499 else if (radix
== 10)
2501 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2502 if (SCM_FALSEP (result
))
2513 /* When returning an inexact zero, make sure it is represented as a
2514 floating point value so that we can change its sign.
2516 if (SCM_EQ_P (result
, SCM_MAKINUM(0)) && *p_exactness
== INEXACT
)
2517 result
= scm_make_real (0.0);
2523 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2526 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2527 unsigned int radix
, enum t_exactness
*p_exactness
)
2551 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2552 if (SCM_FALSEP (ureal
))
2554 /* input must be either +i or -i */
2559 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2565 return scm_make_rectangular (SCM_MAKINUM (0), SCM_MAKINUM (sign
));
2572 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2573 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2582 /* either +<ureal>i or -<ureal>i */
2589 return scm_make_rectangular (SCM_MAKINUM (0), ureal
);
2592 /* polar input: <real>@<real>. */
2617 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2618 if (SCM_FALSEP (angle
))
2623 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2624 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2626 result
= scm_make_polar (ureal
, angle
);
2631 /* expecting input matching <real>[+-]<ureal>?i */
2638 int sign
= (c
== '+') ? 1 : -1;
2639 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2641 if (SCM_FALSEP (imag
))
2642 imag
= SCM_MAKINUM (sign
);
2643 else if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2644 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2648 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2655 return scm_make_rectangular (ureal
, imag
);
2664 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2666 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2669 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2671 unsigned int idx
= 0;
2672 unsigned int radix
= NO_RADIX
;
2673 enum t_exactness forced_x
= NO_EXACTNESS
;
2674 enum t_exactness implicit_x
= EXACT
;
2677 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2678 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2680 switch (mem
[idx
+ 1])
2683 if (radix
!= NO_RADIX
)
2688 if (radix
!= NO_RADIX
)
2693 if (forced_x
!= NO_EXACTNESS
)
2698 if (forced_x
!= NO_EXACTNESS
)
2703 if (radix
!= NO_RADIX
)
2708 if (radix
!= NO_RADIX
)
2718 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2719 if (radix
== NO_RADIX
)
2720 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2722 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2724 if (SCM_FALSEP (result
))
2730 if (SCM_INEXACTP (result
))
2731 return scm_inexact_to_exact (result
);
2735 if (SCM_INEXACTP (result
))
2738 return scm_exact_to_inexact (result
);
2741 if (implicit_x
== INEXACT
)
2743 if (SCM_INEXACTP (result
))
2746 return scm_exact_to_inexact (result
);
2754 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2755 (SCM string
, SCM radix
),
2756 "Return a number of the maximally precise representation\n"
2757 "expressed by the given @var{string}. @var{radix} must be an\n"
2758 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2759 "is a default radix that may be overridden by an explicit radix\n"
2760 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2761 "supplied, then the default radix is 10. If string is not a\n"
2762 "syntactically valid notation for a number, then\n"
2763 "@code{string->number} returns @code{#f}.")
2764 #define FUNC_NAME s_scm_string_to_number
2768 SCM_VALIDATE_STRING (1, string
);
2769 SCM_VALIDATE_INUM_MIN_DEF_COPY (2, radix
,2,10, base
);
2770 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2771 SCM_STRING_LENGTH (string
),
2773 return scm_return_first (answer
, string
);
2778 /*** END strs->nums ***/
2782 scm_make_real (double x
)
2784 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2786 SCM_REAL_VALUE (z
) = x
;
2792 scm_make_complex (double x
, double y
)
2795 return scm_make_real (x
);
2799 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2801 SCM_COMPLEX_REAL (z
) = x
;
2802 SCM_COMPLEX_IMAG (z
) = y
;
2809 scm_bigequal (SCM x
, SCM y
)
2811 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2812 scm_remember_upto_here_2 (x
, y
);
2813 return SCM_BOOL (0 == result
);
2817 scm_real_equalp (SCM x
, SCM y
)
2819 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2823 scm_complex_equalp (SCM x
, SCM y
)
2825 return SCM_BOOL (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2826 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
2830 scm_i_fraction_equalp (SCM x
, SCM y
)
2832 scm_i_fraction_reduce (x
);
2833 scm_i_fraction_reduce (y
);
2834 if (SCM_FALSEP (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
2835 SCM_FRACTION_NUMERATOR (y
)))
2836 || SCM_FALSEP (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
2837 SCM_FRACTION_DENOMINATOR (y
))))
2844 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
2845 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
2846 * "else. Note that the sets of complex, real, rational and\n"
2847 * "integer values form subsets of the set of numbers, i. e. the\n"
2848 * "predicate will be fulfilled for any number."
2850 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
2852 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
2853 "otherwise. Note that the sets of real, rational and integer\n"
2854 "values form subsets of the set of complex numbers, i. e. the\n"
2855 "predicate will also be fulfilled if @var{x} is a real,\n"
2856 "rational or integer number.")
2857 #define FUNC_NAME s_scm_number_p
2859 return SCM_BOOL (SCM_NUMBERP (x
));
2864 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
2866 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
2867 "otherwise. Note that the set of integer values forms a subset of\n"
2868 "the set of real numbers, i. e. the predicate will also be\n"
2869 "fulfilled if @var{x} is an integer number.")
2870 #define FUNC_NAME s_scm_real_p
2872 /* we can't represent irrational numbers. */
2873 return scm_rational_p (x
);
2877 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
2879 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
2880 "otherwise. Note that the set of integer values forms a subset of\n"
2881 "the set of rational numbers, i. e. the predicate will also be\n"
2882 "fulfilled if @var{x} is an integer number.")
2883 #define FUNC_NAME s_scm_rational_p
2887 else if (SCM_IMP (x
))
2889 else if (SCM_BIGP (x
))
2891 else if (SCM_FRACTIONP (x
))
2893 else if (SCM_REALP (x
))
2894 /* due to their limited precision, all floating point numbers are
2895 rational as well. */
2903 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
2905 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
2907 #define FUNC_NAME s_scm_integer_p
2916 if (!SCM_INEXACTP (x
))
2918 if (SCM_COMPLEXP (x
))
2920 r
= SCM_REAL_VALUE (x
);
2928 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
2930 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
2932 #define FUNC_NAME s_scm_inexact_p
2934 if (SCM_INEXACTP (x
))
2936 if (SCM_NUMBERP (x
))
2938 SCM_WRONG_TYPE_ARG (1, x
);
2943 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
2944 /* "Return @code{#t} if all parameters are numerically equal." */
2946 scm_num_eq_p (SCM x
, SCM y
)
2951 long xx
= SCM_INUM (x
);
2954 long yy
= SCM_INUM (y
);
2955 return SCM_BOOL (xx
== yy
);
2957 else if (SCM_BIGP (y
))
2959 else if (SCM_REALP (y
))
2960 return SCM_BOOL ((double) xx
== SCM_REAL_VALUE (y
));
2961 else if (SCM_COMPLEXP (y
))
2962 return SCM_BOOL (((double) xx
== SCM_COMPLEX_REAL (y
))
2963 && (0.0 == SCM_COMPLEX_IMAG (y
)));
2964 else if (SCM_FRACTIONP (y
))
2967 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2969 else if (SCM_BIGP (x
))
2973 else if (SCM_BIGP (y
))
2975 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2976 scm_remember_upto_here_2 (x
, y
);
2977 return SCM_BOOL (0 == cmp
);
2979 else if (SCM_REALP (y
))
2982 if (xisnan (SCM_REAL_VALUE (y
)))
2984 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
2985 scm_remember_upto_here_1 (x
);
2986 return SCM_BOOL (0 == cmp
);
2988 else if (SCM_COMPLEXP (y
))
2991 if (0.0 != SCM_COMPLEX_IMAG (y
))
2993 if (xisnan (SCM_COMPLEX_REAL (y
)))
2995 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
2996 scm_remember_upto_here_1 (x
);
2997 return SCM_BOOL (0 == cmp
);
2999 else if (SCM_FRACTIONP (y
))
3002 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3004 else if (SCM_REALP (x
))
3007 return SCM_BOOL (SCM_REAL_VALUE (x
) == (double) SCM_INUM (y
));
3008 else if (SCM_BIGP (y
))
3011 if (xisnan (SCM_REAL_VALUE (x
)))
3013 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3014 scm_remember_upto_here_1 (y
);
3015 return SCM_BOOL (0 == cmp
);
3017 else if (SCM_REALP (y
))
3018 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3019 else if (SCM_COMPLEXP (y
))
3020 return SCM_BOOL ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3021 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3022 else if (SCM_FRACTIONP (y
))
3024 double xx
= SCM_REAL_VALUE (x
);
3028 return SCM_BOOL (xx
< 0.0);
3029 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3033 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3035 else if (SCM_COMPLEXP (x
))
3038 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == (double) SCM_INUM (y
))
3039 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3040 else if (SCM_BIGP (y
))
3043 if (0.0 != SCM_COMPLEX_IMAG (x
))
3045 if (xisnan (SCM_COMPLEX_REAL (x
)))
3047 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3048 scm_remember_upto_here_1 (y
);
3049 return SCM_BOOL (0 == cmp
);
3051 else if (SCM_REALP (y
))
3052 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3053 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3054 else if (SCM_COMPLEXP (y
))
3055 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3056 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3057 else if (SCM_FRACTIONP (y
))
3060 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3062 xx
= SCM_COMPLEX_REAL (x
);
3066 return SCM_BOOL (xx
< 0.0);
3067 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3071 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3073 else if (SCM_FRACTIONP (x
))
3077 else if (SCM_BIGP (y
))
3079 else if (SCM_REALP (y
))
3081 double yy
= SCM_REAL_VALUE (y
);
3085 return SCM_BOOL (0.0 < yy
);
3086 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3089 else if (SCM_COMPLEXP (y
))
3092 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3094 yy
= SCM_COMPLEX_REAL (y
);
3098 return SCM_BOOL (0.0 < yy
);
3099 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3102 else if (SCM_FRACTIONP (y
))
3103 return scm_i_fraction_equalp (x
, y
);
3105 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3108 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3112 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3113 done are good for inums, but for bignums an answer can almost always be
3114 had by just examining a few high bits of the operands, as done by GMP in
3115 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3116 of the float exponent to take into account. */
3118 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3119 /* "Return @code{#t} if the list of parameters is monotonically\n"
3123 scm_less_p (SCM x
, SCM y
)
3128 long xx
= SCM_INUM (x
);
3131 long yy
= SCM_INUM (y
);
3132 return SCM_BOOL (xx
< yy
);
3134 else if (SCM_BIGP (y
))
3136 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3137 scm_remember_upto_here_1 (y
);
3138 return SCM_BOOL (sgn
> 0);
3140 else if (SCM_REALP (y
))
3141 return SCM_BOOL ((double) xx
< SCM_REAL_VALUE (y
));
3142 else if (SCM_FRACTIONP (y
))
3144 /* "x < a/b" becomes "x*b < a" */
3146 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3147 y
= SCM_FRACTION_NUMERATOR (y
);
3151 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3153 else if (SCM_BIGP (x
))
3157 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3158 scm_remember_upto_here_1 (x
);
3159 return SCM_BOOL (sgn
< 0);
3161 else if (SCM_BIGP (y
))
3163 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3164 scm_remember_upto_here_2 (x
, y
);
3165 return SCM_BOOL (cmp
< 0);
3167 else if (SCM_REALP (y
))
3170 if (xisnan (SCM_REAL_VALUE (y
)))
3172 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3173 scm_remember_upto_here_1 (x
);
3174 return SCM_BOOL (cmp
< 0);
3176 else if (SCM_FRACTIONP (y
))
3179 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3181 else if (SCM_REALP (x
))
3184 return SCM_BOOL (SCM_REAL_VALUE (x
) < (double) SCM_INUM (y
));
3185 else if (SCM_BIGP (y
))
3188 if (xisnan (SCM_REAL_VALUE (x
)))
3190 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3191 scm_remember_upto_here_1 (y
);
3192 return SCM_BOOL (cmp
> 0);
3194 else if (SCM_REALP (y
))
3195 return SCM_BOOL (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3196 else if (SCM_FRACTIONP (y
))
3198 double xx
= SCM_REAL_VALUE (x
);
3202 return SCM_BOOL (xx
< 0.0);
3203 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3207 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3209 else if (SCM_FRACTIONP (x
))
3211 if (SCM_INUMP (y
) || SCM_BIGP (y
))
3213 /* "a/b < y" becomes "a < y*b" */
3214 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3215 x
= SCM_FRACTION_NUMERATOR (x
);
3218 else if (SCM_REALP (y
))
3220 double yy
= SCM_REAL_VALUE (y
);
3224 return SCM_BOOL (0.0 < yy
);
3225 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3228 else if (SCM_FRACTIONP (y
))
3230 /* "a/b < c/d" becomes "a*d < c*b" */
3231 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3232 SCM_FRACTION_DENOMINATOR (y
));
3233 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3234 SCM_FRACTION_DENOMINATOR (x
));
3240 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3243 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3247 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3248 /* "Return @code{#t} if the list of parameters is monotonically\n"
3251 #define FUNC_NAME s_scm_gr_p
3253 scm_gr_p (SCM x
, SCM y
)
3255 if (!SCM_NUMBERP (x
))
3256 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3257 else if (!SCM_NUMBERP (y
))
3258 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3260 return scm_less_p (y
, x
);
3265 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3266 /* "Return @code{#t} if the list of parameters is monotonically\n"
3269 #define FUNC_NAME s_scm_leq_p
3271 scm_leq_p (SCM x
, SCM y
)
3273 if (!SCM_NUMBERP (x
))
3274 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3275 else if (!SCM_NUMBERP (y
))
3276 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3277 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3280 return SCM_BOOL_NOT (scm_less_p (y
, x
));
3285 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3286 /* "Return @code{#t} if the list of parameters is monotonically\n"
3289 #define FUNC_NAME s_scm_geq_p
3291 scm_geq_p (SCM x
, SCM y
)
3293 if (!SCM_NUMBERP (x
))
3294 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3295 else if (!SCM_NUMBERP (y
))
3296 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3297 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3300 return SCM_BOOL_NOT (scm_less_p (x
, y
));
3305 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3306 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3313 return SCM_BOOL (SCM_EQ_P (z
, SCM_INUM0
));
3314 else if (SCM_BIGP (z
))
3316 else if (SCM_REALP (z
))
3317 return SCM_BOOL (SCM_REAL_VALUE (z
) == 0.0);
3318 else if (SCM_COMPLEXP (z
))
3319 return SCM_BOOL (SCM_COMPLEX_REAL (z
) == 0.0
3320 && SCM_COMPLEX_IMAG (z
) == 0.0);
3321 else if (SCM_FRACTIONP (z
))
3324 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3328 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3329 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3333 scm_positive_p (SCM x
)
3336 return SCM_BOOL (SCM_INUM (x
) > 0);
3337 else if (SCM_BIGP (x
))
3339 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3340 scm_remember_upto_here_1 (x
);
3341 return SCM_BOOL (sgn
> 0);
3343 else if (SCM_REALP (x
))
3344 return SCM_BOOL(SCM_REAL_VALUE (x
) > 0.0);
3345 else if (SCM_FRACTIONP (x
))
3346 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3348 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3352 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3353 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3357 scm_negative_p (SCM x
)
3360 return SCM_BOOL (SCM_INUM (x
) < 0);
3361 else if (SCM_BIGP (x
))
3363 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3364 scm_remember_upto_here_1 (x
);
3365 return SCM_BOOL (sgn
< 0);
3367 else if (SCM_REALP (x
))
3368 return SCM_BOOL(SCM_REAL_VALUE (x
) < 0.0);
3369 else if (SCM_FRACTIONP (x
))
3370 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3372 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3376 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3377 /* "Return the maximum of all parameter values."
3380 scm_max (SCM x
, SCM y
)
3385 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3386 else if (SCM_NUMBERP (x
))
3389 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3394 long xx
= SCM_INUM (x
);
3397 long yy
= SCM_INUM (y
);
3398 return (xx
< yy
) ? y
: x
;
3400 else if (SCM_BIGP (y
))
3402 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3403 scm_remember_upto_here_1 (y
);
3404 return (sgn
< 0) ? x
: y
;
3406 else if (SCM_REALP (y
))
3409 /* if y==NaN then ">" is false and we return NaN */
3410 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3412 else if (SCM_FRACTIONP (y
))
3415 return (z
> scm_i_fraction2double (y
)) ? x
: y
;
3418 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3420 else if (SCM_BIGP (x
))
3424 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3425 scm_remember_upto_here_1 (x
);
3426 return (sgn
< 0) ? y
: x
;
3428 else if (SCM_BIGP (y
))
3430 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3431 scm_remember_upto_here_2 (x
, y
);
3432 return (cmp
> 0) ? x
: y
;
3434 else if (SCM_REALP (y
))
3436 double yy
= SCM_REAL_VALUE (y
);
3440 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3441 scm_remember_upto_here_1 (x
);
3442 return (cmp
> 0) ? x
: y
;
3444 else if (SCM_FRACTIONP (y
))
3446 double yy
= scm_i_fraction2double (y
);
3448 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3449 scm_remember_upto_here_1 (x
);
3450 return (cmp
> 0) ? x
: y
;
3453 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3455 else if (SCM_REALP (x
))
3459 double z
= SCM_INUM (y
);
3460 /* if x==NaN then "<" is false and we return NaN */
3461 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3463 else if (SCM_BIGP (y
))
3465 double xx
= SCM_REAL_VALUE (x
);
3469 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3470 scm_remember_upto_here_1 (y
);
3471 return (cmp
< 0) ? x
: y
;
3473 else if (SCM_REALP (y
))
3475 /* if x==NaN then our explicit check means we return NaN
3476 if y==NaN then ">" is false and we return NaN
3477 calling isnan is unavoidable, since it's the only way to know
3478 which of x or y causes any compares to be false */
3479 double xx
= SCM_REAL_VALUE (x
);
3480 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3482 else if (SCM_FRACTIONP (y
))
3484 double yy
= scm_i_fraction2double (y
);
3485 double xx
= SCM_REAL_VALUE (x
);
3486 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3489 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3491 else if (SCM_FRACTIONP (x
))
3495 double z
= SCM_INUM (y
);
3496 return (scm_i_fraction2double (x
) < z
) ? y
: x
;
3498 else if (SCM_BIGP (y
))
3500 double xx
= scm_i_fraction2double (x
);
3502 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3503 scm_remember_upto_here_1 (y
);
3504 return (cmp
< 0) ? x
: y
;
3506 else if (SCM_REALP (y
))
3508 double xx
= scm_i_fraction2double (x
);
3509 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3511 else if (SCM_FRACTIONP (y
))
3513 double yy
= scm_i_fraction2double (y
);
3514 double xx
= scm_i_fraction2double (x
);
3515 return (xx
< yy
) ? y
: x
;
3518 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3521 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3525 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3526 /* "Return the minium of all parameter values."
3529 scm_min (SCM x
, SCM y
)
3534 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3535 else if (SCM_NUMBERP (x
))
3538 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3543 long xx
= SCM_INUM (x
);
3546 long yy
= SCM_INUM (y
);
3547 return (xx
< yy
) ? x
: y
;
3549 else if (SCM_BIGP (y
))
3551 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3552 scm_remember_upto_here_1 (y
);
3553 return (sgn
< 0) ? y
: x
;
3555 else if (SCM_REALP (y
))
3558 /* if y==NaN then "<" is false and we return NaN */
3559 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3561 else if (SCM_FRACTIONP (y
))
3564 return (z
< scm_i_fraction2double (y
)) ? x
: y
;
3567 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3569 else if (SCM_BIGP (x
))
3573 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3574 scm_remember_upto_here_1 (x
);
3575 return (sgn
< 0) ? x
: y
;
3577 else if (SCM_BIGP (y
))
3579 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3580 scm_remember_upto_here_2 (x
, y
);
3581 return (cmp
> 0) ? y
: x
;
3583 else if (SCM_REALP (y
))
3585 double yy
= SCM_REAL_VALUE (y
);
3589 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3590 scm_remember_upto_here_1 (x
);
3591 return (cmp
> 0) ? y
: x
;
3593 else if (SCM_FRACTIONP (y
))
3595 double yy
= scm_i_fraction2double (y
);
3597 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3598 scm_remember_upto_here_1 (x
);
3599 return (cmp
> 0) ? y
: x
;
3602 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3604 else if (SCM_REALP (x
))
3608 double z
= SCM_INUM (y
);
3609 /* if x==NaN then "<" is false and we return NaN */
3610 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3612 else if (SCM_BIGP (y
))
3614 double xx
= SCM_REAL_VALUE (x
);
3618 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3619 scm_remember_upto_here_1 (y
);
3620 return (cmp
< 0) ? y
: x
;
3622 else if (SCM_REALP (y
))
3624 /* if x==NaN then our explicit check means we return NaN
3625 if y==NaN then "<" is false and we return NaN
3626 calling isnan is unavoidable, since it's the only way to know
3627 which of x or y causes any compares to be false */
3628 double xx
= SCM_REAL_VALUE (x
);
3629 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3631 else if (SCM_FRACTIONP (y
))
3633 double yy
= scm_i_fraction2double (y
);
3634 double xx
= SCM_REAL_VALUE (x
);
3635 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3638 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3640 else if (SCM_FRACTIONP (x
))
3644 double z
= SCM_INUM (y
);
3645 return (scm_i_fraction2double (x
) < z
) ? x
: y
;
3647 else if (SCM_BIGP (y
))
3649 double xx
= scm_i_fraction2double (x
);
3651 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3652 scm_remember_upto_here_1 (y
);
3653 return (cmp
< 0) ? y
: x
;
3655 else if (SCM_REALP (y
))
3657 double xx
= scm_i_fraction2double (x
);
3658 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3660 else if (SCM_FRACTIONP (y
))
3662 double yy
= scm_i_fraction2double (y
);
3663 double xx
= scm_i_fraction2double (x
);
3664 return (xx
< yy
) ? x
: y
;
3667 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3670 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3674 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3675 /* "Return the sum of all parameter values. Return 0 if called without\n"
3679 scm_sum (SCM x
, SCM y
)
3683 if (SCM_NUMBERP (x
)) return x
;
3684 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3685 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3692 long xx
= SCM_INUM (x
);
3693 long yy
= SCM_INUM (y
);
3694 long int z
= xx
+ yy
;
3695 return SCM_FIXABLE (z
) ? SCM_MAKINUM (z
) : scm_i_long2big (z
);
3697 else if (SCM_BIGP (y
))
3702 else if (SCM_REALP (y
))
3704 long int xx
= SCM_INUM (x
);
3705 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3707 else if (SCM_COMPLEXP (y
))
3709 long int xx
= SCM_INUM (x
);
3710 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3711 SCM_COMPLEX_IMAG (y
));
3713 else if (SCM_FRACTIONP (y
))
3714 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3715 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3716 SCM_FRACTION_DENOMINATOR (y
));
3718 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3719 } else if (SCM_BIGP (x
))
3726 inum
= SCM_INUM (y
);
3729 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3732 SCM result
= scm_i_mkbig ();
3733 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3734 scm_remember_upto_here_1 (x
);
3735 /* we know the result will have to be a bignum */
3738 return scm_i_normbig (result
);
3742 SCM result
= scm_i_mkbig ();
3743 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3744 scm_remember_upto_here_1 (x
);
3745 /* we know the result will have to be a bignum */
3748 return scm_i_normbig (result
);
3751 else if (SCM_BIGP (y
))
3753 SCM result
= scm_i_mkbig ();
3754 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3755 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3756 mpz_add (SCM_I_BIG_MPZ (result
),
3759 scm_remember_upto_here_2 (x
, y
);
3760 /* we know the result will have to be a bignum */
3763 return scm_i_normbig (result
);
3765 else if (SCM_REALP (y
))
3767 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3768 scm_remember_upto_here_1 (x
);
3769 return scm_make_real (result
);
3771 else if (SCM_COMPLEXP (y
))
3773 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3774 + SCM_COMPLEX_REAL (y
));
3775 scm_remember_upto_here_1 (x
);
3776 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3778 else if (SCM_FRACTIONP (y
))
3779 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3780 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3781 SCM_FRACTION_DENOMINATOR (y
));
3783 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3785 else if (SCM_REALP (x
))
3788 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_INUM (y
));
3789 else if (SCM_BIGP (y
))
3791 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3792 scm_remember_upto_here_1 (y
);
3793 return scm_make_real (result
);
3795 else if (SCM_REALP (y
))
3796 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3797 else if (SCM_COMPLEXP (y
))
3798 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3799 SCM_COMPLEX_IMAG (y
));
3800 else if (SCM_FRACTIONP (y
))
3801 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3803 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3805 else if (SCM_COMPLEXP (x
))
3808 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_INUM (y
),
3809 SCM_COMPLEX_IMAG (x
));
3810 else if (SCM_BIGP (y
))
3812 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3813 + SCM_COMPLEX_REAL (x
));
3814 scm_remember_upto_here_1 (y
);
3815 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3817 else if (SCM_REALP (y
))
3818 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3819 SCM_COMPLEX_IMAG (x
));
3820 else if (SCM_COMPLEXP (y
))
3821 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3822 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3823 else if (SCM_FRACTIONP (y
))
3824 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3825 SCM_COMPLEX_IMAG (x
));
3827 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3829 else if (SCM_FRACTIONP (x
))
3832 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3833 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3834 SCM_FRACTION_DENOMINATOR (x
));
3835 else if (SCM_BIGP (y
))
3836 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3837 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3838 SCM_FRACTION_DENOMINATOR (x
));
3839 else if (SCM_REALP (y
))
3840 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3841 else if (SCM_COMPLEXP (y
))
3842 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3843 SCM_COMPLEX_IMAG (y
));
3844 else if (SCM_FRACTIONP (y
))
3845 /* a/b + c/d = (ad + bc) / bd */
3846 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3847 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3848 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3850 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3853 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
3857 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
3858 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
3859 * the sum of all but the first argument are subtracted from the first
3861 #define FUNC_NAME s_difference
3863 scm_difference (SCM x
, SCM y
)
3868 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
3872 long xx
= -SCM_INUM (x
);
3873 if (SCM_FIXABLE (xx
))
3874 return SCM_MAKINUM (xx
);
3876 return scm_i_long2big (xx
);
3878 else if (SCM_BIGP (x
))
3879 /* FIXME: do we really need to normalize here? */
3880 return scm_i_normbig (scm_i_clonebig (x
, 0));
3881 else if (SCM_REALP (x
))
3882 return scm_make_real (-SCM_REAL_VALUE (x
));
3883 else if (SCM_COMPLEXP (x
))
3884 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
3885 -SCM_COMPLEX_IMAG (x
));
3886 else if (SCM_FRACTIONP (x
))
3887 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
3888 SCM_FRACTION_DENOMINATOR (x
));
3890 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
3897 long int xx
= SCM_INUM (x
);
3898 long int yy
= SCM_INUM (y
);
3899 long int z
= xx
- yy
;
3900 if (SCM_FIXABLE (z
))
3901 return SCM_MAKINUM (z
);
3903 return scm_i_long2big (z
);
3905 else if (SCM_BIGP (y
))
3907 /* inum-x - big-y */
3908 long xx
= SCM_INUM (x
);
3911 return scm_i_clonebig (y
, 0);
3914 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3915 SCM result
= scm_i_mkbig ();
3918 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
3921 /* x - y == -(y + -x) */
3922 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
3923 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
3925 scm_remember_upto_here_1 (y
);
3927 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
3928 /* we know the result will have to be a bignum */
3931 return scm_i_normbig (result
);
3934 else if (SCM_REALP (y
))
3936 long int xx
= SCM_INUM (x
);
3937 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
3939 else if (SCM_COMPLEXP (y
))
3941 long int xx
= SCM_INUM (x
);
3942 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
3943 - SCM_COMPLEX_IMAG (y
));
3945 else if (SCM_FRACTIONP (y
))
3946 /* a - b/c = (ac - b) / c */
3947 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3948 SCM_FRACTION_NUMERATOR (y
)),
3949 SCM_FRACTION_DENOMINATOR (y
));
3951 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3953 else if (SCM_BIGP (x
))
3957 /* big-x - inum-y */
3958 long yy
= SCM_INUM (y
);
3959 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3961 scm_remember_upto_here_1 (x
);
3963 return SCM_FIXABLE (-yy
) ? SCM_MAKINUM (-yy
) : scm_long2num (-yy
);
3966 SCM result
= scm_i_mkbig ();
3969 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
3971 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
3972 scm_remember_upto_here_1 (x
);
3974 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
3975 /* we know the result will have to be a bignum */
3978 return scm_i_normbig (result
);
3981 else if (SCM_BIGP (y
))
3983 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3984 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3985 SCM result
= scm_i_mkbig ();
3986 mpz_sub (SCM_I_BIG_MPZ (result
),
3989 scm_remember_upto_here_2 (x
, y
);
3990 /* we know the result will have to be a bignum */
3991 if ((sgn_x
== 1) && (sgn_y
== -1))
3993 if ((sgn_x
== -1) && (sgn_y
== 1))
3995 return scm_i_normbig (result
);
3997 else if (SCM_REALP (y
))
3999 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4000 scm_remember_upto_here_1 (x
);
4001 return scm_make_real (result
);
4003 else if (SCM_COMPLEXP (y
))
4005 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4006 - SCM_COMPLEX_REAL (y
));
4007 scm_remember_upto_here_1 (x
);
4008 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
4010 else if (SCM_FRACTIONP (y
))
4011 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4012 SCM_FRACTION_NUMERATOR (y
)),
4013 SCM_FRACTION_DENOMINATOR (y
));
4014 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4016 else if (SCM_REALP (x
))
4019 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_INUM (y
));
4020 else if (SCM_BIGP (y
))
4022 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4023 scm_remember_upto_here_1 (x
);
4024 return scm_make_real (result
);
4026 else if (SCM_REALP (y
))
4027 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4028 else if (SCM_COMPLEXP (y
))
4029 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4030 -SCM_COMPLEX_IMAG (y
));
4031 else if (SCM_FRACTIONP (y
))
4032 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4034 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4036 else if (SCM_COMPLEXP (x
))
4039 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_INUM (y
),
4040 SCM_COMPLEX_IMAG (x
));
4041 else if (SCM_BIGP (y
))
4043 double real_part
= (SCM_COMPLEX_REAL (x
)
4044 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4045 scm_remember_upto_here_1 (x
);
4046 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
4048 else if (SCM_REALP (y
))
4049 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4050 SCM_COMPLEX_IMAG (x
));
4051 else if (SCM_COMPLEXP (y
))
4052 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4053 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4054 else if (SCM_FRACTIONP (y
))
4055 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4056 SCM_COMPLEX_IMAG (x
));
4058 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4060 else if (SCM_FRACTIONP (x
))
4063 /* a/b - c = (a - cb) / b */
4064 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4065 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4066 SCM_FRACTION_DENOMINATOR (x
));
4067 else if (SCM_BIGP (y
))
4068 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4069 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4070 SCM_FRACTION_DENOMINATOR (x
));
4071 else if (SCM_REALP (y
))
4072 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4073 else if (SCM_COMPLEXP (y
))
4074 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4075 -SCM_COMPLEX_IMAG (y
));
4076 else if (SCM_FRACTIONP (y
))
4077 /* a/b - c/d = (ad - bc) / bd */
4078 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4079 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4080 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4082 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4085 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4090 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4091 /* "Return the product of all arguments. If called without arguments,\n"
4095 scm_product (SCM x
, SCM y
)
4100 return SCM_MAKINUM (1L);
4101 else if (SCM_NUMBERP (x
))
4104 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4116 case 0: return x
; break;
4117 case 1: return y
; break;
4122 long yy
= SCM_INUM (y
);
4124 SCM k
= SCM_MAKINUM (kk
);
4125 if ((kk
== SCM_INUM (k
)) && (kk
/ xx
== yy
))
4129 SCM result
= scm_i_long2big (xx
);
4130 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4131 return scm_i_normbig (result
);
4134 else if (SCM_BIGP (y
))
4136 SCM result
= scm_i_mkbig ();
4137 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4138 scm_remember_upto_here_1 (y
);
4141 else if (SCM_REALP (y
))
4142 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4143 else if (SCM_COMPLEXP (y
))
4144 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4145 xx
* SCM_COMPLEX_IMAG (y
));
4146 else if (SCM_FRACTIONP (y
))
4147 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4148 SCM_FRACTION_DENOMINATOR (y
));
4150 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4152 else if (SCM_BIGP (x
))
4159 else if (SCM_BIGP (y
))
4161 SCM result
= scm_i_mkbig ();
4162 mpz_mul (SCM_I_BIG_MPZ (result
),
4165 scm_remember_upto_here_2 (x
, y
);
4168 else if (SCM_REALP (y
))
4170 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4171 scm_remember_upto_here_1 (x
);
4172 return scm_make_real (result
);
4174 else if (SCM_COMPLEXP (y
))
4176 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4177 scm_remember_upto_here_1 (x
);
4178 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4179 z
* SCM_COMPLEX_IMAG (y
));
4181 else if (SCM_FRACTIONP (y
))
4182 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4183 SCM_FRACTION_DENOMINATOR (y
));
4185 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4187 else if (SCM_REALP (x
))
4190 return scm_make_real (SCM_INUM (y
) * SCM_REAL_VALUE (x
));
4191 else if (SCM_BIGP (y
))
4193 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4194 scm_remember_upto_here_1 (y
);
4195 return scm_make_real (result
);
4197 else if (SCM_REALP (y
))
4198 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4199 else if (SCM_COMPLEXP (y
))
4200 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4201 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4202 else if (SCM_FRACTIONP (y
))
4203 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4205 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4207 else if (SCM_COMPLEXP (x
))
4210 return scm_make_complex (SCM_INUM (y
) * SCM_COMPLEX_REAL (x
),
4211 SCM_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4212 else if (SCM_BIGP (y
))
4214 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4215 scm_remember_upto_here_1 (y
);
4216 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4217 z
* SCM_COMPLEX_IMAG (x
));
4219 else if (SCM_REALP (y
))
4220 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4221 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4222 else if (SCM_COMPLEXP (y
))
4224 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4225 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4226 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4227 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4229 else if (SCM_FRACTIONP (y
))
4231 double yy
= scm_i_fraction2double (y
);
4232 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4233 yy
* SCM_COMPLEX_IMAG (x
));
4236 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4238 else if (SCM_FRACTIONP (x
))
4241 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4242 SCM_FRACTION_DENOMINATOR (x
));
4243 else if (SCM_BIGP (y
))
4244 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4245 SCM_FRACTION_DENOMINATOR (x
));
4246 else if (SCM_REALP (y
))
4247 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4248 else if (SCM_COMPLEXP (y
))
4250 double xx
= scm_i_fraction2double (x
);
4251 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4252 xx
* SCM_COMPLEX_IMAG (y
));
4254 else if (SCM_FRACTIONP (y
))
4255 /* a/b * c/d = ac / bd */
4256 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4257 SCM_FRACTION_NUMERATOR (y
)),
4258 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4259 SCM_FRACTION_DENOMINATOR (y
)));
4261 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4264 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4268 scm_num2dbl (SCM a
, const char *why
)
4269 #define FUNC_NAME why
4272 return (double) SCM_INUM (a
);
4273 else if (SCM_BIGP (a
))
4275 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4276 scm_remember_upto_here_1 (a
);
4279 else if (SCM_REALP (a
))
4280 return (SCM_REAL_VALUE (a
));
4281 else if (SCM_FRACTIONP (a
))
4282 return scm_i_fraction2double (a
);
4284 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4288 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4289 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4290 #define ALLOW_DIVIDE_BY_ZERO
4291 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4294 /* The code below for complex division is adapted from the GNU
4295 libstdc++, which adapted it from f2c's libF77, and is subject to
4298 /****************************************************************
4299 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4301 Permission to use, copy, modify, and distribute this software
4302 and its documentation for any purpose and without fee is hereby
4303 granted, provided that the above copyright notice appear in all
4304 copies and that both that the copyright notice and this
4305 permission notice and warranty disclaimer appear in supporting
4306 documentation, and that the names of AT&T Bell Laboratories or
4307 Bellcore or any of their entities not be used in advertising or
4308 publicity pertaining to distribution of the software without
4309 specific, written prior permission.
4311 AT&T and Bellcore disclaim all warranties with regard to this
4312 software, including all implied warranties of merchantability
4313 and fitness. In no event shall AT&T or Bellcore be liable for
4314 any special, indirect or consequential damages or any damages
4315 whatsoever resulting from loss of use, data or profits, whether
4316 in an action of contract, negligence or other tortious action,
4317 arising out of or in connection with the use or performance of
4319 ****************************************************************/
4321 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4322 /* Divide the first argument by the product of the remaining
4323 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4325 #define FUNC_NAME s_divide
4327 scm_i_divide (SCM x
, SCM y
, int inexact
)
4334 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4335 else if (SCM_INUMP (x
))
4337 long xx
= SCM_INUM (x
);
4338 if (xx
== 1 || xx
== -1)
4340 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4342 scm_num_overflow (s_divide
);
4347 return scm_make_real (1.0 / (double) xx
);
4348 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4351 else if (SCM_BIGP (x
))
4354 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4355 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4357 else if (SCM_REALP (x
))
4359 double xx
= SCM_REAL_VALUE (x
);
4360 #ifndef ALLOW_DIVIDE_BY_ZERO
4362 scm_num_overflow (s_divide
);
4365 return scm_make_real (1.0 / xx
);
4367 else if (SCM_COMPLEXP (x
))
4369 double r
= SCM_COMPLEX_REAL (x
);
4370 double i
= SCM_COMPLEX_IMAG (x
);
4374 double d
= i
* (1.0 + t
* t
);
4375 return scm_make_complex (t
/ d
, -1.0 / d
);
4380 double d
= r
* (1.0 + t
* t
);
4381 return scm_make_complex (1.0 / d
, -t
/ d
);
4384 else if (SCM_FRACTIONP (x
))
4385 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4386 SCM_FRACTION_NUMERATOR (x
));
4388 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4393 long xx
= SCM_INUM (x
);
4396 long yy
= SCM_INUM (y
);
4399 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4400 scm_num_overflow (s_divide
);
4402 return scm_make_real ((double) xx
/ (double) yy
);
4405 else if (xx
% yy
!= 0)
4408 return scm_make_real ((double) xx
/ (double) yy
);
4409 else return scm_make_ratio (x
, y
);
4414 if (SCM_FIXABLE (z
))
4415 return SCM_MAKINUM (z
);
4417 return scm_i_long2big (z
);
4420 else if (SCM_BIGP (y
))
4423 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4424 else return scm_make_ratio (x
, y
);
4426 else if (SCM_REALP (y
))
4428 double yy
= SCM_REAL_VALUE (y
);
4429 #ifndef ALLOW_DIVIDE_BY_ZERO
4431 scm_num_overflow (s_divide
);
4434 return scm_make_real ((double) xx
/ yy
);
4436 else if (SCM_COMPLEXP (y
))
4439 complex_div
: /* y _must_ be a complex number */
4441 double r
= SCM_COMPLEX_REAL (y
);
4442 double i
= SCM_COMPLEX_IMAG (y
);
4446 double d
= i
* (1.0 + t
* t
);
4447 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4452 double d
= r
* (1.0 + t
* t
);
4453 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4457 else if (SCM_FRACTIONP (y
))
4458 /* a / b/c = ac / b */
4459 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4460 SCM_FRACTION_NUMERATOR (y
));
4462 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4464 else if (SCM_BIGP (x
))
4468 long int yy
= SCM_INUM (y
);
4471 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4472 scm_num_overflow (s_divide
);
4474 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4475 scm_remember_upto_here_1 (x
);
4476 return (sgn
== 0) ? scm_nan () : scm_inf ();
4483 /* FIXME: HMM, what are the relative performance issues here?
4484 We need to test. Is it faster on average to test
4485 divisible_p, then perform whichever operation, or is it
4486 faster to perform the integer div opportunistically and
4487 switch to real if there's a remainder? For now we take the
4488 middle ground: test, then if divisible, use the faster div
4491 long abs_yy
= yy
< 0 ? -yy
: yy
;
4492 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4496 SCM result
= scm_i_mkbig ();
4497 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4498 scm_remember_upto_here_1 (x
);
4500 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4501 return scm_i_normbig (result
);
4506 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4507 else return scm_make_ratio (x
, y
);
4511 else if (SCM_BIGP (y
))
4513 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4516 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4517 scm_num_overflow (s_divide
);
4519 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4520 scm_remember_upto_here_1 (x
);
4521 return (sgn
== 0) ? scm_nan () : scm_inf ();
4527 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4531 SCM result
= scm_i_mkbig ();
4532 mpz_divexact (SCM_I_BIG_MPZ (result
),
4535 scm_remember_upto_here_2 (x
, y
);
4536 return scm_i_normbig (result
);
4542 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4543 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4544 scm_remember_upto_here_2 (x
, y
);
4545 return scm_make_real (dbx
/ dby
);
4547 else return scm_make_ratio (x
, y
);
4551 else if (SCM_REALP (y
))
4553 double yy
= SCM_REAL_VALUE (y
);
4554 #ifndef ALLOW_DIVIDE_BY_ZERO
4556 scm_num_overflow (s_divide
);
4559 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4561 else if (SCM_COMPLEXP (y
))
4563 a
= scm_i_big2dbl (x
);
4566 else if (SCM_FRACTIONP (y
))
4567 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4568 SCM_FRACTION_NUMERATOR (y
));
4570 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4572 else if (SCM_REALP (x
))
4574 double rx
= SCM_REAL_VALUE (x
);
4577 long int yy
= SCM_INUM (y
);
4578 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4580 scm_num_overflow (s_divide
);
4583 return scm_make_real (rx
/ (double) yy
);
4585 else if (SCM_BIGP (y
))
4587 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4588 scm_remember_upto_here_1 (y
);
4589 return scm_make_real (rx
/ dby
);
4591 else if (SCM_REALP (y
))
4593 double yy
= SCM_REAL_VALUE (y
);
4594 #ifndef ALLOW_DIVIDE_BY_ZERO
4596 scm_num_overflow (s_divide
);
4599 return scm_make_real (rx
/ yy
);
4601 else if (SCM_COMPLEXP (y
))
4606 else if (SCM_FRACTIONP (y
))
4607 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4609 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4611 else if (SCM_COMPLEXP (x
))
4613 double rx
= SCM_COMPLEX_REAL (x
);
4614 double ix
= SCM_COMPLEX_IMAG (x
);
4617 long int yy
= SCM_INUM (y
);
4618 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4620 scm_num_overflow (s_divide
);
4625 return scm_make_complex (rx
/ d
, ix
/ d
);
4628 else if (SCM_BIGP (y
))
4630 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4631 scm_remember_upto_here_1 (y
);
4632 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4634 else if (SCM_REALP (y
))
4636 double yy
= SCM_REAL_VALUE (y
);
4637 #ifndef ALLOW_DIVIDE_BY_ZERO
4639 scm_num_overflow (s_divide
);
4642 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4644 else if (SCM_COMPLEXP (y
))
4646 double ry
= SCM_COMPLEX_REAL (y
);
4647 double iy
= SCM_COMPLEX_IMAG (y
);
4651 double d
= iy
* (1.0 + t
* t
);
4652 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4657 double d
= ry
* (1.0 + t
* t
);
4658 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4661 else if (SCM_FRACTIONP (y
))
4663 double yy
= scm_i_fraction2double (y
);
4664 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4667 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4669 else if (SCM_FRACTIONP (x
))
4673 long int yy
= SCM_INUM (y
);
4674 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4676 scm_num_overflow (s_divide
);
4679 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4680 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4682 else if (SCM_BIGP (y
))
4684 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4685 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4687 else if (SCM_REALP (y
))
4689 double yy
= SCM_REAL_VALUE (y
);
4690 #ifndef ALLOW_DIVIDE_BY_ZERO
4692 scm_num_overflow (s_divide
);
4695 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4697 else if (SCM_COMPLEXP (y
))
4699 a
= scm_i_fraction2double (x
);
4702 else if (SCM_FRACTIONP (y
))
4703 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4704 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4706 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4709 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4713 scm_divide (SCM x
, SCM y
)
4715 return scm_i_divide (x
, y
, 0);
4718 static SCM
scm_divide2real (SCM x
, SCM y
)
4720 return scm_i_divide (x
, y
, 1);
4726 scm_asinh (double x
)
4731 #define asinh scm_asinh
4732 return log (x
+ sqrt (x
* x
+ 1));
4735 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4736 /* "Return the inverse hyperbolic sine of @var{x}."
4741 scm_acosh (double x
)
4746 #define acosh scm_acosh
4747 return log (x
+ sqrt (x
* x
- 1));
4750 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4751 /* "Return the inverse hyperbolic cosine of @var{x}."
4756 scm_atanh (double x
)
4761 #define atanh scm_atanh
4762 return 0.5 * log ((1 + x
) / (1 - x
));
4765 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4766 /* "Return the inverse hyperbolic tangent of @var{x}."
4770 /* XXX - eventually, we should remove this definition of scm_round and
4771 rename scm_round_number to scm_round. Likewise for scm_truncate
4772 and scm_truncate_number.
4776 scm_truncate (double x
)
4781 #define trunc scm_truncate
4789 scm_round (double x
)
4791 double plus_half
= x
+ 0.5;
4792 double result
= floor (plus_half
);
4793 /* Adjust so that the scm_round is towards even. */
4794 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4799 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4801 "Round the number @var{x} towards zero.")
4802 #define FUNC_NAME s_scm_truncate_number
4804 if (SCM_FALSEP (scm_negative_p (x
)))
4805 return scm_floor (x
);
4807 return scm_ceiling (x
);
4811 static SCM exactly_one_half
;
4813 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4815 "Round the number @var{x} towards the nearest integer. "
4816 "When it is exactly halfway between two integers, "
4817 "round towards the even one.")
4818 #define FUNC_NAME s_scm_round_number
4820 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4821 SCM result
= scm_floor (plus_half
);
4822 /* Adjust so that the scm_round is towards even. */
4823 if (!SCM_FALSEP (scm_num_eq_p (plus_half
, result
))
4824 && !SCM_FALSEP (scm_odd_p (result
)))
4825 return scm_difference (result
, SCM_MAKINUM (1));
4831 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4833 "Round the number @var{x} towards minus infinity.")
4834 #define FUNC_NAME s_scm_floor
4836 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4838 else if (SCM_REALP (x
))
4839 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
4840 else if (SCM_FRACTIONP (x
))
4842 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4843 SCM_FRACTION_DENOMINATOR (x
));
4844 if (SCM_FALSEP (scm_negative_p (x
)))
4846 /* For positive x, rounding towards zero is correct. */
4851 /* For negative x, we need to return q-1 unless x is an
4852 integer. But fractions are never integer, per our
4854 return scm_difference (q
, SCM_MAKINUM (1));
4858 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
4862 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
4864 "Round the number @var{x} towards infinity.")
4865 #define FUNC_NAME s_scm_ceiling
4867 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4869 else if (SCM_REALP (x
))
4870 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
4871 else if (SCM_FRACTIONP (x
))
4873 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4874 SCM_FRACTION_DENOMINATOR (x
));
4875 if (SCM_FALSEP (scm_positive_p (x
)))
4877 /* For negative x, rounding towards zero is correct. */
4882 /* For positive x, we need to return q+1 unless x is an
4883 integer. But fractions are never integer, per our
4885 return scm_sum (q
, SCM_MAKINUM (1));
4889 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
4893 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
4894 /* "Return the square root of the real number @var{x}."
4896 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
4897 /* "Return the absolute value of the real number @var{x}."
4899 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
4900 /* "Return the @var{x}th power of e."
4902 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
4903 /* "Return the natural logarithm of the real number @var{x}."
4905 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
4906 /* "Return the sine of the real number @var{x}."
4908 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
4909 /* "Return the cosine of the real number @var{x}."
4911 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
4912 /* "Return the tangent of the real number @var{x}."
4914 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
4915 /* "Return the arc sine of the real number @var{x}."
4917 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
4918 /* "Return the arc cosine of the real number @var{x}."
4920 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
4921 /* "Return the arc tangent of the real number @var{x}."
4923 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
4924 /* "Return the hyperbolic sine of the real number @var{x}."
4926 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
4927 /* "Return the hyperbolic cosine of the real number @var{x}."
4929 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
4930 /* "Return the hyperbolic tangent of the real number @var{x}."
4938 static void scm_two_doubles (SCM x
,
4940 const char *sstring
,
4944 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
4947 xy
->x
= SCM_INUM (x
);
4948 else if (SCM_BIGP (x
))
4949 xy
->x
= scm_i_big2dbl (x
);
4950 else if (SCM_REALP (x
))
4951 xy
->x
= SCM_REAL_VALUE (x
);
4952 else if (SCM_FRACTIONP (x
))
4953 xy
->x
= scm_i_fraction2double (x
);
4955 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
4958 xy
->y
= SCM_INUM (y
);
4959 else if (SCM_BIGP (y
))
4960 xy
->y
= scm_i_big2dbl (y
);
4961 else if (SCM_REALP (y
))
4962 xy
->y
= SCM_REAL_VALUE (y
);
4963 else if (SCM_FRACTIONP (y
))
4964 xy
->y
= scm_i_fraction2double (y
);
4966 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
4970 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
4972 "Return @var{x} raised to the power of @var{y}. This\n"
4973 "procedure does not accept complex arguments.")
4974 #define FUNC_NAME s_scm_sys_expt
4977 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4978 return scm_make_real (pow (xy
.x
, xy
.y
));
4983 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
4985 "Return the arc tangent of the two arguments @var{x} and\n"
4986 "@var{y}. This is similar to calculating the arc tangent of\n"
4987 "@var{x} / @var{y}, except that the signs of both arguments\n"
4988 "are used to determine the quadrant of the result. This\n"
4989 "procedure does not accept complex arguments.")
4990 #define FUNC_NAME s_scm_sys_atan2
4993 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4994 return scm_make_real (atan2 (xy
.x
, xy
.y
));
4999 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5000 (SCM real
, SCM imaginary
),
5001 "Return a complex number constructed of the given @var{real} and\n"
5002 "@var{imaginary} parts.")
5003 #define FUNC_NAME s_scm_make_rectangular
5006 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5007 return scm_make_complex (xy
.x
, xy
.y
);
5013 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5015 "Return the complex number @var{x} * e^(i * @var{y}).")
5016 #define FUNC_NAME s_scm_make_polar
5020 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5022 sincos (xy
.y
, &s
, &c
);
5027 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
5032 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5033 /* "Return the real part of the number @var{z}."
5036 scm_real_part (SCM z
)
5040 else if (SCM_BIGP (z
))
5042 else if (SCM_REALP (z
))
5044 else if (SCM_COMPLEXP (z
))
5045 return scm_make_real (SCM_COMPLEX_REAL (z
));
5046 else if (SCM_FRACTIONP (z
))
5049 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5053 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5054 /* "Return the imaginary part of the number @var{z}."
5057 scm_imag_part (SCM z
)
5061 else if (SCM_BIGP (z
))
5063 else if (SCM_REALP (z
))
5065 else if (SCM_COMPLEXP (z
))
5066 return scm_make_real (SCM_COMPLEX_IMAG (z
));
5067 else if (SCM_FRACTIONP (z
))
5070 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5073 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5074 /* "Return the numerator of the number @var{z}."
5077 scm_numerator (SCM z
)
5081 else if (SCM_BIGP (z
))
5083 else if (SCM_FRACTIONP (z
))
5085 scm_i_fraction_reduce (z
);
5086 return SCM_FRACTION_NUMERATOR (z
);
5088 else if (SCM_REALP (z
))
5089 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5091 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5095 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5096 /* "Return the denominator of the number @var{z}."
5099 scm_denominator (SCM z
)
5102 return SCM_MAKINUM (1);
5103 else if (SCM_BIGP (z
))
5104 return SCM_MAKINUM (1);
5105 else if (SCM_FRACTIONP (z
))
5107 scm_i_fraction_reduce (z
);
5108 return SCM_FRACTION_DENOMINATOR (z
);
5110 else if (SCM_REALP (z
))
5111 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5113 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5116 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5117 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5118 * "@code{abs} for real arguments, but also allows complex numbers."
5121 scm_magnitude (SCM z
)
5125 long int zz
= SCM_INUM (z
);
5128 else if (SCM_POSFIXABLE (-zz
))
5129 return SCM_MAKINUM (-zz
);
5131 return scm_i_long2big (-zz
);
5133 else if (SCM_BIGP (z
))
5135 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5136 scm_remember_upto_here_1 (z
);
5138 return scm_i_clonebig (z
, 0);
5142 else if (SCM_REALP (z
))
5143 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5144 else if (SCM_COMPLEXP (z
))
5145 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5146 else if (SCM_FRACTIONP (z
))
5148 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5150 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5151 SCM_FRACTION_DENOMINATOR (z
));
5154 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5158 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5159 /* "Return the angle of the complex number @var{z}."
5164 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5165 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5166 But if atan2 follows the floating point rounding mode, then the value
5167 is not a constant. Maybe it'd be close enough though. */
5170 if (SCM_INUM (z
) >= 0)
5173 return scm_make_real (atan2 (0.0, -1.0));
5175 else if (SCM_BIGP (z
))
5177 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5178 scm_remember_upto_here_1 (z
);
5180 return scm_make_real (atan2 (0.0, -1.0));
5184 else if (SCM_REALP (z
))
5186 if (SCM_REAL_VALUE (z
) >= 0)
5189 return scm_make_real (atan2 (0.0, -1.0));
5191 else if (SCM_COMPLEXP (z
))
5192 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5193 else if (SCM_FRACTIONP (z
))
5195 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5197 else return scm_make_real (atan2 (0.0, -1.0));
5200 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5204 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5205 /* Convert the number @var{x} to its inexact representation.\n"
5208 scm_exact_to_inexact (SCM z
)
5211 return scm_make_real ((double) SCM_INUM (z
));
5212 else if (SCM_BIGP (z
))
5213 return scm_make_real (scm_i_big2dbl (z
));
5214 else if (SCM_FRACTIONP (z
))
5215 return scm_make_real (scm_i_fraction2double (z
));
5216 else if (SCM_INEXACTP (z
))
5219 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5223 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5225 "Return an exact number that is numerically closest to @var{z}.")
5226 #define FUNC_NAME s_scm_inexact_to_exact
5230 else if (SCM_BIGP (z
))
5232 else if (SCM_REALP (z
))
5234 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5235 SCM_OUT_OF_RANGE (1, z
);
5242 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5243 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5244 scm_i_mpz2num (mpq_denref (frac
)));
5246 /* When scm_make_ratio throws, we leak the memory allocated
5253 else if (SCM_FRACTIONP (z
))
5256 SCM_WRONG_TYPE_ARG (1, z
);
5260 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5262 "Return an exact number that is within @var{err} of @var{x}.")
5263 #define FUNC_NAME s_scm_rationalize
5267 else if (SCM_BIGP (x
))
5269 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5271 /* Use continued fractions to find closest ratio. All
5272 arithmetic is done with exact numbers.
5275 SCM ex
= scm_inexact_to_exact (x
);
5276 SCM int_part
= scm_floor (ex
);
5277 SCM tt
= SCM_MAKINUM (1);
5278 SCM a1
= SCM_MAKINUM (0), a2
= SCM_MAKINUM (1), a
= SCM_MAKINUM (0);
5279 SCM b1
= SCM_MAKINUM (1), b2
= SCM_MAKINUM (0), b
= SCM_MAKINUM (0);
5283 if (!SCM_FALSEP (scm_num_eq_p (ex
, int_part
)))
5286 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5287 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5289 /* We stop after a million iterations just to be absolutely sure
5290 that we don't go into an infinite loop. The process normally
5291 converges after less than a dozen iterations.
5294 err
= scm_abs (err
);
5295 while (++i
< 1000000)
5297 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5298 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5299 if (SCM_FALSEP (scm_zero_p (b
)) && /* b != 0 */
5301 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5302 err
))) /* abs(x-a/b) <= err */
5304 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5305 if (SCM_FALSEP (scm_exact_p (x
))
5306 || SCM_FALSEP (scm_exact_p (err
)))
5307 return scm_exact_to_inexact (res
);
5311 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5313 tt
= scm_floor (rx
); /* tt = floor (rx) */
5319 scm_num_overflow (s_scm_rationalize
);
5322 SCM_WRONG_TYPE_ARG (1, x
);
5326 /* if you need to change this, change test-num2integral.c as well */
5327 #if SCM_SIZEOF_LONG_LONG != 0
5329 # define ULLONG_MAX ((unsigned long long) (-1))
5330 # define LLONG_MAX ((long long) (ULLONG_MAX >> 1))
5331 # define LLONG_MIN (~LLONG_MAX)
5335 /* Parameters for creating integer conversion routines.
5337 Define the following preprocessor macros before including
5338 "libguile/num2integral.i.c":
5340 NUM2INTEGRAL - the name of the function for converting from a
5341 Scheme object to the integral type. This function will be
5342 defined when including "num2integral.i.c".
5344 INTEGRAL2NUM - the name of the function for converting from the
5345 integral type to a Scheme object. This function will be defined.
5347 INTEGRAL2BIG - the name of an internal function that createas a
5348 bignum from the integral type. This function will be defined.
5349 The name should start with "scm_i_".
5351 ITYPE - the name of the integral type.
5353 UNSIGNED - Define this to 1 when ITYPE is an unsigned type. Define
5356 UNSIGNED_ITYPE - the name of the the unsigned variant of the
5357 integral type. If you don't define this, it defaults to
5358 "unsigned ITYPE" for signed types and simply "ITYPE" for unsigned
5361 SIZEOF_ITYPE - an expression giving the size of the integral type
5362 in bytes. This expression must be computable by the
5363 preprocessor. (SIZEOF_FOO values are calculated by configure.in
5368 #define NUM2INTEGRAL scm_num2short
5369 #define INTEGRAL2NUM scm_short2num
5370 #define INTEGRAL2BIG scm_i_short2big
5373 #define SIZEOF_ITYPE SIZEOF_SHORT
5374 #include "libguile/num2integral.i.c"
5376 #define NUM2INTEGRAL scm_num2ushort
5377 #define INTEGRAL2NUM scm_ushort2num
5378 #define INTEGRAL2BIG scm_i_ushort2big
5380 #define ITYPE unsigned short
5381 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_SHORT
5382 #include "libguile/num2integral.i.c"
5384 #define NUM2INTEGRAL scm_num2int
5385 #define INTEGRAL2NUM scm_int2num
5386 #define INTEGRAL2BIG scm_i_int2big
5389 #define SIZEOF_ITYPE SIZEOF_INT
5390 #include "libguile/num2integral.i.c"
5392 #define NUM2INTEGRAL scm_num2uint
5393 #define INTEGRAL2NUM scm_uint2num
5394 #define INTEGRAL2BIG scm_i_uint2big
5396 #define ITYPE unsigned int
5397 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_INT
5398 #include "libguile/num2integral.i.c"
5400 #define NUM2INTEGRAL scm_num2long
5401 #define INTEGRAL2NUM scm_long2num
5402 #define INTEGRAL2BIG scm_i_long2big
5405 #define SIZEOF_ITYPE SIZEOF_LONG
5406 #include "libguile/num2integral.i.c"
5408 #define NUM2INTEGRAL scm_num2ulong
5409 #define INTEGRAL2NUM scm_ulong2num
5410 #define INTEGRAL2BIG scm_i_ulong2big
5412 #define ITYPE unsigned long
5413 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG
5414 #include "libguile/num2integral.i.c"
5416 #define NUM2INTEGRAL scm_num2ptrdiff
5417 #define INTEGRAL2NUM scm_ptrdiff2num
5418 #define INTEGRAL2BIG scm_i_ptrdiff2big
5420 #define ITYPE scm_t_ptrdiff
5421 #define UNSIGNED_ITYPE size_t
5422 #define SIZEOF_ITYPE SCM_SIZEOF_SCM_T_PTRDIFF
5423 #include "libguile/num2integral.i.c"
5425 #define NUM2INTEGRAL scm_num2size
5426 #define INTEGRAL2NUM scm_size2num
5427 #define INTEGRAL2BIG scm_i_size2big
5429 #define ITYPE size_t
5430 #define SIZEOF_ITYPE SIZEOF_SIZE_T
5431 #include "libguile/num2integral.i.c"
5433 #if SCM_SIZEOF_LONG_LONG != 0
5435 #ifndef ULONG_LONG_MAX
5436 #define ULONG_LONG_MAX (~0ULL)
5439 #define NUM2INTEGRAL scm_num2long_long
5440 #define INTEGRAL2NUM scm_long_long2num
5441 #define INTEGRAL2BIG scm_i_long_long2big
5443 #define ITYPE long long
5444 #define SIZEOF_ITYPE SIZEOF_LONG_LONG
5445 #include "libguile/num2integral.i.c"
5447 #define NUM2INTEGRAL scm_num2ulong_long
5448 #define INTEGRAL2NUM scm_ulong_long2num
5449 #define INTEGRAL2BIG scm_i_ulong_long2big
5451 #define ITYPE unsigned long long
5452 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG_LONG
5453 #include "libguile/num2integral.i.c"
5455 #endif /* SCM_SIZEOF_LONG_LONG != 0 */
5457 #define NUM2FLOAT scm_num2float
5458 #define FLOAT2NUM scm_float2num
5460 #include "libguile/num2float.i.c"
5462 #define NUM2FLOAT scm_num2double
5463 #define FLOAT2NUM scm_double2num
5464 #define FTYPE double
5465 #include "libguile/num2float.i.c"
5470 #define SIZE_MAX ((size_t) (-1))
5473 #define PTRDIFF_MIN \
5474 ((scm_t_ptrdiff) ((scm_t_ptrdiff) 1 \
5475 << ((sizeof (scm_t_ptrdiff) * SCM_CHAR_BIT) - 1)))
5478 #define PTRDIFF_MAX (~ PTRDIFF_MIN)
5481 #define CHECK(type, v) \
5484 if ((v) != scm_num2##type (scm_##type##2num (v), 1, "check_sanity")) \
5504 CHECK (ptrdiff
, -1);
5506 CHECK (short, SHRT_MAX
);
5507 CHECK (short, SHRT_MIN
);
5508 CHECK (ushort
, USHRT_MAX
);
5509 CHECK (int, INT_MAX
);
5510 CHECK (int, INT_MIN
);
5511 CHECK (uint
, UINT_MAX
);
5512 CHECK (long, LONG_MAX
);
5513 CHECK (long, LONG_MIN
);
5514 CHECK (ulong
, ULONG_MAX
);
5515 CHECK (size
, SIZE_MAX
);
5516 CHECK (ptrdiff
, PTRDIFF_MAX
);
5517 CHECK (ptrdiff
, PTRDIFF_MIN
);
5519 #if SCM_SIZEOF_LONG_LONG != 0
5520 CHECK (long_long
, 0LL);
5521 CHECK (ulong_long
, 0ULL);
5522 CHECK (long_long
, -1LL);
5523 CHECK (long_long
, LLONG_MAX
);
5524 CHECK (long_long
, LLONG_MIN
);
5525 CHECK (ulong_long
, ULLONG_MAX
);
5532 scm_internal_catch (SCM_BOOL_T, check_body, &data, check_handler, &data); \
5533 if (!SCM_FALSEP (data)) abort();
5536 check_body (void *data
)
5538 SCM num
= *(SCM
*) data
;
5539 scm_num2ulong (num
, 1, NULL
);
5541 return SCM_UNSPECIFIED
;
5545 check_handler (void *data
, SCM tag
, SCM throw_args
)
5547 SCM
*num
= (SCM
*) data
;
5550 return SCM_UNSPECIFIED
;
5553 SCM_DEFINE (scm_sys_check_number_conversions
, "%check-number-conversions", 0, 0, 0,
5555 "Number conversion sanity checking.")
5556 #define FUNC_NAME s_scm_sys_check_number_conversions
5558 SCM data
= SCM_MAKINUM (-1);
5560 data
= scm_int2num (INT_MIN
);
5562 data
= scm_ulong2num (ULONG_MAX
);
5563 data
= scm_difference (SCM_INUM0
, data
);
5565 data
= scm_ulong2num (ULONG_MAX
);
5566 data
= scm_sum (SCM_MAKINUM (1), data
); data
= scm_difference (SCM_INUM0
, data
);
5568 data
= scm_int2num (-10000); data
= scm_product (data
, data
); data
= scm_product (data
, data
);
5571 return SCM_UNSPECIFIED
;
5580 abs_most_negative_fixnum
= scm_i_long2big (- SCM_MOST_NEGATIVE_FIXNUM
);
5581 scm_permanent_object (abs_most_negative_fixnum
);
5583 mpz_init_set_si (z_negative_one
, -1);
5585 /* It may be possible to tune the performance of some algorithms by using
5586 * the following constants to avoid the creation of bignums. Please, before
5587 * using these values, remember the two rules of program optimization:
5588 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5589 scm_c_define ("most-positive-fixnum",
5590 SCM_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5591 scm_c_define ("most-negative-fixnum",
5592 SCM_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5594 scm_add_feature ("complex");
5595 scm_add_feature ("inexact");
5596 scm_flo0
= scm_make_real (0.0);
5598 scm_dblprec
= (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5600 { /* determine floating point precision */
5602 double fsum
= 1.0 + f
;
5605 if (++scm_dblprec
> 20)
5613 scm_dblprec
= scm_dblprec
- 1;
5615 #endif /* DBL_DIG */
5621 exactly_one_half
= scm_permanent_object (scm_divide (SCM_MAKINUM (1),
5623 #include "libguile/numbers.x"