1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010, 2011 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 License
9 * as published by the Free Software Foundation; either version 3 of
10 * the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful, but
13 * 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., 51 Franklin Street, Fifth Floor, Boston, MA
24 /* General assumptions:
25 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
26 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
27 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
28 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
29 * All objects satisfying SCM_FRACTIONP are never an integer.
34 - see if special casing bignums and reals in integer-exponent when
35 possible (to use mpz_pow and mpf_pow_ui) is faster.
37 - look in to better short-circuiting of common cases in
38 integer-expt and elsewhere.
40 - see if direct mpz operations can help in ash and elsewhere.
57 #include "libguile/_scm.h"
58 #include "libguile/feature.h"
59 #include "libguile/ports.h"
60 #include "libguile/root.h"
61 #include "libguile/smob.h"
62 #include "libguile/strings.h"
63 #include "libguile/bdw-gc.h"
65 #include "libguile/validate.h"
66 #include "libguile/numbers.h"
67 #include "libguile/deprecation.h"
69 #include "libguile/eq.h"
71 /* values per glibc, if not already defined */
73 #define M_LOG10E 0.43429448190325182765
76 #define M_PI 3.14159265358979323846
79 typedef scm_t_signed_bits scm_t_inum
;
80 #define scm_from_inum(x) (scm_from_signed_integer (x))
82 /* Tests to see if a C double is neither infinite nor a NaN.
83 TODO: if it's available, use C99's isfinite(x) instead */
84 #define DOUBLE_IS_FINITE(x) (!isinf(x) && !isnan(x))
89 Wonder if this might be faster for some of our code? A switch on
90 the numtag would jump directly to the right case, and the
91 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
93 #define SCM_I_NUMTAG_NOTNUM 0
94 #define SCM_I_NUMTAG_INUM 1
95 #define SCM_I_NUMTAG_BIG scm_tc16_big
96 #define SCM_I_NUMTAG_REAL scm_tc16_real
97 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
98 #define SCM_I_NUMTAG(x) \
99 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
100 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
101 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
102 : SCM_I_NUMTAG_NOTNUM)))
104 /* the macro above will not work as is with fractions */
109 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
111 /* FLOBUFLEN is the maximum number of characters neccessary for the
112 * printed or scm_string representation of an inexact number.
114 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
117 #if !defined (HAVE_ASINH)
118 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
120 #if !defined (HAVE_ACOSH)
121 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
123 #if !defined (HAVE_ATANH)
124 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
127 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
128 an explicit check. In some future gmp (don't know what version number),
129 mpz_cmp_d is supposed to do this itself. */
131 #define xmpz_cmp_d(z, d) \
132 (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
134 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
138 #if defined (GUILE_I)
139 #if HAVE_COMPLEX_DOUBLE
141 /* For an SCM object Z which is a complex number (ie. satisfies
142 SCM_COMPLEXP), return its value as a C level "complex double". */
143 #define SCM_COMPLEX_VALUE(z) \
144 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
146 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
148 /* Convert a C "complex double" to an SCM value. */
150 scm_from_complex_double (complex double z
)
152 return scm_c_make_rectangular (creal (z
), cimag (z
));
155 #endif /* HAVE_COMPLEX_DOUBLE */
160 static mpz_t z_negative_one
;
163 /* Clear the `mpz_t' embedded in bignum PTR. */
165 finalize_bignum (GC_PTR ptr
, GC_PTR data
)
169 bignum
= PTR2SCM (ptr
);
170 mpz_clear (SCM_I_BIG_MPZ (bignum
));
173 /* Return a new uninitialized bignum. */
178 GC_finalization_proc prev_finalizer
;
179 GC_PTR prev_finalizer_data
;
181 /* Allocate one word for the type tag and enough room for an `mpz_t'. */
182 p
= scm_gc_malloc_pointerless (sizeof (scm_t_bits
) + sizeof (mpz_t
),
186 GC_REGISTER_FINALIZER_NO_ORDER (p
, finalize_bignum
, NULL
,
188 &prev_finalizer_data
);
197 /* Return a newly created bignum. */
198 SCM z
= make_bignum ();
199 mpz_init (SCM_I_BIG_MPZ (z
));
204 scm_i_inum2big (scm_t_inum x
)
206 /* Return a newly created bignum initialized to X. */
207 SCM z
= make_bignum ();
208 #if SIZEOF_VOID_P == SIZEOF_LONG
209 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
211 /* Note that in this case, you'll also have to check all mpz_*_ui and
212 mpz_*_si invocations in Guile. */
213 #error creation of mpz not implemented for this inum size
219 scm_i_long2big (long x
)
221 /* Return a newly created bignum initialized to X. */
222 SCM z
= make_bignum ();
223 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
228 scm_i_ulong2big (unsigned long x
)
230 /* Return a newly created bignum initialized to X. */
231 SCM z
= make_bignum ();
232 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
237 scm_i_clonebig (SCM src_big
, int same_sign_p
)
239 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
240 SCM z
= make_bignum ();
241 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
243 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
248 scm_i_bigcmp (SCM x
, SCM y
)
250 /* Return neg if x < y, pos if x > y, and 0 if x == y */
251 /* presume we already know x and y are bignums */
252 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
253 scm_remember_upto_here_2 (x
, y
);
258 scm_i_dbl2big (double d
)
260 /* results are only defined if d is an integer */
261 SCM z
= make_bignum ();
262 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
266 /* Convert a integer in double representation to a SCM number. */
269 scm_i_dbl2num (double u
)
271 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
272 powers of 2, so there's no rounding when making "double" values
273 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
274 get rounded on a 64-bit machine, hence the "+1".
276 The use of floor() to force to an integer value ensures we get a
277 "numerically closest" value without depending on how a
278 double->long cast or how mpz_set_d will round. For reference,
279 double->long probably follows the hardware rounding mode,
280 mpz_set_d truncates towards zero. */
282 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
283 representable as a double? */
285 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
286 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
287 return SCM_I_MAKINUM ((scm_t_inum
) u
);
289 return scm_i_dbl2big (u
);
292 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
293 with R5RS exact->inexact.
295 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
296 (ie. truncate towards zero), then adjust to get the closest double by
297 examining the next lower bit and adding 1 (to the absolute value) if
300 Bignums exactly half way between representable doubles are rounded to the
301 next higher absolute value (ie. away from zero). This seems like an
302 adequate interpretation of R5RS "numerically closest", and it's easier
303 and faster than a full "nearest-even" style.
305 The bit test must be done on the absolute value of the mpz_t, which means
306 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
307 negatives as twos complement.
309 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
310 following the hardware rounding mode, but applied to the absolute value
311 of the mpz_t operand. This is not what we want so we put the high
312 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
313 mpz_get_d is supposed to always truncate towards zero.
315 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
316 is a slowdown. It'd be faster to pick out the relevant high bits with
317 mpz_getlimbn if we could be bothered coding that, and if the new
318 truncating gmp doesn't come out. */
321 scm_i_big2dbl (SCM b
)
326 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
330 /* Current GMP, eg. 4.1.3, force truncation towards zero */
332 if (bits
> DBL_MANT_DIG
)
334 size_t shift
= bits
- DBL_MANT_DIG
;
335 mpz_init2 (tmp
, DBL_MANT_DIG
);
336 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
337 result
= ldexp (mpz_get_d (tmp
), shift
);
342 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
347 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
350 if (bits
> DBL_MANT_DIG
)
352 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
353 /* test bit number "pos" in absolute value */
354 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
355 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
357 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
361 scm_remember_upto_here_1 (b
);
366 scm_i_normbig (SCM b
)
368 /* convert a big back to a fixnum if it'll fit */
369 /* presume b is a bignum */
370 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
372 scm_t_inum val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
373 if (SCM_FIXABLE (val
))
374 b
= SCM_I_MAKINUM (val
);
379 static SCM_C_INLINE_KEYWORD SCM
380 scm_i_mpz2num (mpz_t b
)
382 /* convert a mpz number to a SCM number. */
383 if (mpz_fits_slong_p (b
))
385 scm_t_inum val
= mpz_get_si (b
);
386 if (SCM_FIXABLE (val
))
387 return SCM_I_MAKINUM (val
);
391 SCM z
= make_bignum ();
392 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
397 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
398 static SCM
scm_divide2real (SCM x
, SCM y
);
401 scm_i_make_ratio (SCM numerator
, SCM denominator
)
402 #define FUNC_NAME "make-ratio"
404 /* First make sure the arguments are proper.
406 if (SCM_I_INUMP (denominator
))
408 if (scm_is_eq (denominator
, SCM_INUM0
))
409 scm_num_overflow ("make-ratio");
410 if (scm_is_eq (denominator
, SCM_INUM1
))
415 if (!(SCM_BIGP(denominator
)))
416 SCM_WRONG_TYPE_ARG (2, denominator
);
418 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
419 SCM_WRONG_TYPE_ARG (1, numerator
);
421 /* Then flip signs so that the denominator is positive.
423 if (scm_is_true (scm_negative_p (denominator
)))
425 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
426 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
429 /* Now consider for each of the four fixnum/bignum combinations
430 whether the rational number is really an integer.
432 if (SCM_I_INUMP (numerator
))
434 scm_t_inum x
= SCM_I_INUM (numerator
);
435 if (scm_is_eq (numerator
, SCM_INUM0
))
437 if (SCM_I_INUMP (denominator
))
440 y
= SCM_I_INUM (denominator
);
444 return SCM_I_MAKINUM (x
/ y
);
448 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
449 of that value for the denominator, as a bignum. Apart from
450 that case, abs(bignum) > abs(inum) so inum/bignum is not an
452 if (x
== SCM_MOST_NEGATIVE_FIXNUM
453 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
454 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
455 return SCM_I_MAKINUM(-1);
458 else if (SCM_BIGP (numerator
))
460 if (SCM_I_INUMP (denominator
))
462 scm_t_inum yy
= SCM_I_INUM (denominator
);
463 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
464 return scm_divide (numerator
, denominator
);
468 if (scm_is_eq (numerator
, denominator
))
470 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
471 SCM_I_BIG_MPZ (denominator
)))
472 return scm_divide(numerator
, denominator
);
476 /* No, it's a proper fraction.
479 SCM divisor
= scm_gcd (numerator
, denominator
);
480 if (!(scm_is_eq (divisor
, SCM_INUM1
)))
482 numerator
= scm_divide (numerator
, divisor
);
483 denominator
= scm_divide (denominator
, divisor
);
486 return scm_double_cell (scm_tc16_fraction
,
487 SCM_UNPACK (numerator
),
488 SCM_UNPACK (denominator
), 0);
494 scm_i_fraction2double (SCM z
)
496 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
497 SCM_FRACTION_DENOMINATOR (z
)));
500 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
502 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
504 #define FUNC_NAME s_scm_exact_p
510 if (SCM_FRACTIONP (x
))
514 SCM_WRONG_TYPE_ARG (1, x
);
519 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
521 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
523 #define FUNC_NAME s_scm_odd_p
527 scm_t_inum val
= SCM_I_INUM (n
);
528 return scm_from_bool ((val
& 1L) != 0);
530 else if (SCM_BIGP (n
))
532 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
533 scm_remember_upto_here_1 (n
);
534 return scm_from_bool (odd_p
);
536 else if (scm_is_true (scm_inf_p (n
)))
537 SCM_WRONG_TYPE_ARG (1, n
);
538 else if (SCM_REALP (n
))
540 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
546 SCM_WRONG_TYPE_ARG (1, n
);
549 SCM_WRONG_TYPE_ARG (1, n
);
554 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
556 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
558 #define FUNC_NAME s_scm_even_p
562 scm_t_inum val
= SCM_I_INUM (n
);
563 return scm_from_bool ((val
& 1L) == 0);
565 else if (SCM_BIGP (n
))
567 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
568 scm_remember_upto_here_1 (n
);
569 return scm_from_bool (even_p
);
571 else if (scm_is_true (scm_inf_p (n
)))
572 SCM_WRONG_TYPE_ARG (1, n
);
573 else if (SCM_REALP (n
))
575 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
581 SCM_WRONG_TYPE_ARG (1, n
);
584 SCM_WRONG_TYPE_ARG (1, n
);
588 SCM_DEFINE (scm_finite_p
, "finite?", 1, 0, 0,
590 "Return @code{#t} if @var{x} is neither infinite\n"
591 "nor a NaN, @code{#f} otherwise.")
592 #define FUNC_NAME s_scm_finite_p
595 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
596 else if (SCM_COMPLEXP (x
))
597 return scm_from_bool (DOUBLE_IS_FINITE (SCM_COMPLEX_REAL (x
))
598 && DOUBLE_IS_FINITE (SCM_COMPLEX_IMAG (x
)));
599 else if (SCM_NUMBERP (x
))
602 SCM_WRONG_TYPE_ARG (1, x
);
606 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
608 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
609 "or @samp{-inf.0}, @code{#f} otherwise.")
610 #define FUNC_NAME s_scm_inf_p
613 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
614 else if (SCM_COMPLEXP (x
))
615 return scm_from_bool (isinf (SCM_COMPLEX_REAL (x
))
616 || isinf (SCM_COMPLEX_IMAG (x
)));
622 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
624 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
626 #define FUNC_NAME s_scm_nan_p
629 return scm_from_bool (isnan (SCM_REAL_VALUE (n
)));
630 else if (SCM_COMPLEXP (n
))
631 return scm_from_bool (isnan (SCM_COMPLEX_REAL (n
))
632 || isnan (SCM_COMPLEX_IMAG (n
)));
638 /* Guile's idea of infinity. */
639 static double guile_Inf
;
641 /* Guile's idea of not a number. */
642 static double guile_NaN
;
645 guile_ieee_init (void)
647 /* Some version of gcc on some old version of Linux used to crash when
648 trying to make Inf and NaN. */
651 /* C99 INFINITY, when available.
652 FIXME: The standard allows for INFINITY to be something that overflows
653 at compile time. We ought to have a configure test to check for that
654 before trying to use it. (But in practice we believe this is not a
655 problem on any system guile is likely to target.) */
656 guile_Inf
= INFINITY
;
657 #elif defined HAVE_DINFINITY
659 extern unsigned int DINFINITY
[2];
660 guile_Inf
= (*((double *) (DINFINITY
)));
667 if (guile_Inf
== tmp
)
674 /* C99 NAN, when available */
676 #elif defined HAVE_DQNAN
679 extern unsigned int DQNAN
[2];
680 guile_NaN
= (*((double *)(DQNAN
)));
683 guile_NaN
= guile_Inf
/ guile_Inf
;
687 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
690 #define FUNC_NAME s_scm_inf
692 static int initialized
= 0;
698 return scm_from_double (guile_Inf
);
702 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
705 #define FUNC_NAME s_scm_nan
707 static int initialized
= 0;
713 return scm_from_double (guile_NaN
);
718 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
720 "Return the absolute value of @var{x}.")
725 scm_t_inum xx
= SCM_I_INUM (x
);
728 else if (SCM_POSFIXABLE (-xx
))
729 return SCM_I_MAKINUM (-xx
);
731 return scm_i_inum2big (-xx
);
733 else if (SCM_BIGP (x
))
735 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
737 return scm_i_clonebig (x
, 0);
741 else if (SCM_REALP (x
))
743 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
744 double xx
= SCM_REAL_VALUE (x
);
746 return scm_from_double (-xx
);
750 else if (SCM_FRACTIONP (x
))
752 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
754 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
755 SCM_FRACTION_DENOMINATOR (x
));
758 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
763 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
764 /* "Return the quotient of the numbers @var{x} and @var{y}."
767 scm_quotient (SCM x
, SCM y
)
771 scm_t_inum xx
= SCM_I_INUM (x
);
774 scm_t_inum yy
= SCM_I_INUM (y
);
776 scm_num_overflow (s_quotient
);
779 scm_t_inum z
= xx
/ yy
;
781 return SCM_I_MAKINUM (z
);
783 return scm_i_inum2big (z
);
786 else if (SCM_BIGP (y
))
788 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
789 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
790 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
792 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
793 scm_remember_upto_here_1 (y
);
794 return SCM_I_MAKINUM (-1);
800 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
802 else if (SCM_BIGP (x
))
806 scm_t_inum yy
= SCM_I_INUM (y
);
808 scm_num_overflow (s_quotient
);
813 SCM result
= scm_i_mkbig ();
816 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
819 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
822 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
823 scm_remember_upto_here_1 (x
);
824 return scm_i_normbig (result
);
827 else if (SCM_BIGP (y
))
829 SCM result
= scm_i_mkbig ();
830 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
833 scm_remember_upto_here_2 (x
, y
);
834 return scm_i_normbig (result
);
837 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
840 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
843 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
844 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
846 * "(remainder 13 4) @result{} 1\n"
847 * "(remainder -13 4) @result{} -1\n"
851 scm_remainder (SCM x
, SCM y
)
857 scm_t_inum yy
= SCM_I_INUM (y
);
859 scm_num_overflow (s_remainder
);
862 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
863 return SCM_I_MAKINUM (z
);
866 else if (SCM_BIGP (y
))
868 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
869 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
870 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
872 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
873 scm_remember_upto_here_1 (y
);
880 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
882 else if (SCM_BIGP (x
))
886 scm_t_inum yy
= SCM_I_INUM (y
);
888 scm_num_overflow (s_remainder
);
891 SCM result
= scm_i_mkbig ();
894 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
895 scm_remember_upto_here_1 (x
);
896 return scm_i_normbig (result
);
899 else if (SCM_BIGP (y
))
901 SCM result
= scm_i_mkbig ();
902 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
905 scm_remember_upto_here_2 (x
, y
);
906 return scm_i_normbig (result
);
909 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
912 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
916 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
917 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
919 * "(modulo 13 4) @result{} 1\n"
920 * "(modulo -13 4) @result{} 3\n"
924 scm_modulo (SCM x
, SCM y
)
928 scm_t_inum xx
= SCM_I_INUM (x
);
931 scm_t_inum yy
= SCM_I_INUM (y
);
933 scm_num_overflow (s_modulo
);
936 /* C99 specifies that "%" is the remainder corresponding to a
937 quotient rounded towards zero, and that's also traditional
938 for machine division, so z here should be well defined. */
939 scm_t_inum z
= xx
% yy
;
956 return SCM_I_MAKINUM (result
);
959 else if (SCM_BIGP (y
))
961 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
968 SCM pos_y
= scm_i_clonebig (y
, 0);
969 /* do this after the last scm_op */
970 mpz_init_set_si (z_x
, xx
);
971 result
= pos_y
; /* re-use this bignum */
972 mpz_mod (SCM_I_BIG_MPZ (result
),
974 SCM_I_BIG_MPZ (pos_y
));
975 scm_remember_upto_here_1 (pos_y
);
979 result
= scm_i_mkbig ();
980 /* do this after the last scm_op */
981 mpz_init_set_si (z_x
, xx
);
982 mpz_mod (SCM_I_BIG_MPZ (result
),
985 scm_remember_upto_here_1 (y
);
988 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
989 mpz_add (SCM_I_BIG_MPZ (result
),
991 SCM_I_BIG_MPZ (result
));
992 scm_remember_upto_here_1 (y
);
993 /* and do this before the next one */
995 return scm_i_normbig (result
);
999 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1001 else if (SCM_BIGP (x
))
1003 if (SCM_I_INUMP (y
))
1005 scm_t_inum yy
= SCM_I_INUM (y
);
1007 scm_num_overflow (s_modulo
);
1010 SCM result
= scm_i_mkbig ();
1011 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1013 (yy
< 0) ? - yy
: yy
);
1014 scm_remember_upto_here_1 (x
);
1015 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1016 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1017 SCM_I_BIG_MPZ (result
),
1019 return scm_i_normbig (result
);
1022 else if (SCM_BIGP (y
))
1025 SCM result
= scm_i_mkbig ();
1026 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1027 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1028 mpz_mod (SCM_I_BIG_MPZ (result
),
1030 SCM_I_BIG_MPZ (pos_y
));
1032 scm_remember_upto_here_1 (x
);
1033 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1034 mpz_add (SCM_I_BIG_MPZ (result
),
1036 SCM_I_BIG_MPZ (result
));
1037 scm_remember_upto_here_2 (y
, pos_y
);
1038 return scm_i_normbig (result
);
1042 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1045 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1048 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1049 (SCM x
, SCM y
, SCM rest
),
1050 "Return the greatest common divisor of all parameter values.\n"
1051 "If called without arguments, 0 is returned.")
1052 #define FUNC_NAME s_scm_i_gcd
1054 while (!scm_is_null (rest
))
1055 { x
= scm_gcd (x
, y
);
1057 rest
= scm_cdr (rest
);
1059 return scm_gcd (x
, y
);
1063 #define s_gcd s_scm_i_gcd
1064 #define g_gcd g_scm_i_gcd
1067 scm_gcd (SCM x
, SCM y
)
1070 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1072 if (SCM_I_INUMP (x
))
1074 if (SCM_I_INUMP (y
))
1076 scm_t_inum xx
= SCM_I_INUM (x
);
1077 scm_t_inum yy
= SCM_I_INUM (y
);
1078 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
1079 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
1089 /* Determine a common factor 2^k */
1090 while (!(1 & (u
| v
)))
1096 /* Now, any factor 2^n can be eliminated */
1116 return (SCM_POSFIXABLE (result
)
1117 ? SCM_I_MAKINUM (result
)
1118 : scm_i_inum2big (result
));
1120 else if (SCM_BIGP (y
))
1126 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1128 else if (SCM_BIGP (x
))
1130 if (SCM_I_INUMP (y
))
1135 yy
= SCM_I_INUM (y
);
1140 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1141 scm_remember_upto_here_1 (x
);
1142 return (SCM_POSFIXABLE (result
)
1143 ? SCM_I_MAKINUM (result
)
1144 : scm_from_unsigned_integer (result
));
1146 else if (SCM_BIGP (y
))
1148 SCM result
= scm_i_mkbig ();
1149 mpz_gcd (SCM_I_BIG_MPZ (result
),
1152 scm_remember_upto_here_2 (x
, y
);
1153 return scm_i_normbig (result
);
1156 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1159 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1162 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1163 (SCM x
, SCM y
, SCM rest
),
1164 "Return the least common multiple of the arguments.\n"
1165 "If called without arguments, 1 is returned.")
1166 #define FUNC_NAME s_scm_i_lcm
1168 while (!scm_is_null (rest
))
1169 { x
= scm_lcm (x
, y
);
1171 rest
= scm_cdr (rest
);
1173 return scm_lcm (x
, y
);
1177 #define s_lcm s_scm_i_lcm
1178 #define g_lcm g_scm_i_lcm
1181 scm_lcm (SCM n1
, SCM n2
)
1183 if (SCM_UNBNDP (n2
))
1185 if (SCM_UNBNDP (n1
))
1186 return SCM_I_MAKINUM (1L);
1187 n2
= SCM_I_MAKINUM (1L);
1190 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1191 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1192 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1193 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1195 if (SCM_I_INUMP (n1
))
1197 if (SCM_I_INUMP (n2
))
1199 SCM d
= scm_gcd (n1
, n2
);
1200 if (scm_is_eq (d
, SCM_INUM0
))
1203 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1207 /* inum n1, big n2 */
1210 SCM result
= scm_i_mkbig ();
1211 scm_t_inum nn1
= SCM_I_INUM (n1
);
1212 if (nn1
== 0) return SCM_INUM0
;
1213 if (nn1
< 0) nn1
= - nn1
;
1214 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1215 scm_remember_upto_here_1 (n2
);
1223 if (SCM_I_INUMP (n2
))
1230 SCM result
= scm_i_mkbig ();
1231 mpz_lcm(SCM_I_BIG_MPZ (result
),
1233 SCM_I_BIG_MPZ (n2
));
1234 scm_remember_upto_here_2(n1
, n2
);
1235 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1241 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1246 + + + x (map digit:logand X Y)
1247 + - + x (map digit:logand X (lognot (+ -1 Y)))
1248 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1249 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1254 + + + (map digit:logior X Y)
1255 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1256 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1257 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1262 + + + (map digit:logxor X Y)
1263 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1264 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1265 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1270 + + (any digit:logand X Y)
1271 + - (any digit:logand X (lognot (+ -1 Y)))
1272 - + (any digit:logand (lognot (+ -1 X)) Y)
1277 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1278 (SCM x
, SCM y
, SCM rest
),
1279 "Return the bitwise AND of the integer arguments.\n\n"
1281 "(logand) @result{} -1\n"
1282 "(logand 7) @result{} 7\n"
1283 "(logand #b111 #b011 #b001) @result{} 1\n"
1285 #define FUNC_NAME s_scm_i_logand
1287 while (!scm_is_null (rest
))
1288 { x
= scm_logand (x
, y
);
1290 rest
= scm_cdr (rest
);
1292 return scm_logand (x
, y
);
1296 #define s_scm_logand s_scm_i_logand
1298 SCM
scm_logand (SCM n1
, SCM n2
)
1299 #define FUNC_NAME s_scm_logand
1303 if (SCM_UNBNDP (n2
))
1305 if (SCM_UNBNDP (n1
))
1306 return SCM_I_MAKINUM (-1);
1307 else if (!SCM_NUMBERP (n1
))
1308 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1309 else if (SCM_NUMBERP (n1
))
1312 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1315 if (SCM_I_INUMP (n1
))
1317 nn1
= SCM_I_INUM (n1
);
1318 if (SCM_I_INUMP (n2
))
1320 scm_t_inum nn2
= SCM_I_INUM (n2
);
1321 return SCM_I_MAKINUM (nn1
& nn2
);
1323 else if SCM_BIGP (n2
)
1329 SCM result_z
= scm_i_mkbig ();
1331 mpz_init_set_si (nn1_z
, nn1
);
1332 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1333 scm_remember_upto_here_1 (n2
);
1335 return scm_i_normbig (result_z
);
1339 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1341 else if (SCM_BIGP (n1
))
1343 if (SCM_I_INUMP (n2
))
1346 nn1
= SCM_I_INUM (n1
);
1349 else if (SCM_BIGP (n2
))
1351 SCM result_z
= scm_i_mkbig ();
1352 mpz_and (SCM_I_BIG_MPZ (result_z
),
1354 SCM_I_BIG_MPZ (n2
));
1355 scm_remember_upto_here_2 (n1
, n2
);
1356 return scm_i_normbig (result_z
);
1359 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1362 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1367 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1368 (SCM x
, SCM y
, SCM rest
),
1369 "Return the bitwise OR of the integer arguments.\n\n"
1371 "(logior) @result{} 0\n"
1372 "(logior 7) @result{} 7\n"
1373 "(logior #b000 #b001 #b011) @result{} 3\n"
1375 #define FUNC_NAME s_scm_i_logior
1377 while (!scm_is_null (rest
))
1378 { x
= scm_logior (x
, y
);
1380 rest
= scm_cdr (rest
);
1382 return scm_logior (x
, y
);
1386 #define s_scm_logior s_scm_i_logior
1388 SCM
scm_logior (SCM n1
, SCM n2
)
1389 #define FUNC_NAME s_scm_logior
1393 if (SCM_UNBNDP (n2
))
1395 if (SCM_UNBNDP (n1
))
1397 else if (SCM_NUMBERP (n1
))
1400 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1403 if (SCM_I_INUMP (n1
))
1405 nn1
= SCM_I_INUM (n1
);
1406 if (SCM_I_INUMP (n2
))
1408 long nn2
= SCM_I_INUM (n2
);
1409 return SCM_I_MAKINUM (nn1
| nn2
);
1411 else if (SCM_BIGP (n2
))
1417 SCM result_z
= scm_i_mkbig ();
1419 mpz_init_set_si (nn1_z
, nn1
);
1420 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1421 scm_remember_upto_here_1 (n2
);
1423 return scm_i_normbig (result_z
);
1427 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1429 else if (SCM_BIGP (n1
))
1431 if (SCM_I_INUMP (n2
))
1434 nn1
= SCM_I_INUM (n1
);
1437 else if (SCM_BIGP (n2
))
1439 SCM result_z
= scm_i_mkbig ();
1440 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1442 SCM_I_BIG_MPZ (n2
));
1443 scm_remember_upto_here_2 (n1
, n2
);
1444 return scm_i_normbig (result_z
);
1447 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1450 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1455 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1456 (SCM x
, SCM y
, SCM rest
),
1457 "Return the bitwise XOR of the integer arguments. A bit is\n"
1458 "set in the result if it is set in an odd number of arguments.\n"
1460 "(logxor) @result{} 0\n"
1461 "(logxor 7) @result{} 7\n"
1462 "(logxor #b000 #b001 #b011) @result{} 2\n"
1463 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1465 #define FUNC_NAME s_scm_i_logxor
1467 while (!scm_is_null (rest
))
1468 { x
= scm_logxor (x
, y
);
1470 rest
= scm_cdr (rest
);
1472 return scm_logxor (x
, y
);
1476 #define s_scm_logxor s_scm_i_logxor
1478 SCM
scm_logxor (SCM n1
, SCM n2
)
1479 #define FUNC_NAME s_scm_logxor
1483 if (SCM_UNBNDP (n2
))
1485 if (SCM_UNBNDP (n1
))
1487 else if (SCM_NUMBERP (n1
))
1490 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1493 if (SCM_I_INUMP (n1
))
1495 nn1
= SCM_I_INUM (n1
);
1496 if (SCM_I_INUMP (n2
))
1498 scm_t_inum nn2
= SCM_I_INUM (n2
);
1499 return SCM_I_MAKINUM (nn1
^ nn2
);
1501 else if (SCM_BIGP (n2
))
1505 SCM result_z
= scm_i_mkbig ();
1507 mpz_init_set_si (nn1_z
, nn1
);
1508 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1509 scm_remember_upto_here_1 (n2
);
1511 return scm_i_normbig (result_z
);
1515 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1517 else if (SCM_BIGP (n1
))
1519 if (SCM_I_INUMP (n2
))
1522 nn1
= SCM_I_INUM (n1
);
1525 else if (SCM_BIGP (n2
))
1527 SCM result_z
= scm_i_mkbig ();
1528 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1530 SCM_I_BIG_MPZ (n2
));
1531 scm_remember_upto_here_2 (n1
, n2
);
1532 return scm_i_normbig (result_z
);
1535 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1538 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1543 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1545 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1546 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1547 "without actually calculating the @code{logand}, just testing\n"
1551 "(logtest #b0100 #b1011) @result{} #f\n"
1552 "(logtest #b0100 #b0111) @result{} #t\n"
1554 #define FUNC_NAME s_scm_logtest
1558 if (SCM_I_INUMP (j
))
1560 nj
= SCM_I_INUM (j
);
1561 if (SCM_I_INUMP (k
))
1563 scm_t_inum nk
= SCM_I_INUM (k
);
1564 return scm_from_bool (nj
& nk
);
1566 else if (SCM_BIGP (k
))
1574 mpz_init_set_si (nj_z
, nj
);
1575 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1576 scm_remember_upto_here_1 (k
);
1577 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1583 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1585 else if (SCM_BIGP (j
))
1587 if (SCM_I_INUMP (k
))
1590 nj
= SCM_I_INUM (j
);
1593 else if (SCM_BIGP (k
))
1597 mpz_init (result_z
);
1601 scm_remember_upto_here_2 (j
, k
);
1602 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1603 mpz_clear (result_z
);
1607 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1610 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1615 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1617 "Test whether bit number @var{index} in @var{j} is set.\n"
1618 "@var{index} starts from 0 for the least significant bit.\n"
1621 "(logbit? 0 #b1101) @result{} #t\n"
1622 "(logbit? 1 #b1101) @result{} #f\n"
1623 "(logbit? 2 #b1101) @result{} #t\n"
1624 "(logbit? 3 #b1101) @result{} #t\n"
1625 "(logbit? 4 #b1101) @result{} #f\n"
1627 #define FUNC_NAME s_scm_logbit_p
1629 unsigned long int iindex
;
1630 iindex
= scm_to_ulong (index
);
1632 if (SCM_I_INUMP (j
))
1634 /* bits above what's in an inum follow the sign bit */
1635 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1636 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1638 else if (SCM_BIGP (j
))
1640 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1641 scm_remember_upto_here_1 (j
);
1642 return scm_from_bool (val
);
1645 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1650 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1652 "Return the integer which is the ones-complement of the integer\n"
1656 "(number->string (lognot #b10000000) 2)\n"
1657 " @result{} \"-10000001\"\n"
1658 "(number->string (lognot #b0) 2)\n"
1659 " @result{} \"-1\"\n"
1661 #define FUNC_NAME s_scm_lognot
1663 if (SCM_I_INUMP (n
)) {
1664 /* No overflow here, just need to toggle all the bits making up the inum.
1665 Enhancement: No need to strip the tag and add it back, could just xor
1666 a block of 1 bits, if that worked with the various debug versions of
1668 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1670 } else if (SCM_BIGP (n
)) {
1671 SCM result
= scm_i_mkbig ();
1672 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1673 scm_remember_upto_here_1 (n
);
1677 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1682 /* returns 0 if IN is not an integer. OUT must already be
1685 coerce_to_big (SCM in
, mpz_t out
)
1688 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1689 else if (SCM_I_INUMP (in
))
1690 mpz_set_si (out
, SCM_I_INUM (in
));
1697 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1698 (SCM n
, SCM k
, SCM m
),
1699 "Return @var{n} raised to the integer exponent\n"
1700 "@var{k}, modulo @var{m}.\n"
1703 "(modulo-expt 2 3 5)\n"
1706 #define FUNC_NAME s_scm_modulo_expt
1712 /* There are two classes of error we might encounter --
1713 1) Math errors, which we'll report by calling scm_num_overflow,
1715 2) wrong-type errors, which of course we'll report by calling
1717 We don't report those errors immediately, however; instead we do
1718 some cleanup first. These variables tell us which error (if
1719 any) we should report after cleaning up.
1721 int report_overflow
= 0;
1723 int position_of_wrong_type
= 0;
1724 SCM value_of_wrong_type
= SCM_INUM0
;
1726 SCM result
= SCM_UNDEFINED
;
1732 if (scm_is_eq (m
, SCM_INUM0
))
1734 report_overflow
= 1;
1738 if (!coerce_to_big (n
, n_tmp
))
1740 value_of_wrong_type
= n
;
1741 position_of_wrong_type
= 1;
1745 if (!coerce_to_big (k
, k_tmp
))
1747 value_of_wrong_type
= k
;
1748 position_of_wrong_type
= 2;
1752 if (!coerce_to_big (m
, m_tmp
))
1754 value_of_wrong_type
= m
;
1755 position_of_wrong_type
= 3;
1759 /* if the exponent K is negative, and we simply call mpz_powm, we
1760 will get a divide-by-zero exception when an inverse 1/n mod m
1761 doesn't exist (or is not unique). Since exceptions are hard to
1762 handle, we'll attempt the inversion "by hand" -- that way, we get
1763 a simple failure code, which is easy to handle. */
1765 if (-1 == mpz_sgn (k_tmp
))
1767 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1769 report_overflow
= 1;
1772 mpz_neg (k_tmp
, k_tmp
);
1775 result
= scm_i_mkbig ();
1776 mpz_powm (SCM_I_BIG_MPZ (result
),
1781 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1782 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1789 if (report_overflow
)
1790 scm_num_overflow (FUNC_NAME
);
1792 if (position_of_wrong_type
)
1793 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1794 value_of_wrong_type
);
1796 return scm_i_normbig (result
);
1800 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1802 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1803 "exact integer, @var{n} can be any number.\n"
1805 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1806 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1807 "includes @math{0^0} is 1.\n"
1810 "(integer-expt 2 5) @result{} 32\n"
1811 "(integer-expt -3 3) @result{} -27\n"
1812 "(integer-expt 5 -3) @result{} 1/125\n"
1813 "(integer-expt 0 0) @result{} 1\n"
1815 #define FUNC_NAME s_scm_integer_expt
1818 SCM z_i2
= SCM_BOOL_F
;
1820 SCM acc
= SCM_I_MAKINUM (1L);
1822 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1823 if (!SCM_I_INUMP (k
) && !SCM_BIGP (k
))
1824 SCM_WRONG_TYPE_ARG (2, k
);
1826 if (scm_is_true (scm_zero_p (n
)))
1828 if (scm_is_true (scm_zero_p (k
))) /* 0^0 == 1 per R5RS */
1829 return acc
; /* return exact 1, regardless of n */
1830 else if (scm_is_true (scm_positive_p (k
)))
1832 else /* return NaN for (0 ^ k) for negative k per R6RS */
1835 else if (scm_is_eq (n
, acc
))
1837 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1838 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1840 if (SCM_I_INUMP (k
))
1841 i2
= SCM_I_INUM (k
);
1842 else if (SCM_BIGP (k
))
1844 z_i2
= scm_i_clonebig (k
, 1);
1845 scm_remember_upto_here_1 (k
);
1849 SCM_WRONG_TYPE_ARG (2, k
);
1853 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1855 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1856 n
= scm_divide (n
, SCM_UNDEFINED
);
1860 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1864 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1866 return scm_product (acc
, n
);
1868 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1869 acc
= scm_product (acc
, n
);
1870 n
= scm_product (n
, n
);
1871 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1879 n
= scm_divide (n
, SCM_UNDEFINED
);
1886 return scm_product (acc
, n
);
1888 acc
= scm_product (acc
, n
);
1889 n
= scm_product (n
, n
);
1896 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1898 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1899 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1901 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1902 "@var{cnt} is negative it's a division, rounded towards negative\n"
1903 "infinity. (Note that this is not the same rounding as\n"
1904 "@code{quotient} does.)\n"
1906 "With @var{n} viewed as an infinite precision twos complement,\n"
1907 "@code{ash} means a left shift introducing zero bits, or a right\n"
1908 "shift dropping bits.\n"
1911 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1912 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1914 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1915 "(ash -23 -2) @result{} -6\n"
1917 #define FUNC_NAME s_scm_ash
1920 bits_to_shift
= scm_to_long (cnt
);
1922 if (SCM_I_INUMP (n
))
1924 scm_t_inum nn
= SCM_I_INUM (n
);
1926 if (bits_to_shift
> 0)
1928 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1929 overflow a non-zero fixnum. For smaller shifts we check the
1930 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1931 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1932 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1938 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1940 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1943 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1947 SCM result
= scm_i_inum2big (nn
);
1948 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1955 bits_to_shift
= -bits_to_shift
;
1956 if (bits_to_shift
>= SCM_LONG_BIT
)
1957 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
1959 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1963 else if (SCM_BIGP (n
))
1967 if (bits_to_shift
== 0)
1970 result
= scm_i_mkbig ();
1971 if (bits_to_shift
>= 0)
1973 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1979 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1980 we have to allocate a bignum even if the result is going to be a
1982 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1984 return scm_i_normbig (result
);
1990 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1996 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1997 (SCM n
, SCM start
, SCM end
),
1998 "Return the integer composed of the @var{start} (inclusive)\n"
1999 "through @var{end} (exclusive) bits of @var{n}. The\n"
2000 "@var{start}th bit becomes the 0-th bit in the result.\n"
2003 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
2004 " @result{} \"1010\"\n"
2005 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
2006 " @result{} \"10110\"\n"
2008 #define FUNC_NAME s_scm_bit_extract
2010 unsigned long int istart
, iend
, bits
;
2011 istart
= scm_to_ulong (start
);
2012 iend
= scm_to_ulong (end
);
2013 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
2015 /* how many bits to keep */
2016 bits
= iend
- istart
;
2018 if (SCM_I_INUMP (n
))
2020 scm_t_inum in
= SCM_I_INUM (n
);
2022 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
2023 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
2024 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
2026 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
2028 /* Since we emulate two's complement encoded numbers, this
2029 * special case requires us to produce a result that has
2030 * more bits than can be stored in a fixnum.
2032 SCM result
= scm_i_inum2big (in
);
2033 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2038 /* mask down to requisite bits */
2039 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2040 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2042 else if (SCM_BIGP (n
))
2047 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2051 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2052 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2053 such bits into a ulong. */
2054 result
= scm_i_mkbig ();
2055 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2056 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2057 result
= scm_i_normbig (result
);
2059 scm_remember_upto_here_1 (n
);
2063 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2068 static const char scm_logtab
[] = {
2069 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2072 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2074 "Return the number of bits in integer @var{n}. If integer is\n"
2075 "positive, the 1-bits in its binary representation are counted.\n"
2076 "If negative, the 0-bits in its two's-complement binary\n"
2077 "representation are counted. If 0, 0 is returned.\n"
2080 "(logcount #b10101010)\n"
2087 #define FUNC_NAME s_scm_logcount
2089 if (SCM_I_INUMP (n
))
2091 unsigned long c
= 0;
2092 scm_t_inum nn
= SCM_I_INUM (n
);
2097 c
+= scm_logtab
[15 & nn
];
2100 return SCM_I_MAKINUM (c
);
2102 else if (SCM_BIGP (n
))
2104 unsigned long count
;
2105 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2106 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2108 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2109 scm_remember_upto_here_1 (n
);
2110 return SCM_I_MAKINUM (count
);
2113 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2118 static const char scm_ilentab
[] = {
2119 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2123 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2125 "Return the number of bits necessary to represent @var{n}.\n"
2128 "(integer-length #b10101010)\n"
2130 "(integer-length 0)\n"
2132 "(integer-length #b1111)\n"
2135 #define FUNC_NAME s_scm_integer_length
2137 if (SCM_I_INUMP (n
))
2139 unsigned long c
= 0;
2141 scm_t_inum nn
= SCM_I_INUM (n
);
2147 l
= scm_ilentab
[15 & nn
];
2150 return SCM_I_MAKINUM (c
- 4 + l
);
2152 else if (SCM_BIGP (n
))
2154 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2155 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2156 1 too big, so check for that and adjust. */
2157 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2158 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2159 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2160 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2162 scm_remember_upto_here_1 (n
);
2163 return SCM_I_MAKINUM (size
);
2166 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2170 /*** NUMBERS -> STRINGS ***/
2171 #define SCM_MAX_DBL_PREC 60
2172 #define SCM_MAX_DBL_RADIX 36
2174 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2175 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2176 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2179 void init_dblprec(int *prec
, int radix
) {
2180 /* determine floating point precision by adding successively
2181 smaller increments to 1.0 until it is considered == 1.0 */
2182 double f
= ((double)1.0)/radix
;
2183 double fsum
= 1.0 + f
;
2188 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2200 void init_fx_radix(double *fx_list
, int radix
)
2202 /* initialize a per-radix list of tolerances. When added
2203 to a number < 1.0, we can determine if we should raund
2204 up and quit converting a number to a string. */
2208 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2209 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2212 /* use this array as a way to generate a single digit */
2213 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
2216 idbl2str (double f
, char *a
, int radix
)
2218 int efmt
, dpt
, d
, i
, wp
;
2220 #ifdef DBL_MIN_10_EXP
2223 #endif /* DBL_MIN_10_EXP */
2228 radix
> SCM_MAX_DBL_RADIX
)
2230 /* revert to existing behavior */
2234 wp
= scm_dblprec
[radix
-2];
2235 fx
= fx_per_radix
[radix
-2];
2239 #ifdef HAVE_COPYSIGN
2240 double sgn
= copysign (1.0, f
);
2245 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2251 strcpy (a
, "-inf.0");
2253 strcpy (a
, "+inf.0");
2258 strcpy (a
, "+nan.0");
2268 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2269 make-uniform-vector, from causing infinite loops. */
2270 /* just do the checking...if it passes, we do the conversion for our
2271 radix again below */
2278 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2286 while (f_cpy
> 10.0)
2289 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2310 if (f
+ fx
[wp
] >= radix
)
2317 /* adding 9999 makes this equivalent to abs(x) % 3 */
2318 dpt
= (exp
+ 9999) % 3;
2322 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2344 a
[ch
++] = number_chars
[d
];
2347 if (f
+ fx
[wp
] >= 1.0)
2349 a
[ch
- 1] = number_chars
[d
+1];
2361 if ((dpt
> 4) && (exp
> 6))
2363 d
= (a
[0] == '-' ? 2 : 1);
2364 for (i
= ch
++; i
> d
; i
--)
2377 if (a
[ch
- 1] == '.')
2378 a
[ch
++] = '0'; /* trailing zero */
2387 for (i
= radix
; i
<= exp
; i
*= radix
);
2388 for (i
/= radix
; i
; i
/= radix
)
2390 a
[ch
++] = number_chars
[exp
/ i
];
2399 icmplx2str (double real
, double imag
, char *str
, int radix
)
2403 i
= idbl2str (real
, str
, radix
);
2406 /* Don't output a '+' for negative numbers or for Inf and
2407 NaN. They will provide their own sign. */
2408 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
2410 i
+= idbl2str (imag
, &str
[i
], radix
);
2417 iflo2str (SCM flt
, char *str
, int radix
)
2420 if (SCM_REALP (flt
))
2421 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2423 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2428 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2429 characters in the result.
2431 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2433 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2438 return scm_iuint2str (-num
, rad
, p
) + 1;
2441 return scm_iuint2str (num
, rad
, p
);
2444 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2445 characters in the result.
2447 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2449 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2453 scm_t_uintmax n
= num
;
2455 if (rad
< 2 || rad
> 36)
2456 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
2458 for (n
/= rad
; n
> 0; n
/= rad
)
2468 p
[i
] = number_chars
[d
];
2473 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2475 "Return a string holding the external representation of the\n"
2476 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2477 "inexact, a radix of 10 will be used.")
2478 #define FUNC_NAME s_scm_number_to_string
2482 if (SCM_UNBNDP (radix
))
2485 base
= scm_to_signed_integer (radix
, 2, 36);
2487 if (SCM_I_INUMP (n
))
2489 char num_buf
[SCM_INTBUFLEN
];
2490 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2491 return scm_from_locale_stringn (num_buf
, length
);
2493 else if (SCM_BIGP (n
))
2495 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2496 scm_remember_upto_here_1 (n
);
2497 return scm_take_locale_string (str
);
2499 else if (SCM_FRACTIONP (n
))
2501 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2502 scm_from_locale_string ("/"),
2503 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2505 else if (SCM_INEXACTP (n
))
2507 char num_buf
[FLOBUFLEN
];
2508 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2511 SCM_WRONG_TYPE_ARG (1, n
);
2516 /* These print routines used to be stubbed here so that scm_repl.c
2517 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2520 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2522 char num_buf
[FLOBUFLEN
];
2523 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2528 scm_i_print_double (double val
, SCM port
)
2530 char num_buf
[FLOBUFLEN
];
2531 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2535 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2538 char num_buf
[FLOBUFLEN
];
2539 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2544 scm_i_print_complex (double real
, double imag
, SCM port
)
2546 char num_buf
[FLOBUFLEN
];
2547 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2551 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2554 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2555 scm_display (str
, port
);
2556 scm_remember_upto_here_1 (str
);
2561 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2563 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2564 scm_remember_upto_here_1 (exp
);
2565 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2569 /*** END nums->strs ***/
2572 /*** STRINGS -> NUMBERS ***/
2574 /* The following functions implement the conversion from strings to numbers.
2575 * The implementation somehow follows the grammar for numbers as it is given
2576 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2577 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2578 * points should be noted about the implementation:
2579 * * Each function keeps a local index variable 'idx' that points at the
2580 * current position within the parsed string. The global index is only
2581 * updated if the function could parse the corresponding syntactic unit
2583 * * Similarly, the functions keep track of indicators of inexactness ('#',
2584 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2585 * global exactness information is only updated after each part has been
2586 * successfully parsed.
2587 * * Sequences of digits are parsed into temporary variables holding fixnums.
2588 * Only if these fixnums would overflow, the result variables are updated
2589 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2590 * the temporary variables holding the fixnums are cleared, and the process
2591 * starts over again. If for example fixnums were able to store five decimal
2592 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2593 * and the result was computed as 12345 * 100000 + 67890. In other words,
2594 * only every five digits two bignum operations were performed.
2597 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2599 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2601 /* Caller is responsible for checking that the return value is in range
2602 for the given radix, which should be <= 36. */
2604 char_decimal_value (scm_t_uint32 c
)
2606 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
2607 that's certainly above any valid decimal, so we take advantage of
2608 that to elide some tests. */
2609 unsigned int d
= (unsigned int) uc_decimal_value (c
);
2611 /* If that failed, try extended hexadecimals, then. Only accept ascii
2616 if (c
>= (scm_t_uint32
) 'a')
2617 d
= c
- (scm_t_uint32
)'a' + 10U;
2623 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2624 unsigned int radix
, enum t_exactness
*p_exactness
)
2626 unsigned int idx
= *p_idx
;
2627 unsigned int hash_seen
= 0;
2628 scm_t_bits shift
= 1;
2630 unsigned int digit_value
;
2633 size_t len
= scm_i_string_length (mem
);
2638 c
= scm_i_string_ref (mem
, idx
);
2639 digit_value
= char_decimal_value (c
);
2640 if (digit_value
>= radix
)
2644 result
= SCM_I_MAKINUM (digit_value
);
2647 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2657 digit_value
= char_decimal_value (c
);
2658 /* This check catches non-decimals in addition to out-of-range
2660 if (digit_value
>= radix
)
2665 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2667 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2669 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2676 shift
= shift
* radix
;
2677 add
= add
* radix
+ digit_value
;
2682 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2684 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2688 *p_exactness
= INEXACT
;
2694 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2695 * covers the parts of the rules that start at a potential point. The value
2696 * of the digits up to the point have been parsed by the caller and are given
2697 * in variable result. The content of *p_exactness indicates, whether a hash
2698 * has already been seen in the digits before the point.
2701 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2704 mem2decimal_from_point (SCM result
, SCM mem
,
2705 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2707 unsigned int idx
= *p_idx
;
2708 enum t_exactness x
= *p_exactness
;
2709 size_t len
= scm_i_string_length (mem
);
2714 if (scm_i_string_ref (mem
, idx
) == '.')
2716 scm_t_bits shift
= 1;
2718 unsigned int digit_value
;
2719 SCM big_shift
= SCM_INUM1
;
2724 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2725 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2730 digit_value
= DIGIT2UINT (c
);
2741 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2743 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2744 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2746 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2754 add
= add
* 10 + digit_value
;
2760 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2761 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2762 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2765 result
= scm_divide (result
, big_shift
);
2767 /* We've seen a decimal point, thus the value is implicitly inexact. */
2779 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2781 switch (scm_i_string_ref (mem
, idx
))
2793 c
= scm_i_string_ref (mem
, idx
);
2801 c
= scm_i_string_ref (mem
, idx
);
2810 c
= scm_i_string_ref (mem
, idx
);
2815 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2819 exponent
= DIGIT2UINT (c
);
2822 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2823 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2826 if (exponent
<= SCM_MAXEXP
)
2827 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2833 if (exponent
> SCM_MAXEXP
)
2835 size_t exp_len
= idx
- start
;
2836 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2837 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2838 scm_out_of_range ("string->number", exp_num
);
2841 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2843 result
= scm_product (result
, e
);
2845 result
= scm_divide2real (result
, e
);
2847 /* We've seen an exponent, thus the value is implicitly inexact. */
2865 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2868 mem2ureal (SCM mem
, unsigned int *p_idx
,
2869 unsigned int radix
, enum t_exactness
*p_exactness
)
2871 unsigned int idx
= *p_idx
;
2873 size_t len
= scm_i_string_length (mem
);
2875 /* Start off believing that the number will be exact. This changes
2876 to INEXACT if we see a decimal point or a hash. */
2877 enum t_exactness x
= EXACT
;
2882 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2888 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2890 /* Cobble up the fractional part. We might want to set the
2891 NaN's mantissa from it. */
2893 mem2uinteger (mem
, &idx
, 10, &x
);
2898 if (scm_i_string_ref (mem
, idx
) == '.')
2902 else if (idx
+ 1 == len
)
2904 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2907 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
2914 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2915 if (scm_is_false (uinteger
))
2920 else if (scm_i_string_ref (mem
, idx
) == '/')
2928 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2929 if (scm_is_false (divisor
))
2932 /* both are int/big here, I assume */
2933 result
= scm_i_make_ratio (uinteger
, divisor
);
2935 else if (radix
== 10)
2937 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2938 if (scm_is_false (result
))
2947 /* Update *p_exactness if the number just read was inexact. This is
2948 important for complex numbers, so that a complex number is
2949 treated as inexact overall if either its real or imaginary part
2955 /* When returning an inexact zero, make sure it is represented as a
2956 floating point value so that we can change its sign.
2958 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
2959 result
= scm_from_double (0.0);
2965 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2968 mem2complex (SCM mem
, unsigned int idx
,
2969 unsigned int radix
, enum t_exactness
*p_exactness
)
2974 size_t len
= scm_i_string_length (mem
);
2979 c
= scm_i_string_ref (mem
, idx
);
2994 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2995 if (scm_is_false (ureal
))
2997 /* input must be either +i or -i */
3002 if (scm_i_string_ref (mem
, idx
) == 'i'
3003 || scm_i_string_ref (mem
, idx
) == 'I')
3009 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
3016 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3017 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
3022 c
= scm_i_string_ref (mem
, idx
);
3026 /* either +<ureal>i or -<ureal>i */
3033 return scm_make_rectangular (SCM_INUM0
, ureal
);
3036 /* polar input: <real>@<real>. */
3047 c
= scm_i_string_ref (mem
, idx
);
3065 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3066 if (scm_is_false (angle
))
3071 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3072 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3074 result
= scm_make_polar (ureal
, angle
);
3079 /* expecting input matching <real>[+-]<ureal>?i */
3086 int sign
= (c
== '+') ? 1 : -1;
3087 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3089 if (scm_is_false (imag
))
3090 imag
= SCM_I_MAKINUM (sign
);
3091 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
3092 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3096 if (scm_i_string_ref (mem
, idx
) != 'i'
3097 && scm_i_string_ref (mem
, idx
) != 'I')
3104 return scm_make_rectangular (ureal
, imag
);
3113 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3115 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3118 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3120 unsigned int idx
= 0;
3121 unsigned int radix
= NO_RADIX
;
3122 enum t_exactness forced_x
= NO_EXACTNESS
;
3123 enum t_exactness implicit_x
= EXACT
;
3125 size_t len
= scm_i_string_length (mem
);
3127 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3128 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3130 switch (scm_i_string_ref (mem
, idx
+ 1))
3133 if (radix
!= NO_RADIX
)
3138 if (radix
!= NO_RADIX
)
3143 if (forced_x
!= NO_EXACTNESS
)
3148 if (forced_x
!= NO_EXACTNESS
)
3153 if (radix
!= NO_RADIX
)
3158 if (radix
!= NO_RADIX
)
3168 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3169 if (radix
== NO_RADIX
)
3170 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3172 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3174 if (scm_is_false (result
))
3180 if (SCM_INEXACTP (result
))
3181 return scm_inexact_to_exact (result
);
3185 if (SCM_INEXACTP (result
))
3188 return scm_exact_to_inexact (result
);
3191 if (implicit_x
== INEXACT
)
3193 if (SCM_INEXACTP (result
))
3196 return scm_exact_to_inexact (result
);
3204 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3205 unsigned int default_radix
)
3207 SCM str
= scm_from_locale_stringn (mem
, len
);
3209 return scm_i_string_to_number (str
, default_radix
);
3213 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3214 (SCM string
, SCM radix
),
3215 "Return a number of the maximally precise representation\n"
3216 "expressed by the given @var{string}. @var{radix} must be an\n"
3217 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3218 "is a default radix that may be overridden by an explicit radix\n"
3219 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3220 "supplied, then the default radix is 10. If string is not a\n"
3221 "syntactically valid notation for a number, then\n"
3222 "@code{string->number} returns @code{#f}.")
3223 #define FUNC_NAME s_scm_string_to_number
3227 SCM_VALIDATE_STRING (1, string
);
3229 if (SCM_UNBNDP (radix
))
3232 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3234 answer
= scm_i_string_to_number (string
, base
);
3235 scm_remember_upto_here_1 (string
);
3241 /*** END strs->nums ***/
3245 scm_bigequal (SCM x
, SCM y
)
3247 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3248 scm_remember_upto_here_2 (x
, y
);
3249 return scm_from_bool (0 == result
);
3253 scm_real_equalp (SCM x
, SCM y
)
3255 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3259 scm_complex_equalp (SCM x
, SCM y
)
3261 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3262 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3266 scm_i_fraction_equalp (SCM x
, SCM y
)
3268 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3269 SCM_FRACTION_NUMERATOR (y
)))
3270 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3271 SCM_FRACTION_DENOMINATOR (y
))))
3278 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3280 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3282 #define FUNC_NAME s_scm_number_p
3284 return scm_from_bool (SCM_NUMBERP (x
));
3288 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3290 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3291 "otherwise. Note that the sets of real, rational and integer\n"
3292 "values form subsets of the set of complex numbers, i. e. the\n"
3293 "predicate will also be fulfilled if @var{x} is a real,\n"
3294 "rational or integer number.")
3295 #define FUNC_NAME s_scm_complex_p
3297 /* all numbers are complex. */
3298 return scm_number_p (x
);
3302 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3304 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3305 "otherwise. Note that the set of integer values forms a subset of\n"
3306 "the set of real numbers, i. e. the predicate will also be\n"
3307 "fulfilled if @var{x} is an integer number.")
3308 #define FUNC_NAME s_scm_real_p
3310 /* we can't represent irrational numbers. */
3311 return scm_rational_p (x
);
3315 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3317 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3318 "otherwise. Note that the set of integer values forms a subset of\n"
3319 "the set of rational numbers, i. e. the predicate will also be\n"
3320 "fulfilled if @var{x} is an integer number.")
3321 #define FUNC_NAME s_scm_rational_p
3323 if (SCM_I_INUMP (x
))
3325 else if (SCM_IMP (x
))
3327 else if (SCM_BIGP (x
))
3329 else if (SCM_FRACTIONP (x
))
3331 else if (SCM_REALP (x
))
3332 /* due to their limited precision, all floating point numbers are
3333 rational as well. */
3340 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3342 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3344 #define FUNC_NAME s_scm_integer_p
3347 if (SCM_I_INUMP (x
))
3353 if (!SCM_INEXACTP (x
))
3355 if (SCM_COMPLEXP (x
))
3357 r
= SCM_REAL_VALUE (x
);
3367 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3369 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3371 #define FUNC_NAME s_scm_inexact_p
3373 if (SCM_INEXACTP (x
))
3375 if (SCM_NUMBERP (x
))
3377 SCM_WRONG_TYPE_ARG (1, x
);
3382 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3383 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3384 (SCM x
, SCM y
, SCM rest
),
3385 "Return @code{#t} if all parameters are numerically equal.")
3386 #define FUNC_NAME s_scm_i_num_eq_p
3388 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3390 while (!scm_is_null (rest
))
3392 if (scm_is_false (scm_num_eq_p (x
, y
)))
3396 rest
= scm_cdr (rest
);
3398 return scm_num_eq_p (x
, y
);
3402 scm_num_eq_p (SCM x
, SCM y
)
3405 if (SCM_I_INUMP (x
))
3407 scm_t_signed_bits xx
= SCM_I_INUM (x
);
3408 if (SCM_I_INUMP (y
))
3410 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3411 return scm_from_bool (xx
== yy
);
3413 else if (SCM_BIGP (y
))
3415 else if (SCM_REALP (y
))
3417 /* On a 32-bit system an inum fits a double, we can cast the inum
3418 to a double and compare.
3420 But on a 64-bit system an inum is bigger than a double and
3421 casting it to a double (call that dxx) will round. dxx is at
3422 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3423 an integer and fits a long. So we cast yy to a long and
3424 compare with plain xx.
3426 An alternative (for any size system actually) would be to check
3427 yy is an integer (with floor) and is in range of an inum
3428 (compare against appropriate powers of 2) then test
3429 xx==(scm_t_signed_bits)yy. It's just a matter of which
3430 casts/comparisons might be fastest or easiest for the cpu. */
3432 double yy
= SCM_REAL_VALUE (y
);
3433 return scm_from_bool ((double) xx
== yy
3434 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3435 || xx
== (scm_t_signed_bits
) yy
));
3437 else if (SCM_COMPLEXP (y
))
3438 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3439 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3440 else if (SCM_FRACTIONP (y
))
3443 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3445 else if (SCM_BIGP (x
))
3447 if (SCM_I_INUMP (y
))
3449 else if (SCM_BIGP (y
))
3451 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3452 scm_remember_upto_here_2 (x
, y
);
3453 return scm_from_bool (0 == cmp
);
3455 else if (SCM_REALP (y
))
3458 if (isnan (SCM_REAL_VALUE (y
)))
3460 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3461 scm_remember_upto_here_1 (x
);
3462 return scm_from_bool (0 == cmp
);
3464 else if (SCM_COMPLEXP (y
))
3467 if (0.0 != SCM_COMPLEX_IMAG (y
))
3469 if (isnan (SCM_COMPLEX_REAL (y
)))
3471 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3472 scm_remember_upto_here_1 (x
);
3473 return scm_from_bool (0 == cmp
);
3475 else if (SCM_FRACTIONP (y
))
3478 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3480 else if (SCM_REALP (x
))
3482 double xx
= SCM_REAL_VALUE (x
);
3483 if (SCM_I_INUMP (y
))
3485 /* see comments with inum/real above */
3486 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3487 return scm_from_bool (xx
== (double) yy
3488 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3489 || (scm_t_signed_bits
) xx
== yy
));
3491 else if (SCM_BIGP (y
))
3494 if (isnan (SCM_REAL_VALUE (x
)))
3496 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3497 scm_remember_upto_here_1 (y
);
3498 return scm_from_bool (0 == cmp
);
3500 else if (SCM_REALP (y
))
3501 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3502 else if (SCM_COMPLEXP (y
))
3503 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3504 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3505 else if (SCM_FRACTIONP (y
))
3507 double xx
= SCM_REAL_VALUE (x
);
3511 return scm_from_bool (xx
< 0.0);
3512 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3516 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3518 else if (SCM_COMPLEXP (x
))
3520 if (SCM_I_INUMP (y
))
3521 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3522 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3523 else if (SCM_BIGP (y
))
3526 if (0.0 != SCM_COMPLEX_IMAG (x
))
3528 if (isnan (SCM_COMPLEX_REAL (x
)))
3530 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3531 scm_remember_upto_here_1 (y
);
3532 return scm_from_bool (0 == cmp
);
3534 else if (SCM_REALP (y
))
3535 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3536 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3537 else if (SCM_COMPLEXP (y
))
3538 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3539 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3540 else if (SCM_FRACTIONP (y
))
3543 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3545 xx
= SCM_COMPLEX_REAL (x
);
3549 return scm_from_bool (xx
< 0.0);
3550 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3554 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3556 else if (SCM_FRACTIONP (x
))
3558 if (SCM_I_INUMP (y
))
3560 else if (SCM_BIGP (y
))
3562 else if (SCM_REALP (y
))
3564 double yy
= SCM_REAL_VALUE (y
);
3568 return scm_from_bool (0.0 < yy
);
3569 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3572 else if (SCM_COMPLEXP (y
))
3575 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3577 yy
= SCM_COMPLEX_REAL (y
);
3581 return scm_from_bool (0.0 < yy
);
3582 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3585 else if (SCM_FRACTIONP (y
))
3586 return scm_i_fraction_equalp (x
, y
);
3588 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3591 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3595 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3596 done are good for inums, but for bignums an answer can almost always be
3597 had by just examining a few high bits of the operands, as done by GMP in
3598 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3599 of the float exponent to take into account. */
3601 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3602 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3603 (SCM x
, SCM y
, SCM rest
),
3604 "Return @code{#t} if the list of parameters is monotonically\n"
3606 #define FUNC_NAME s_scm_i_num_less_p
3608 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3610 while (!scm_is_null (rest
))
3612 if (scm_is_false (scm_less_p (x
, y
)))
3616 rest
= scm_cdr (rest
);
3618 return scm_less_p (x
, y
);
3622 scm_less_p (SCM x
, SCM y
)
3625 if (SCM_I_INUMP (x
))
3627 scm_t_inum xx
= SCM_I_INUM (x
);
3628 if (SCM_I_INUMP (y
))
3630 scm_t_inum yy
= SCM_I_INUM (y
);
3631 return scm_from_bool (xx
< yy
);
3633 else if (SCM_BIGP (y
))
3635 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3636 scm_remember_upto_here_1 (y
);
3637 return scm_from_bool (sgn
> 0);
3639 else if (SCM_REALP (y
))
3640 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3641 else if (SCM_FRACTIONP (y
))
3643 /* "x < a/b" becomes "x*b < a" */
3645 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3646 y
= SCM_FRACTION_NUMERATOR (y
);
3650 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3652 else if (SCM_BIGP (x
))
3654 if (SCM_I_INUMP (y
))
3656 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3657 scm_remember_upto_here_1 (x
);
3658 return scm_from_bool (sgn
< 0);
3660 else if (SCM_BIGP (y
))
3662 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3663 scm_remember_upto_here_2 (x
, y
);
3664 return scm_from_bool (cmp
< 0);
3666 else if (SCM_REALP (y
))
3669 if (isnan (SCM_REAL_VALUE (y
)))
3671 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3672 scm_remember_upto_here_1 (x
);
3673 return scm_from_bool (cmp
< 0);
3675 else if (SCM_FRACTIONP (y
))
3678 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3680 else if (SCM_REALP (x
))
3682 if (SCM_I_INUMP (y
))
3683 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3684 else if (SCM_BIGP (y
))
3687 if (isnan (SCM_REAL_VALUE (x
)))
3689 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3690 scm_remember_upto_here_1 (y
);
3691 return scm_from_bool (cmp
> 0);
3693 else if (SCM_REALP (y
))
3694 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3695 else if (SCM_FRACTIONP (y
))
3697 double xx
= SCM_REAL_VALUE (x
);
3701 return scm_from_bool (xx
< 0.0);
3702 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3706 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3708 else if (SCM_FRACTIONP (x
))
3710 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3712 /* "a/b < y" becomes "a < y*b" */
3713 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3714 x
= SCM_FRACTION_NUMERATOR (x
);
3717 else if (SCM_REALP (y
))
3719 double yy
= SCM_REAL_VALUE (y
);
3723 return scm_from_bool (0.0 < yy
);
3724 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3727 else if (SCM_FRACTIONP (y
))
3729 /* "a/b < c/d" becomes "a*d < c*b" */
3730 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3731 SCM_FRACTION_DENOMINATOR (y
));
3732 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3733 SCM_FRACTION_DENOMINATOR (x
));
3739 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3742 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3746 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3747 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3748 (SCM x
, SCM y
, SCM rest
),
3749 "Return @code{#t} if the list of parameters is monotonically\n"
3751 #define FUNC_NAME s_scm_i_num_gr_p
3753 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3755 while (!scm_is_null (rest
))
3757 if (scm_is_false (scm_gr_p (x
, y
)))
3761 rest
= scm_cdr (rest
);
3763 return scm_gr_p (x
, y
);
3766 #define FUNC_NAME s_scm_i_num_gr_p
3768 scm_gr_p (SCM x
, SCM y
)
3770 if (!SCM_NUMBERP (x
))
3771 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3772 else if (!SCM_NUMBERP (y
))
3773 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3775 return scm_less_p (y
, x
);
3780 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3781 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3782 (SCM x
, SCM y
, SCM rest
),
3783 "Return @code{#t} if the list of parameters is monotonically\n"
3785 #define FUNC_NAME s_scm_i_num_leq_p
3787 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3789 while (!scm_is_null (rest
))
3791 if (scm_is_false (scm_leq_p (x
, y
)))
3795 rest
= scm_cdr (rest
);
3797 return scm_leq_p (x
, y
);
3800 #define FUNC_NAME s_scm_i_num_leq_p
3802 scm_leq_p (SCM x
, SCM y
)
3804 if (!SCM_NUMBERP (x
))
3805 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3806 else if (!SCM_NUMBERP (y
))
3807 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3808 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3811 return scm_not (scm_less_p (y
, x
));
3816 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3817 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3818 (SCM x
, SCM y
, SCM rest
),
3819 "Return @code{#t} if the list of parameters is monotonically\n"
3821 #define FUNC_NAME s_scm_i_num_geq_p
3823 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3825 while (!scm_is_null (rest
))
3827 if (scm_is_false (scm_geq_p (x
, y
)))
3831 rest
= scm_cdr (rest
);
3833 return scm_geq_p (x
, y
);
3836 #define FUNC_NAME s_scm_i_num_geq_p
3838 scm_geq_p (SCM x
, SCM y
)
3840 if (!SCM_NUMBERP (x
))
3841 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3842 else if (!SCM_NUMBERP (y
))
3843 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3844 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3847 return scm_not (scm_less_p (x
, y
));
3852 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3853 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3859 if (SCM_I_INUMP (z
))
3860 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3861 else if (SCM_BIGP (z
))
3863 else if (SCM_REALP (z
))
3864 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3865 else if (SCM_COMPLEXP (z
))
3866 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3867 && SCM_COMPLEX_IMAG (z
) == 0.0);
3868 else if (SCM_FRACTIONP (z
))
3871 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3875 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3876 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3880 scm_positive_p (SCM x
)
3882 if (SCM_I_INUMP (x
))
3883 return scm_from_bool (SCM_I_INUM (x
) > 0);
3884 else if (SCM_BIGP (x
))
3886 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3887 scm_remember_upto_here_1 (x
);
3888 return scm_from_bool (sgn
> 0);
3890 else if (SCM_REALP (x
))
3891 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3892 else if (SCM_FRACTIONP (x
))
3893 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3895 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3899 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3900 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3904 scm_negative_p (SCM x
)
3906 if (SCM_I_INUMP (x
))
3907 return scm_from_bool (SCM_I_INUM (x
) < 0);
3908 else if (SCM_BIGP (x
))
3910 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3911 scm_remember_upto_here_1 (x
);
3912 return scm_from_bool (sgn
< 0);
3914 else if (SCM_REALP (x
))
3915 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3916 else if (SCM_FRACTIONP (x
))
3917 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3919 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3923 /* scm_min and scm_max return an inexact when either argument is inexact, as
3924 required by r5rs. On that basis, for exact/inexact combinations the
3925 exact is converted to inexact to compare and possibly return. This is
3926 unlike scm_less_p above which takes some trouble to preserve all bits in
3927 its test, such trouble is not required for min and max. */
3929 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3930 (SCM x
, SCM y
, SCM rest
),
3931 "Return the maximum of all parameter values.")
3932 #define FUNC_NAME s_scm_i_max
3934 while (!scm_is_null (rest
))
3935 { x
= scm_max (x
, y
);
3937 rest
= scm_cdr (rest
);
3939 return scm_max (x
, y
);
3943 #define s_max s_scm_i_max
3944 #define g_max g_scm_i_max
3947 scm_max (SCM x
, SCM y
)
3952 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3953 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3956 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3959 if (SCM_I_INUMP (x
))
3961 scm_t_inum xx
= SCM_I_INUM (x
);
3962 if (SCM_I_INUMP (y
))
3964 scm_t_inum yy
= SCM_I_INUM (y
);
3965 return (xx
< yy
) ? y
: x
;
3967 else if (SCM_BIGP (y
))
3969 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3970 scm_remember_upto_here_1 (y
);
3971 return (sgn
< 0) ? x
: y
;
3973 else if (SCM_REALP (y
))
3976 /* if y==NaN then ">" is false and we return NaN */
3977 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3979 else if (SCM_FRACTIONP (y
))
3982 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3985 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3987 else if (SCM_BIGP (x
))
3989 if (SCM_I_INUMP (y
))
3991 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3992 scm_remember_upto_here_1 (x
);
3993 return (sgn
< 0) ? y
: x
;
3995 else if (SCM_BIGP (y
))
3997 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3998 scm_remember_upto_here_2 (x
, y
);
3999 return (cmp
> 0) ? x
: y
;
4001 else if (SCM_REALP (y
))
4003 /* if y==NaN then xx>yy is false, so we return the NaN y */
4006 xx
= scm_i_big2dbl (x
);
4007 yy
= SCM_REAL_VALUE (y
);
4008 return (xx
> yy
? scm_from_double (xx
) : y
);
4010 else if (SCM_FRACTIONP (y
))
4015 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4017 else if (SCM_REALP (x
))
4019 if (SCM_I_INUMP (y
))
4021 double z
= SCM_I_INUM (y
);
4022 /* if x==NaN then "<" is false and we return NaN */
4023 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
4025 else if (SCM_BIGP (y
))
4030 else if (SCM_REALP (y
))
4032 /* if x==NaN then our explicit check means we return NaN
4033 if y==NaN then ">" is false and we return NaN
4034 calling isnan is unavoidable, since it's the only way to know
4035 which of x or y causes any compares to be false */
4036 double xx
= SCM_REAL_VALUE (x
);
4037 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
4039 else if (SCM_FRACTIONP (y
))
4041 double yy
= scm_i_fraction2double (y
);
4042 double xx
= SCM_REAL_VALUE (x
);
4043 return (xx
< yy
) ? scm_from_double (yy
) : x
;
4046 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4048 else if (SCM_FRACTIONP (x
))
4050 if (SCM_I_INUMP (y
))
4054 else if (SCM_BIGP (y
))
4058 else if (SCM_REALP (y
))
4060 double xx
= scm_i_fraction2double (x
);
4061 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4063 else if (SCM_FRACTIONP (y
))
4068 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4071 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4075 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4076 (SCM x
, SCM y
, SCM rest
),
4077 "Return the minimum of all parameter values.")
4078 #define FUNC_NAME s_scm_i_min
4080 while (!scm_is_null (rest
))
4081 { x
= scm_min (x
, y
);
4083 rest
= scm_cdr (rest
);
4085 return scm_min (x
, y
);
4089 #define s_min s_scm_i_min
4090 #define g_min g_scm_i_min
4093 scm_min (SCM x
, SCM y
)
4098 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4099 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4102 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4105 if (SCM_I_INUMP (x
))
4107 scm_t_inum xx
= SCM_I_INUM (x
);
4108 if (SCM_I_INUMP (y
))
4110 scm_t_inum yy
= SCM_I_INUM (y
);
4111 return (xx
< yy
) ? x
: y
;
4113 else if (SCM_BIGP (y
))
4115 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4116 scm_remember_upto_here_1 (y
);
4117 return (sgn
< 0) ? y
: x
;
4119 else if (SCM_REALP (y
))
4122 /* if y==NaN then "<" is false and we return NaN */
4123 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4125 else if (SCM_FRACTIONP (y
))
4128 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4131 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4133 else if (SCM_BIGP (x
))
4135 if (SCM_I_INUMP (y
))
4137 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4138 scm_remember_upto_here_1 (x
);
4139 return (sgn
< 0) ? x
: y
;
4141 else if (SCM_BIGP (y
))
4143 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4144 scm_remember_upto_here_2 (x
, y
);
4145 return (cmp
> 0) ? y
: x
;
4147 else if (SCM_REALP (y
))
4149 /* if y==NaN then xx<yy is false, so we return the NaN y */
4152 xx
= scm_i_big2dbl (x
);
4153 yy
= SCM_REAL_VALUE (y
);
4154 return (xx
< yy
? scm_from_double (xx
) : y
);
4156 else if (SCM_FRACTIONP (y
))
4161 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4163 else if (SCM_REALP (x
))
4165 if (SCM_I_INUMP (y
))
4167 double z
= SCM_I_INUM (y
);
4168 /* if x==NaN then "<" is false and we return NaN */
4169 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4171 else if (SCM_BIGP (y
))
4176 else if (SCM_REALP (y
))
4178 /* if x==NaN then our explicit check means we return NaN
4179 if y==NaN then "<" is false and we return NaN
4180 calling isnan is unavoidable, since it's the only way to know
4181 which of x or y causes any compares to be false */
4182 double xx
= SCM_REAL_VALUE (x
);
4183 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4185 else if (SCM_FRACTIONP (y
))
4187 double yy
= scm_i_fraction2double (y
);
4188 double xx
= SCM_REAL_VALUE (x
);
4189 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4192 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4194 else if (SCM_FRACTIONP (x
))
4196 if (SCM_I_INUMP (y
))
4200 else if (SCM_BIGP (y
))
4204 else if (SCM_REALP (y
))
4206 double xx
= scm_i_fraction2double (x
);
4207 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4209 else if (SCM_FRACTIONP (y
))
4214 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4217 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4221 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4222 (SCM x
, SCM y
, SCM rest
),
4223 "Return the sum of all parameter values. Return 0 if called without\n"
4225 #define FUNC_NAME s_scm_i_sum
4227 while (!scm_is_null (rest
))
4228 { x
= scm_sum (x
, y
);
4230 rest
= scm_cdr (rest
);
4232 return scm_sum (x
, y
);
4236 #define s_sum s_scm_i_sum
4237 #define g_sum g_scm_i_sum
4240 scm_sum (SCM x
, SCM y
)
4242 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4244 if (SCM_NUMBERP (x
)) return x
;
4245 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4246 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4249 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4251 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4253 scm_t_inum xx
= SCM_I_INUM (x
);
4254 scm_t_inum yy
= SCM_I_INUM (y
);
4255 scm_t_inum z
= xx
+ yy
;
4256 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
4258 else if (SCM_BIGP (y
))
4263 else if (SCM_REALP (y
))
4265 scm_t_inum xx
= SCM_I_INUM (x
);
4266 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4268 else if (SCM_COMPLEXP (y
))
4270 scm_t_inum xx
= SCM_I_INUM (x
);
4271 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4272 SCM_COMPLEX_IMAG (y
));
4274 else if (SCM_FRACTIONP (y
))
4275 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4276 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4277 SCM_FRACTION_DENOMINATOR (y
));
4279 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4280 } else if (SCM_BIGP (x
))
4282 if (SCM_I_INUMP (y
))
4287 inum
= SCM_I_INUM (y
);
4290 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4293 SCM result
= scm_i_mkbig ();
4294 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4295 scm_remember_upto_here_1 (x
);
4296 /* we know the result will have to be a bignum */
4299 return scm_i_normbig (result
);
4303 SCM result
= scm_i_mkbig ();
4304 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4305 scm_remember_upto_here_1 (x
);
4306 /* we know the result will have to be a bignum */
4309 return scm_i_normbig (result
);
4312 else if (SCM_BIGP (y
))
4314 SCM result
= scm_i_mkbig ();
4315 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4316 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4317 mpz_add (SCM_I_BIG_MPZ (result
),
4320 scm_remember_upto_here_2 (x
, y
);
4321 /* we know the result will have to be a bignum */
4324 return scm_i_normbig (result
);
4326 else if (SCM_REALP (y
))
4328 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4329 scm_remember_upto_here_1 (x
);
4330 return scm_from_double (result
);
4332 else if (SCM_COMPLEXP (y
))
4334 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4335 + SCM_COMPLEX_REAL (y
));
4336 scm_remember_upto_here_1 (x
);
4337 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4339 else if (SCM_FRACTIONP (y
))
4340 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4341 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4342 SCM_FRACTION_DENOMINATOR (y
));
4344 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4346 else if (SCM_REALP (x
))
4348 if (SCM_I_INUMP (y
))
4349 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4350 else if (SCM_BIGP (y
))
4352 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4353 scm_remember_upto_here_1 (y
);
4354 return scm_from_double (result
);
4356 else if (SCM_REALP (y
))
4357 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4358 else if (SCM_COMPLEXP (y
))
4359 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4360 SCM_COMPLEX_IMAG (y
));
4361 else if (SCM_FRACTIONP (y
))
4362 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4364 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4366 else if (SCM_COMPLEXP (x
))
4368 if (SCM_I_INUMP (y
))
4369 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4370 SCM_COMPLEX_IMAG (x
));
4371 else if (SCM_BIGP (y
))
4373 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4374 + SCM_COMPLEX_REAL (x
));
4375 scm_remember_upto_here_1 (y
);
4376 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4378 else if (SCM_REALP (y
))
4379 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4380 SCM_COMPLEX_IMAG (x
));
4381 else if (SCM_COMPLEXP (y
))
4382 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4383 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4384 else if (SCM_FRACTIONP (y
))
4385 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4386 SCM_COMPLEX_IMAG (x
));
4388 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4390 else if (SCM_FRACTIONP (x
))
4392 if (SCM_I_INUMP (y
))
4393 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4394 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4395 SCM_FRACTION_DENOMINATOR (x
));
4396 else if (SCM_BIGP (y
))
4397 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4398 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4399 SCM_FRACTION_DENOMINATOR (x
));
4400 else if (SCM_REALP (y
))
4401 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4402 else if (SCM_COMPLEXP (y
))
4403 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4404 SCM_COMPLEX_IMAG (y
));
4405 else if (SCM_FRACTIONP (y
))
4406 /* a/b + c/d = (ad + bc) / bd */
4407 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4408 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4409 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4411 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4414 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4418 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4420 "Return @math{@var{x}+1}.")
4421 #define FUNC_NAME s_scm_oneplus
4423 return scm_sum (x
, SCM_INUM1
);
4428 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4429 (SCM x
, SCM y
, SCM rest
),
4430 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4431 "the sum of all but the first argument are subtracted from the first\n"
4433 #define FUNC_NAME s_scm_i_difference
4435 while (!scm_is_null (rest
))
4436 { x
= scm_difference (x
, y
);
4438 rest
= scm_cdr (rest
);
4440 return scm_difference (x
, y
);
4444 #define s_difference s_scm_i_difference
4445 #define g_difference g_scm_i_difference
4448 scm_difference (SCM x
, SCM y
)
4449 #define FUNC_NAME s_difference
4451 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4454 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4456 if (SCM_I_INUMP (x
))
4458 scm_t_inum xx
= -SCM_I_INUM (x
);
4459 if (SCM_FIXABLE (xx
))
4460 return SCM_I_MAKINUM (xx
);
4462 return scm_i_inum2big (xx
);
4464 else if (SCM_BIGP (x
))
4465 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4466 bignum, but negating that gives a fixnum. */
4467 return scm_i_normbig (scm_i_clonebig (x
, 0));
4468 else if (SCM_REALP (x
))
4469 return scm_from_double (-SCM_REAL_VALUE (x
));
4470 else if (SCM_COMPLEXP (x
))
4471 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4472 -SCM_COMPLEX_IMAG (x
));
4473 else if (SCM_FRACTIONP (x
))
4474 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4475 SCM_FRACTION_DENOMINATOR (x
));
4477 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4480 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4482 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4484 scm_t_inum xx
= SCM_I_INUM (x
);
4485 scm_t_inum yy
= SCM_I_INUM (y
);
4486 scm_t_inum z
= xx
- yy
;
4487 if (SCM_FIXABLE (z
))
4488 return SCM_I_MAKINUM (z
);
4490 return scm_i_inum2big (z
);
4492 else if (SCM_BIGP (y
))
4494 /* inum-x - big-y */
4495 scm_t_inum xx
= SCM_I_INUM (x
);
4498 return scm_i_clonebig (y
, 0);
4501 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4502 SCM result
= scm_i_mkbig ();
4505 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4508 /* x - y == -(y + -x) */
4509 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4510 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4512 scm_remember_upto_here_1 (y
);
4514 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4515 /* we know the result will have to be a bignum */
4518 return scm_i_normbig (result
);
4521 else if (SCM_REALP (y
))
4523 scm_t_inum xx
= SCM_I_INUM (x
);
4524 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4526 else if (SCM_COMPLEXP (y
))
4528 scm_t_inum xx
= SCM_I_INUM (x
);
4529 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4530 - SCM_COMPLEX_IMAG (y
));
4532 else if (SCM_FRACTIONP (y
))
4533 /* a - b/c = (ac - b) / c */
4534 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4535 SCM_FRACTION_NUMERATOR (y
)),
4536 SCM_FRACTION_DENOMINATOR (y
));
4538 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4540 else if (SCM_BIGP (x
))
4542 if (SCM_I_INUMP (y
))
4544 /* big-x - inum-y */
4545 scm_t_inum yy
= SCM_I_INUM (y
);
4546 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4548 scm_remember_upto_here_1 (x
);
4550 return (SCM_FIXABLE (-yy
) ?
4551 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
4554 SCM result
= scm_i_mkbig ();
4557 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4559 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4560 scm_remember_upto_here_1 (x
);
4562 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4563 /* we know the result will have to be a bignum */
4566 return scm_i_normbig (result
);
4569 else if (SCM_BIGP (y
))
4571 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4572 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4573 SCM result
= scm_i_mkbig ();
4574 mpz_sub (SCM_I_BIG_MPZ (result
),
4577 scm_remember_upto_here_2 (x
, y
);
4578 /* we know the result will have to be a bignum */
4579 if ((sgn_x
== 1) && (sgn_y
== -1))
4581 if ((sgn_x
== -1) && (sgn_y
== 1))
4583 return scm_i_normbig (result
);
4585 else if (SCM_REALP (y
))
4587 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4588 scm_remember_upto_here_1 (x
);
4589 return scm_from_double (result
);
4591 else if (SCM_COMPLEXP (y
))
4593 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4594 - SCM_COMPLEX_REAL (y
));
4595 scm_remember_upto_here_1 (x
);
4596 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4598 else if (SCM_FRACTIONP (y
))
4599 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4600 SCM_FRACTION_NUMERATOR (y
)),
4601 SCM_FRACTION_DENOMINATOR (y
));
4602 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4604 else if (SCM_REALP (x
))
4606 if (SCM_I_INUMP (y
))
4607 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4608 else if (SCM_BIGP (y
))
4610 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4611 scm_remember_upto_here_1 (x
);
4612 return scm_from_double (result
);
4614 else if (SCM_REALP (y
))
4615 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4616 else if (SCM_COMPLEXP (y
))
4617 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4618 -SCM_COMPLEX_IMAG (y
));
4619 else if (SCM_FRACTIONP (y
))
4620 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4622 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4624 else if (SCM_COMPLEXP (x
))
4626 if (SCM_I_INUMP (y
))
4627 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4628 SCM_COMPLEX_IMAG (x
));
4629 else if (SCM_BIGP (y
))
4631 double real_part
= (SCM_COMPLEX_REAL (x
)
4632 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4633 scm_remember_upto_here_1 (x
);
4634 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4636 else if (SCM_REALP (y
))
4637 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4638 SCM_COMPLEX_IMAG (x
));
4639 else if (SCM_COMPLEXP (y
))
4640 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4641 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4642 else if (SCM_FRACTIONP (y
))
4643 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4644 SCM_COMPLEX_IMAG (x
));
4646 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4648 else if (SCM_FRACTIONP (x
))
4650 if (SCM_I_INUMP (y
))
4651 /* a/b - c = (a - cb) / b */
4652 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4653 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4654 SCM_FRACTION_DENOMINATOR (x
));
4655 else if (SCM_BIGP (y
))
4656 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4657 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4658 SCM_FRACTION_DENOMINATOR (x
));
4659 else if (SCM_REALP (y
))
4660 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4661 else if (SCM_COMPLEXP (y
))
4662 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4663 -SCM_COMPLEX_IMAG (y
));
4664 else if (SCM_FRACTIONP (y
))
4665 /* a/b - c/d = (ad - bc) / bd */
4666 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4667 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4668 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4670 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4673 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4678 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4680 "Return @math{@var{x}-1}.")
4681 #define FUNC_NAME s_scm_oneminus
4683 return scm_difference (x
, SCM_INUM1
);
4688 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4689 (SCM x
, SCM y
, SCM rest
),
4690 "Return the product of all arguments. If called without arguments,\n"
4692 #define FUNC_NAME s_scm_i_product
4694 while (!scm_is_null (rest
))
4695 { x
= scm_product (x
, y
);
4697 rest
= scm_cdr (rest
);
4699 return scm_product (x
, y
);
4703 #define s_product s_scm_i_product
4704 #define g_product g_scm_i_product
4707 scm_product (SCM x
, SCM y
)
4709 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4712 return SCM_I_MAKINUM (1L);
4713 else if (SCM_NUMBERP (x
))
4716 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4719 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4724 xx
= SCM_I_INUM (x
);
4728 case 0: return x
; break;
4729 case 1: return y
; break;
4732 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4734 scm_t_inum yy
= SCM_I_INUM (y
);
4735 scm_t_inum kk
= xx
* yy
;
4736 SCM k
= SCM_I_MAKINUM (kk
);
4737 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4741 SCM result
= scm_i_inum2big (xx
);
4742 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4743 return scm_i_normbig (result
);
4746 else if (SCM_BIGP (y
))
4748 SCM result
= scm_i_mkbig ();
4749 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4750 scm_remember_upto_here_1 (y
);
4753 else if (SCM_REALP (y
))
4754 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4755 else if (SCM_COMPLEXP (y
))
4756 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4757 xx
* SCM_COMPLEX_IMAG (y
));
4758 else if (SCM_FRACTIONP (y
))
4759 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4760 SCM_FRACTION_DENOMINATOR (y
));
4762 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4764 else if (SCM_BIGP (x
))
4766 if (SCM_I_INUMP (y
))
4771 else if (SCM_BIGP (y
))
4773 SCM result
= scm_i_mkbig ();
4774 mpz_mul (SCM_I_BIG_MPZ (result
),
4777 scm_remember_upto_here_2 (x
, y
);
4780 else if (SCM_REALP (y
))
4782 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4783 scm_remember_upto_here_1 (x
);
4784 return scm_from_double (result
);
4786 else if (SCM_COMPLEXP (y
))
4788 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4789 scm_remember_upto_here_1 (x
);
4790 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4791 z
* SCM_COMPLEX_IMAG (y
));
4793 else if (SCM_FRACTIONP (y
))
4794 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4795 SCM_FRACTION_DENOMINATOR (y
));
4797 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4799 else if (SCM_REALP (x
))
4801 if (SCM_I_INUMP (y
))
4803 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4804 if (scm_is_eq (y
, SCM_INUM0
))
4806 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4808 else if (SCM_BIGP (y
))
4810 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4811 scm_remember_upto_here_1 (y
);
4812 return scm_from_double (result
);
4814 else if (SCM_REALP (y
))
4815 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4816 else if (SCM_COMPLEXP (y
))
4817 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4818 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4819 else if (SCM_FRACTIONP (y
))
4820 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4822 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4824 else if (SCM_COMPLEXP (x
))
4826 if (SCM_I_INUMP (y
))
4828 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4829 if (scm_is_eq (y
, SCM_INUM0
))
4831 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4832 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4834 else if (SCM_BIGP (y
))
4836 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4837 scm_remember_upto_here_1 (y
);
4838 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4839 z
* SCM_COMPLEX_IMAG (x
));
4841 else if (SCM_REALP (y
))
4842 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4843 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4844 else if (SCM_COMPLEXP (y
))
4846 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4847 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4848 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4849 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4851 else if (SCM_FRACTIONP (y
))
4853 double yy
= scm_i_fraction2double (y
);
4854 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4855 yy
* SCM_COMPLEX_IMAG (x
));
4858 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4860 else if (SCM_FRACTIONP (x
))
4862 if (SCM_I_INUMP (y
))
4863 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4864 SCM_FRACTION_DENOMINATOR (x
));
4865 else if (SCM_BIGP (y
))
4866 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4867 SCM_FRACTION_DENOMINATOR (x
));
4868 else if (SCM_REALP (y
))
4869 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4870 else if (SCM_COMPLEXP (y
))
4872 double xx
= scm_i_fraction2double (x
);
4873 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4874 xx
* SCM_COMPLEX_IMAG (y
));
4876 else if (SCM_FRACTIONP (y
))
4877 /* a/b * c/d = ac / bd */
4878 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4879 SCM_FRACTION_NUMERATOR (y
)),
4880 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4881 SCM_FRACTION_DENOMINATOR (y
)));
4883 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4886 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4889 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4890 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4891 #define ALLOW_DIVIDE_BY_ZERO
4892 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4895 /* The code below for complex division is adapted from the GNU
4896 libstdc++, which adapted it from f2c's libF77, and is subject to
4899 /****************************************************************
4900 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4902 Permission to use, copy, modify, and distribute this software
4903 and its documentation for any purpose and without fee is hereby
4904 granted, provided that the above copyright notice appear in all
4905 copies and that both that the copyright notice and this
4906 permission notice and warranty disclaimer appear in supporting
4907 documentation, and that the names of AT&T Bell Laboratories or
4908 Bellcore or any of their entities not be used in advertising or
4909 publicity pertaining to distribution of the software without
4910 specific, written prior permission.
4912 AT&T and Bellcore disclaim all warranties with regard to this
4913 software, including all implied warranties of merchantability
4914 and fitness. In no event shall AT&T or Bellcore be liable for
4915 any special, indirect or consequential damages or any damages
4916 whatsoever resulting from loss of use, data or profits, whether
4917 in an action of contract, negligence or other tortious action,
4918 arising out of or in connection with the use or performance of
4920 ****************************************************************/
4922 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4923 (SCM x
, SCM y
, SCM rest
),
4924 "Divide the first argument by the product of the remaining\n"
4925 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4927 #define FUNC_NAME s_scm_i_divide
4929 while (!scm_is_null (rest
))
4930 { x
= scm_divide (x
, y
);
4932 rest
= scm_cdr (rest
);
4934 return scm_divide (x
, y
);
4938 #define s_divide s_scm_i_divide
4939 #define g_divide g_scm_i_divide
4942 do_divide (SCM x
, SCM y
, int inexact
)
4943 #define FUNC_NAME s_divide
4947 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4950 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4951 else if (SCM_I_INUMP (x
))
4953 scm_t_inum xx
= SCM_I_INUM (x
);
4954 if (xx
== 1 || xx
== -1)
4956 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4958 scm_num_overflow (s_divide
);
4963 return scm_from_double (1.0 / (double) xx
);
4964 else return scm_i_make_ratio (SCM_INUM1
, x
);
4967 else if (SCM_BIGP (x
))
4970 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4971 else return scm_i_make_ratio (SCM_INUM1
, x
);
4973 else if (SCM_REALP (x
))
4975 double xx
= SCM_REAL_VALUE (x
);
4976 #ifndef ALLOW_DIVIDE_BY_ZERO
4978 scm_num_overflow (s_divide
);
4981 return scm_from_double (1.0 / xx
);
4983 else if (SCM_COMPLEXP (x
))
4985 double r
= SCM_COMPLEX_REAL (x
);
4986 double i
= SCM_COMPLEX_IMAG (x
);
4987 if (fabs(r
) <= fabs(i
))
4990 double d
= i
* (1.0 + t
* t
);
4991 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4996 double d
= r
* (1.0 + t
* t
);
4997 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
5000 else if (SCM_FRACTIONP (x
))
5001 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
5002 SCM_FRACTION_NUMERATOR (x
));
5004 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
5007 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5009 scm_t_inum xx
= SCM_I_INUM (x
);
5010 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5012 scm_t_inum yy
= SCM_I_INUM (y
);
5015 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5016 scm_num_overflow (s_divide
);
5018 return scm_from_double ((double) xx
/ (double) yy
);
5021 else if (xx
% yy
!= 0)
5024 return scm_from_double ((double) xx
/ (double) yy
);
5025 else return scm_i_make_ratio (x
, y
);
5029 scm_t_inum z
= xx
/ yy
;
5030 if (SCM_FIXABLE (z
))
5031 return SCM_I_MAKINUM (z
);
5033 return scm_i_inum2big (z
);
5036 else if (SCM_BIGP (y
))
5039 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
5040 else return scm_i_make_ratio (x
, y
);
5042 else if (SCM_REALP (y
))
5044 double yy
= SCM_REAL_VALUE (y
);
5045 #ifndef ALLOW_DIVIDE_BY_ZERO
5047 scm_num_overflow (s_divide
);
5050 return scm_from_double ((double) xx
/ yy
);
5052 else if (SCM_COMPLEXP (y
))
5055 complex_div
: /* y _must_ be a complex number */
5057 double r
= SCM_COMPLEX_REAL (y
);
5058 double i
= SCM_COMPLEX_IMAG (y
);
5059 if (fabs(r
) <= fabs(i
))
5062 double d
= i
* (1.0 + t
* t
);
5063 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5068 double d
= r
* (1.0 + t
* t
);
5069 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5073 else if (SCM_FRACTIONP (y
))
5074 /* a / b/c = ac / b */
5075 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5076 SCM_FRACTION_NUMERATOR (y
));
5078 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5080 else if (SCM_BIGP (x
))
5082 if (SCM_I_INUMP (y
))
5084 scm_t_inum yy
= SCM_I_INUM (y
);
5087 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5088 scm_num_overflow (s_divide
);
5090 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5091 scm_remember_upto_here_1 (x
);
5092 return (sgn
== 0) ? scm_nan () : scm_inf ();
5099 /* FIXME: HMM, what are the relative performance issues here?
5100 We need to test. Is it faster on average to test
5101 divisible_p, then perform whichever operation, or is it
5102 faster to perform the integer div opportunistically and
5103 switch to real if there's a remainder? For now we take the
5104 middle ground: test, then if divisible, use the faster div
5107 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
5108 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5112 SCM result
= scm_i_mkbig ();
5113 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5114 scm_remember_upto_here_1 (x
);
5116 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5117 return scm_i_normbig (result
);
5122 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5123 else return scm_i_make_ratio (x
, y
);
5127 else if (SCM_BIGP (y
))
5129 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5132 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5133 scm_num_overflow (s_divide
);
5135 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5136 scm_remember_upto_here_1 (x
);
5137 return (sgn
== 0) ? scm_nan () : scm_inf ();
5145 /* It's easily possible for the ratio x/y to fit a double
5146 but one or both x and y be too big to fit a double,
5147 hence the use of mpq_get_d rather than converting and
5150 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5151 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5152 return scm_from_double (mpq_get_d (q
));
5156 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5160 SCM result
= scm_i_mkbig ();
5161 mpz_divexact (SCM_I_BIG_MPZ (result
),
5164 scm_remember_upto_here_2 (x
, y
);
5165 return scm_i_normbig (result
);
5168 return scm_i_make_ratio (x
, y
);
5172 else if (SCM_REALP (y
))
5174 double yy
= SCM_REAL_VALUE (y
);
5175 #ifndef ALLOW_DIVIDE_BY_ZERO
5177 scm_num_overflow (s_divide
);
5180 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5182 else if (SCM_COMPLEXP (y
))
5184 a
= scm_i_big2dbl (x
);
5187 else if (SCM_FRACTIONP (y
))
5188 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5189 SCM_FRACTION_NUMERATOR (y
));
5191 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5193 else if (SCM_REALP (x
))
5195 double rx
= SCM_REAL_VALUE (x
);
5196 if (SCM_I_INUMP (y
))
5198 scm_t_inum yy
= SCM_I_INUM (y
);
5199 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5201 scm_num_overflow (s_divide
);
5204 return scm_from_double (rx
/ (double) yy
);
5206 else if (SCM_BIGP (y
))
5208 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5209 scm_remember_upto_here_1 (y
);
5210 return scm_from_double (rx
/ dby
);
5212 else if (SCM_REALP (y
))
5214 double yy
= SCM_REAL_VALUE (y
);
5215 #ifndef ALLOW_DIVIDE_BY_ZERO
5217 scm_num_overflow (s_divide
);
5220 return scm_from_double (rx
/ yy
);
5222 else if (SCM_COMPLEXP (y
))
5227 else if (SCM_FRACTIONP (y
))
5228 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5230 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5232 else if (SCM_COMPLEXP (x
))
5234 double rx
= SCM_COMPLEX_REAL (x
);
5235 double ix
= SCM_COMPLEX_IMAG (x
);
5236 if (SCM_I_INUMP (y
))
5238 scm_t_inum yy
= SCM_I_INUM (y
);
5239 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5241 scm_num_overflow (s_divide
);
5246 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5249 else if (SCM_BIGP (y
))
5251 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5252 scm_remember_upto_here_1 (y
);
5253 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5255 else if (SCM_REALP (y
))
5257 double yy
= SCM_REAL_VALUE (y
);
5258 #ifndef ALLOW_DIVIDE_BY_ZERO
5260 scm_num_overflow (s_divide
);
5263 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5265 else if (SCM_COMPLEXP (y
))
5267 double ry
= SCM_COMPLEX_REAL (y
);
5268 double iy
= SCM_COMPLEX_IMAG (y
);
5269 if (fabs(ry
) <= fabs(iy
))
5272 double d
= iy
* (1.0 + t
* t
);
5273 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5278 double d
= ry
* (1.0 + t
* t
);
5279 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5282 else if (SCM_FRACTIONP (y
))
5284 double yy
= scm_i_fraction2double (y
);
5285 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5288 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5290 else if (SCM_FRACTIONP (x
))
5292 if (SCM_I_INUMP (y
))
5294 scm_t_inum yy
= SCM_I_INUM (y
);
5295 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5297 scm_num_overflow (s_divide
);
5300 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5301 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5303 else if (SCM_BIGP (y
))
5305 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5306 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5308 else if (SCM_REALP (y
))
5310 double yy
= SCM_REAL_VALUE (y
);
5311 #ifndef ALLOW_DIVIDE_BY_ZERO
5313 scm_num_overflow (s_divide
);
5316 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5318 else if (SCM_COMPLEXP (y
))
5320 a
= scm_i_fraction2double (x
);
5323 else if (SCM_FRACTIONP (y
))
5324 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5325 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5327 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5330 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5334 scm_divide (SCM x
, SCM y
)
5336 return do_divide (x
, y
, 0);
5339 static SCM
scm_divide2real (SCM x
, SCM y
)
5341 return do_divide (x
, y
, 1);
5347 scm_c_truncate (double x
)
5358 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5359 half-way case (ie. when x is an integer plus 0.5) going upwards.
5360 Then half-way cases are identified and adjusted down if the
5361 round-upwards didn't give the desired even integer.
5363 "plus_half == result" identifies a half-way case. If plus_half, which is
5364 x + 0.5, is an integer then x must be an integer plus 0.5.
5366 An odd "result" value is identified with result/2 != floor(result/2).
5367 This is done with plus_half, since that value is ready for use sooner in
5368 a pipelined cpu, and we're already requiring plus_half == result.
5370 Note however that we need to be careful when x is big and already an
5371 integer. In that case "x+0.5" may round to an adjacent integer, causing
5372 us to return such a value, incorrectly. For instance if the hardware is
5373 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5374 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5375 returned. Or if the hardware is in round-upwards mode, then other bigger
5376 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5377 representable value, 2^128+2^76 (or whatever), again incorrect.
5379 These bad roundings of x+0.5 are avoided by testing at the start whether
5380 x is already an integer. If it is then clearly that's the desired result
5381 already. And if it's not then the exponent must be small enough to allow
5382 an 0.5 to be represented, and hence added without a bad rounding. */
5385 scm_c_round (double x
)
5387 double plus_half
, result
;
5392 plus_half
= x
+ 0.5;
5393 result
= floor (plus_half
);
5394 /* Adjust so that the rounding is towards even. */
5395 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5400 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5402 "Round the number @var{x} towards zero.")
5403 #define FUNC_NAME s_scm_truncate_number
5405 if (scm_is_false (scm_negative_p (x
)))
5406 return scm_floor (x
);
5408 return scm_ceiling (x
);
5412 static SCM exactly_one_half
;
5414 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5416 "Round the number @var{x} towards the nearest integer. "
5417 "When it is exactly halfway between two integers, "
5418 "round towards the even one.")
5419 #define FUNC_NAME s_scm_round_number
5421 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5423 else if (SCM_REALP (x
))
5424 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5427 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5428 single quotient+remainder division then examining to see which way
5429 the rounding should go. */
5430 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5431 SCM result
= scm_floor (plus_half
);
5432 /* Adjust so that the rounding is towards even. */
5433 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5434 && scm_is_true (scm_odd_p (result
)))
5435 return scm_difference (result
, SCM_INUM1
);
5442 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5444 "Round the number @var{x} towards minus infinity.")
5445 #define FUNC_NAME s_scm_floor
5447 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5449 else if (SCM_REALP (x
))
5450 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5451 else if (SCM_FRACTIONP (x
))
5453 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5454 SCM_FRACTION_DENOMINATOR (x
));
5455 if (scm_is_false (scm_negative_p (x
)))
5457 /* For positive x, rounding towards zero is correct. */
5462 /* For negative x, we need to return q-1 unless x is an
5463 integer. But fractions are never integer, per our
5465 return scm_difference (q
, SCM_INUM1
);
5469 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5473 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5475 "Round the number @var{x} towards infinity.")
5476 #define FUNC_NAME s_scm_ceiling
5478 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5480 else if (SCM_REALP (x
))
5481 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5482 else if (SCM_FRACTIONP (x
))
5484 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5485 SCM_FRACTION_DENOMINATOR (x
));
5486 if (scm_is_false (scm_positive_p (x
)))
5488 /* For negative x, rounding towards zero is correct. */
5493 /* For positive x, we need to return q+1 unless x is an
5494 integer. But fractions are never integer, per our
5496 return scm_sum (q
, SCM_INUM1
);
5500 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5504 /* sin/cos/tan/asin/acos/atan
5505 sinh/cosh/tanh/asinh/acosh/atanh
5506 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5507 Written by Jerry D. Hedden, (C) FSF.
5508 See the file `COPYING' for terms applying to this program. */
5510 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5512 "Return @var{x} raised to the power of @var{y}.")
5513 #define FUNC_NAME s_scm_expt
5515 if (scm_is_integer (y
))
5517 if (scm_is_true (scm_exact_p (y
)))
5518 return scm_integer_expt (x
, y
);
5521 /* Here we handle the case where the exponent is an inexact
5522 integer. We make the exponent exact in order to use
5523 scm_integer_expt, and thus avoid the spurious imaginary
5524 parts that may result from round-off errors in the general
5525 e^(y log x) method below (for example when squaring a large
5526 negative number). In this case, we must return an inexact
5527 result for correctness. We also make the base inexact so
5528 that scm_integer_expt will use fast inexact arithmetic
5529 internally. Note that making the base inexact is not
5530 sufficient to guarantee an inexact result, because
5531 scm_integer_expt will return an exact 1 when the exponent
5532 is 0, even if the base is inexact. */
5533 return scm_exact_to_inexact
5534 (scm_integer_expt (scm_exact_to_inexact (x
),
5535 scm_inexact_to_exact (y
)));
5538 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5540 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5543 return scm_exp (scm_product (scm_log (x
), y
));
5547 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5549 "Compute the sine of @var{z}.")
5550 #define FUNC_NAME s_scm_sin
5552 if (scm_is_real (z
))
5553 return scm_from_double (sin (scm_to_double (z
)));
5554 else if (SCM_COMPLEXP (z
))
5556 x
= SCM_COMPLEX_REAL (z
);
5557 y
= SCM_COMPLEX_IMAG (z
);
5558 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5559 cos (x
) * sinh (y
));
5562 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5566 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5568 "Compute the cosine of @var{z}.")
5569 #define FUNC_NAME s_scm_cos
5571 if (scm_is_real (z
))
5572 return scm_from_double (cos (scm_to_double (z
)));
5573 else if (SCM_COMPLEXP (z
))
5575 x
= SCM_COMPLEX_REAL (z
);
5576 y
= SCM_COMPLEX_IMAG (z
);
5577 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5578 -sin (x
) * sinh (y
));
5581 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5585 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5587 "Compute the tangent of @var{z}.")
5588 #define FUNC_NAME s_scm_tan
5590 if (scm_is_real (z
))
5591 return scm_from_double (tan (scm_to_double (z
)));
5592 else if (SCM_COMPLEXP (z
))
5594 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5595 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5596 w
= cos (x
) + cosh (y
);
5597 #ifndef ALLOW_DIVIDE_BY_ZERO
5599 scm_num_overflow (s_scm_tan
);
5601 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5604 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5608 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5610 "Compute the hyperbolic sine of @var{z}.")
5611 #define FUNC_NAME s_scm_sinh
5613 if (scm_is_real (z
))
5614 return scm_from_double (sinh (scm_to_double (z
)));
5615 else if (SCM_COMPLEXP (z
))
5617 x
= SCM_COMPLEX_REAL (z
);
5618 y
= SCM_COMPLEX_IMAG (z
);
5619 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5620 cosh (x
) * sin (y
));
5623 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5627 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5629 "Compute the hyperbolic cosine of @var{z}.")
5630 #define FUNC_NAME s_scm_cosh
5632 if (scm_is_real (z
))
5633 return scm_from_double (cosh (scm_to_double (z
)));
5634 else if (SCM_COMPLEXP (z
))
5636 x
= SCM_COMPLEX_REAL (z
);
5637 y
= SCM_COMPLEX_IMAG (z
);
5638 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5639 sinh (x
) * sin (y
));
5642 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5646 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5648 "Compute the hyperbolic tangent of @var{z}.")
5649 #define FUNC_NAME s_scm_tanh
5651 if (scm_is_real (z
))
5652 return scm_from_double (tanh (scm_to_double (z
)));
5653 else if (SCM_COMPLEXP (z
))
5655 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5656 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5657 w
= cosh (x
) + cos (y
);
5658 #ifndef ALLOW_DIVIDE_BY_ZERO
5660 scm_num_overflow (s_scm_tanh
);
5662 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5665 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5669 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5671 "Compute the arc sine of @var{z}.")
5672 #define FUNC_NAME s_scm_asin
5674 if (scm_is_real (z
))
5676 double w
= scm_to_double (z
);
5677 if (w
>= -1.0 && w
<= 1.0)
5678 return scm_from_double (asin (w
));
5680 return scm_product (scm_c_make_rectangular (0, -1),
5681 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5683 else if (SCM_COMPLEXP (z
))
5685 x
= SCM_COMPLEX_REAL (z
);
5686 y
= SCM_COMPLEX_IMAG (z
);
5687 return scm_product (scm_c_make_rectangular (0, -1),
5688 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5691 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5695 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5697 "Compute the arc cosine of @var{z}.")
5698 #define FUNC_NAME s_scm_acos
5700 if (scm_is_real (z
))
5702 double w
= scm_to_double (z
);
5703 if (w
>= -1.0 && w
<= 1.0)
5704 return scm_from_double (acos (w
));
5706 return scm_sum (scm_from_double (acos (0.0)),
5707 scm_product (scm_c_make_rectangular (0, 1),
5708 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5710 else if (SCM_COMPLEXP (z
))
5712 x
= SCM_COMPLEX_REAL (z
);
5713 y
= SCM_COMPLEX_IMAG (z
);
5714 return scm_sum (scm_from_double (acos (0.0)),
5715 scm_product (scm_c_make_rectangular (0, 1),
5716 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5719 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5723 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5725 "With one argument, compute the arc tangent of @var{z}.\n"
5726 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5727 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5728 #define FUNC_NAME s_scm_atan
5732 if (scm_is_real (z
))
5733 return scm_from_double (atan (scm_to_double (z
)));
5734 else if (SCM_COMPLEXP (z
))
5737 v
= SCM_COMPLEX_REAL (z
);
5738 w
= SCM_COMPLEX_IMAG (z
);
5739 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5740 scm_c_make_rectangular (v
, w
+ 1.0))),
5741 scm_c_make_rectangular (0, 2));
5744 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5746 else if (scm_is_real (z
))
5748 if (scm_is_real (y
))
5749 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5751 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5754 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5758 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5760 "Compute the inverse hyperbolic sine of @var{z}.")
5761 #define FUNC_NAME s_scm_sys_asinh
5763 if (scm_is_real (z
))
5764 return scm_from_double (asinh (scm_to_double (z
)));
5765 else if (scm_is_number (z
))
5766 return scm_log (scm_sum (z
,
5767 scm_sqrt (scm_sum (scm_product (z
, z
),
5770 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5774 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5776 "Compute the inverse hyperbolic cosine of @var{z}.")
5777 #define FUNC_NAME s_scm_sys_acosh
5779 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5780 return scm_from_double (acosh (scm_to_double (z
)));
5781 else if (scm_is_number (z
))
5782 return scm_log (scm_sum (z
,
5783 scm_sqrt (scm_difference (scm_product (z
, z
),
5786 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5790 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5792 "Compute the inverse hyperbolic tangent of @var{z}.")
5793 #define FUNC_NAME s_scm_sys_atanh
5795 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5796 return scm_from_double (atanh (scm_to_double (z
)));
5797 else if (scm_is_number (z
))
5798 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
5799 scm_difference (SCM_INUM1
, z
))),
5802 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5807 scm_c_make_rectangular (double re
, double im
)
5810 return scm_from_double (re
);
5815 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5817 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
5818 SCM_COMPLEX_REAL (z
) = re
;
5819 SCM_COMPLEX_IMAG (z
) = im
;
5824 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5825 (SCM real_part
, SCM imaginary_part
),
5826 "Return a complex number constructed of the given @var{real-part} "
5827 "and @var{imaginary-part} parts.")
5828 #define FUNC_NAME s_scm_make_rectangular
5830 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5831 SCM_ARG1
, FUNC_NAME
, "real");
5832 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5833 SCM_ARG2
, FUNC_NAME
, "real");
5834 return scm_c_make_rectangular (scm_to_double (real_part
),
5835 scm_to_double (imaginary_part
));
5840 scm_c_make_polar (double mag
, double ang
)
5844 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5845 use it on Glibc-based systems that have it (it's a GNU extension). See
5846 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5848 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5849 sincos (ang
, &s
, &c
);
5854 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5857 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5859 "Return the complex number @var{x} * e^(i * @var{y}).")
5860 #define FUNC_NAME s_scm_make_polar
5862 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5863 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5864 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5869 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5870 /* "Return the real part of the number @var{z}."
5873 scm_real_part (SCM z
)
5875 if (SCM_I_INUMP (z
))
5877 else if (SCM_BIGP (z
))
5879 else if (SCM_REALP (z
))
5881 else if (SCM_COMPLEXP (z
))
5882 return scm_from_double (SCM_COMPLEX_REAL (z
));
5883 else if (SCM_FRACTIONP (z
))
5886 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5890 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5891 /* "Return the imaginary part of the number @var{z}."
5894 scm_imag_part (SCM z
)
5896 if (SCM_I_INUMP (z
))
5898 else if (SCM_BIGP (z
))
5900 else if (SCM_REALP (z
))
5902 else if (SCM_COMPLEXP (z
))
5903 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5904 else if (SCM_FRACTIONP (z
))
5907 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5910 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5911 /* "Return the numerator of the number @var{z}."
5914 scm_numerator (SCM z
)
5916 if (SCM_I_INUMP (z
))
5918 else if (SCM_BIGP (z
))
5920 else if (SCM_FRACTIONP (z
))
5921 return SCM_FRACTION_NUMERATOR (z
);
5922 else if (SCM_REALP (z
))
5923 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5925 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5929 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5930 /* "Return the denominator of the number @var{z}."
5933 scm_denominator (SCM z
)
5935 if (SCM_I_INUMP (z
))
5937 else if (SCM_BIGP (z
))
5939 else if (SCM_FRACTIONP (z
))
5940 return SCM_FRACTION_DENOMINATOR (z
);
5941 else if (SCM_REALP (z
))
5942 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5944 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5947 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5948 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5949 * "@code{abs} for real arguments, but also allows complex numbers."
5952 scm_magnitude (SCM z
)
5954 if (SCM_I_INUMP (z
))
5956 scm_t_inum zz
= SCM_I_INUM (z
);
5959 else if (SCM_POSFIXABLE (-zz
))
5960 return SCM_I_MAKINUM (-zz
);
5962 return scm_i_inum2big (-zz
);
5964 else if (SCM_BIGP (z
))
5966 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5967 scm_remember_upto_here_1 (z
);
5969 return scm_i_clonebig (z
, 0);
5973 else if (SCM_REALP (z
))
5974 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5975 else if (SCM_COMPLEXP (z
))
5976 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5977 else if (SCM_FRACTIONP (z
))
5979 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5981 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5982 SCM_FRACTION_DENOMINATOR (z
));
5985 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5989 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5990 /* "Return the angle of the complex number @var{z}."
5995 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5996 flo0 to save allocating a new flonum with scm_from_double each time.
5997 But if atan2 follows the floating point rounding mode, then the value
5998 is not a constant. Maybe it'd be close enough though. */
5999 if (SCM_I_INUMP (z
))
6001 if (SCM_I_INUM (z
) >= 0)
6004 return scm_from_double (atan2 (0.0, -1.0));
6006 else if (SCM_BIGP (z
))
6008 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
6009 scm_remember_upto_here_1 (z
);
6011 return scm_from_double (atan2 (0.0, -1.0));
6015 else if (SCM_REALP (z
))
6017 if (SCM_REAL_VALUE (z
) >= 0)
6020 return scm_from_double (atan2 (0.0, -1.0));
6022 else if (SCM_COMPLEXP (z
))
6023 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
6024 else if (SCM_FRACTIONP (z
))
6026 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
6028 else return scm_from_double (atan2 (0.0, -1.0));
6031 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
6035 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
6036 /* Convert the number @var{x} to its inexact representation.\n"
6039 scm_exact_to_inexact (SCM z
)
6041 if (SCM_I_INUMP (z
))
6042 return scm_from_double ((double) SCM_I_INUM (z
));
6043 else if (SCM_BIGP (z
))
6044 return scm_from_double (scm_i_big2dbl (z
));
6045 else if (SCM_FRACTIONP (z
))
6046 return scm_from_double (scm_i_fraction2double (z
));
6047 else if (SCM_INEXACTP (z
))
6050 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
6054 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
6056 "Return an exact number that is numerically closest to @var{z}.")
6057 #define FUNC_NAME s_scm_inexact_to_exact
6059 if (SCM_I_INUMP (z
))
6061 else if (SCM_BIGP (z
))
6063 else if (SCM_REALP (z
))
6065 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
6066 SCM_OUT_OF_RANGE (1, z
);
6073 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6074 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6075 scm_i_mpz2num (mpq_denref (frac
)));
6077 /* When scm_i_make_ratio throws, we leak the memory allocated
6084 else if (SCM_FRACTIONP (z
))
6087 SCM_WRONG_TYPE_ARG (1, z
);
6091 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6093 "Returns the @emph{simplest} rational number differing\n"
6094 "from @var{x} by no more than @var{eps}.\n"
6096 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6097 "exact result when both its arguments are exact. Thus, you might need\n"
6098 "to use @code{inexact->exact} on the arguments.\n"
6101 "(rationalize (inexact->exact 1.2) 1/100)\n"
6104 #define FUNC_NAME s_scm_rationalize
6106 if (SCM_I_INUMP (x
))
6108 else if (SCM_BIGP (x
))
6110 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6112 /* Use continued fractions to find closest ratio. All
6113 arithmetic is done with exact numbers.
6116 SCM ex
= scm_inexact_to_exact (x
);
6117 SCM int_part
= scm_floor (ex
);
6119 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
6120 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
6124 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6127 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6128 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6130 /* We stop after a million iterations just to be absolutely sure
6131 that we don't go into an infinite loop. The process normally
6132 converges after less than a dozen iterations.
6135 eps
= scm_abs (eps
);
6136 while (++i
< 1000000)
6138 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6139 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6140 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6142 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6143 eps
))) /* abs(x-a/b) <= eps */
6145 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6146 if (scm_is_false (scm_exact_p (x
))
6147 || scm_is_false (scm_exact_p (eps
)))
6148 return scm_exact_to_inexact (res
);
6152 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6154 tt
= scm_floor (rx
); /* tt = floor (rx) */
6160 scm_num_overflow (s_scm_rationalize
);
6163 SCM_WRONG_TYPE_ARG (1, x
);
6167 /* conversion functions */
6170 scm_is_integer (SCM val
)
6172 return scm_is_true (scm_integer_p (val
));
6176 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6178 if (SCM_I_INUMP (val
))
6180 scm_t_signed_bits n
= SCM_I_INUM (val
);
6181 return n
>= min
&& n
<= max
;
6183 else if (SCM_BIGP (val
))
6185 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6187 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6189 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6191 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6192 return n
>= min
&& n
<= max
;
6202 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6203 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6206 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6207 SCM_I_BIG_MPZ (val
));
6209 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6221 return n
>= min
&& n
<= max
;
6229 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6231 if (SCM_I_INUMP (val
))
6233 scm_t_signed_bits n
= SCM_I_INUM (val
);
6234 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6236 else if (SCM_BIGP (val
))
6238 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6240 else if (max
<= ULONG_MAX
)
6242 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6244 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6245 return n
>= min
&& n
<= max
;
6255 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6258 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6259 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6262 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6263 SCM_I_BIG_MPZ (val
));
6265 return n
>= min
&& n
<= max
;
6273 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6275 scm_error (scm_out_of_range_key
,
6277 "Value out of range ~S to ~S: ~S",
6278 scm_list_3 (min
, max
, bad_val
),
6279 scm_list_1 (bad_val
));
6282 #define TYPE scm_t_intmax
6283 #define TYPE_MIN min
6284 #define TYPE_MAX max
6285 #define SIZEOF_TYPE 0
6286 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6287 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6288 #include "libguile/conv-integer.i.c"
6290 #define TYPE scm_t_uintmax
6291 #define TYPE_MIN min
6292 #define TYPE_MAX max
6293 #define SIZEOF_TYPE 0
6294 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6295 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6296 #include "libguile/conv-uinteger.i.c"
6298 #define TYPE scm_t_int8
6299 #define TYPE_MIN SCM_T_INT8_MIN
6300 #define TYPE_MAX SCM_T_INT8_MAX
6301 #define SIZEOF_TYPE 1
6302 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6303 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6304 #include "libguile/conv-integer.i.c"
6306 #define TYPE scm_t_uint8
6308 #define TYPE_MAX SCM_T_UINT8_MAX
6309 #define SIZEOF_TYPE 1
6310 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6311 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6312 #include "libguile/conv-uinteger.i.c"
6314 #define TYPE scm_t_int16
6315 #define TYPE_MIN SCM_T_INT16_MIN
6316 #define TYPE_MAX SCM_T_INT16_MAX
6317 #define SIZEOF_TYPE 2
6318 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6319 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6320 #include "libguile/conv-integer.i.c"
6322 #define TYPE scm_t_uint16
6324 #define TYPE_MAX SCM_T_UINT16_MAX
6325 #define SIZEOF_TYPE 2
6326 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6327 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6328 #include "libguile/conv-uinteger.i.c"
6330 #define TYPE scm_t_int32
6331 #define TYPE_MIN SCM_T_INT32_MIN
6332 #define TYPE_MAX SCM_T_INT32_MAX
6333 #define SIZEOF_TYPE 4
6334 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6335 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6336 #include "libguile/conv-integer.i.c"
6338 #define TYPE scm_t_uint32
6340 #define TYPE_MAX SCM_T_UINT32_MAX
6341 #define SIZEOF_TYPE 4
6342 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6343 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6344 #include "libguile/conv-uinteger.i.c"
6346 #define TYPE scm_t_wchar
6347 #define TYPE_MIN (scm_t_int32)-1
6348 #define TYPE_MAX (scm_t_int32)0x10ffff
6349 #define SIZEOF_TYPE 4
6350 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6351 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6352 #include "libguile/conv-integer.i.c"
6354 #define TYPE scm_t_int64
6355 #define TYPE_MIN SCM_T_INT64_MIN
6356 #define TYPE_MAX SCM_T_INT64_MAX
6357 #define SIZEOF_TYPE 8
6358 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6359 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6360 #include "libguile/conv-integer.i.c"
6362 #define TYPE scm_t_uint64
6364 #define TYPE_MAX SCM_T_UINT64_MAX
6365 #define SIZEOF_TYPE 8
6366 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6367 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6368 #include "libguile/conv-uinteger.i.c"
6371 scm_to_mpz (SCM val
, mpz_t rop
)
6373 if (SCM_I_INUMP (val
))
6374 mpz_set_si (rop
, SCM_I_INUM (val
));
6375 else if (SCM_BIGP (val
))
6376 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6378 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6382 scm_from_mpz (mpz_t val
)
6384 return scm_i_mpz2num (val
);
6388 scm_is_real (SCM val
)
6390 return scm_is_true (scm_real_p (val
));
6394 scm_is_rational (SCM val
)
6396 return scm_is_true (scm_rational_p (val
));
6400 scm_to_double (SCM val
)
6402 if (SCM_I_INUMP (val
))
6403 return SCM_I_INUM (val
);
6404 else if (SCM_BIGP (val
))
6405 return scm_i_big2dbl (val
);
6406 else if (SCM_FRACTIONP (val
))
6407 return scm_i_fraction2double (val
);
6408 else if (SCM_REALP (val
))
6409 return SCM_REAL_VALUE (val
);
6411 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6415 scm_from_double (double val
)
6419 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
6421 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
6422 SCM_REAL_VALUE (z
) = val
;
6427 #if SCM_ENABLE_DEPRECATED == 1
6430 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
6432 scm_c_issue_deprecation_warning
6433 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6437 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6441 scm_out_of_range (NULL
, num
);
6444 return scm_to_double (num
);
6448 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
6450 scm_c_issue_deprecation_warning
6451 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6455 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6459 scm_out_of_range (NULL
, num
);
6462 return scm_to_double (num
);
6468 scm_is_complex (SCM val
)
6470 return scm_is_true (scm_complex_p (val
));
6474 scm_c_real_part (SCM z
)
6476 if (SCM_COMPLEXP (z
))
6477 return SCM_COMPLEX_REAL (z
);
6480 /* Use the scm_real_part to get proper error checking and
6483 return scm_to_double (scm_real_part (z
));
6488 scm_c_imag_part (SCM z
)
6490 if (SCM_COMPLEXP (z
))
6491 return SCM_COMPLEX_IMAG (z
);
6494 /* Use the scm_imag_part to get proper error checking and
6495 dispatching. The result will almost always be 0.0, but not
6498 return scm_to_double (scm_imag_part (z
));
6503 scm_c_magnitude (SCM z
)
6505 return scm_to_double (scm_magnitude (z
));
6511 return scm_to_double (scm_angle (z
));
6515 scm_is_number (SCM z
)
6517 return scm_is_true (scm_number_p (z
));
6521 /* In the following functions we dispatch to the real-arg funcs like log()
6522 when we know the arg is real, instead of just handing everything to
6523 clog() for instance. This is in case clog() doesn't optimize for a
6524 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6525 well use it to go straight to the applicable C func. */
6527 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6529 "Return the natural logarithm of @var{z}.")
6530 #define FUNC_NAME s_scm_log
6532 if (SCM_COMPLEXP (z
))
6534 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6535 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6537 double re
= SCM_COMPLEX_REAL (z
);
6538 double im
= SCM_COMPLEX_IMAG (z
);
6539 return scm_c_make_rectangular (log (hypot (re
, im
)),
6545 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6546 although the value itself overflows. */
6547 double re
= scm_to_double (z
);
6548 double l
= log (fabs (re
));
6550 return scm_from_double (l
);
6552 return scm_c_make_rectangular (l
, M_PI
);
6558 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6560 "Return the base 10 logarithm of @var{z}.")
6561 #define FUNC_NAME s_scm_log10
6563 if (SCM_COMPLEXP (z
))
6565 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6566 clog() and a multiply by M_LOG10E, rather than the fallback
6567 log10+hypot+atan2.) */
6568 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6569 && defined SCM_COMPLEX_VALUE
6570 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6572 double re
= SCM_COMPLEX_REAL (z
);
6573 double im
= SCM_COMPLEX_IMAG (z
);
6574 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6575 M_LOG10E
* atan2 (im
, re
));
6580 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6581 although the value itself overflows. */
6582 double re
= scm_to_double (z
);
6583 double l
= log10 (fabs (re
));
6585 return scm_from_double (l
);
6587 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6593 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6595 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6596 "base of natural logarithms (2.71828@dots{}).")
6597 #define FUNC_NAME s_scm_exp
6599 if (SCM_COMPLEXP (z
))
6601 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6602 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6604 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6605 SCM_COMPLEX_IMAG (z
));
6610 /* When z is a negative bignum the conversion to double overflows,
6611 giving -infinity, but that's ok, the exp is still 0.0. */
6612 return scm_from_double (exp (scm_to_double (z
)));
6618 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6620 "Return the square root of @var{z}. Of the two possible roots\n"
6621 "(positive and negative), the one with the a positive real part\n"
6622 "is returned, or if that's zero then a positive imaginary part.\n"
6626 "(sqrt 9.0) @result{} 3.0\n"
6627 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6628 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6629 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6631 #define FUNC_NAME s_scm_sqrt
6633 if (SCM_COMPLEXP (x
))
6635 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6636 && defined SCM_COMPLEX_VALUE
6637 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6639 double re
= SCM_COMPLEX_REAL (x
);
6640 double im
= SCM_COMPLEX_IMAG (x
);
6641 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6642 0.5 * atan2 (im
, re
));
6647 double xx
= scm_to_double (x
);
6649 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6651 return scm_from_double (sqrt (xx
));
6663 mpz_init_set_si (z_negative_one
, -1);
6665 /* It may be possible to tune the performance of some algorithms by using
6666 * the following constants to avoid the creation of bignums. Please, before
6667 * using these values, remember the two rules of program optimization:
6668 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6669 scm_c_define ("most-positive-fixnum",
6670 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6671 scm_c_define ("most-negative-fixnum",
6672 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6674 scm_add_feature ("complex");
6675 scm_add_feature ("inexact");
6676 flo0
= scm_from_double (0.0);
6678 /* determine floating point precision */
6679 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6681 init_dblprec(&scm_dblprec
[i
-2],i
);
6682 init_fx_radix(fx_per_radix
[i
-2],i
);
6685 /* hard code precision for base 10 if the preprocessor tells us to... */
6686 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6689 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
6690 #include "libguile/numbers.x"