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
506 if (SCM_INEXACTP (x
))
508 else if (SCM_NUMBERP (x
))
511 SCM_WRONG_TYPE_ARG (1, x
);
516 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
518 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
520 #define FUNC_NAME s_scm_inexact_p
522 if (SCM_INEXACTP (x
))
524 else if (SCM_NUMBERP (x
))
527 SCM_WRONG_TYPE_ARG (1, x
);
532 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
534 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
536 #define FUNC_NAME s_scm_odd_p
540 scm_t_inum val
= SCM_I_INUM (n
);
541 return scm_from_bool ((val
& 1L) != 0);
543 else if (SCM_BIGP (n
))
545 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
546 scm_remember_upto_here_1 (n
);
547 return scm_from_bool (odd_p
);
549 else if (scm_is_true (scm_inf_p (n
)))
550 SCM_WRONG_TYPE_ARG (1, n
);
551 else if (SCM_REALP (n
))
553 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
559 SCM_WRONG_TYPE_ARG (1, n
);
562 SCM_WRONG_TYPE_ARG (1, n
);
567 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
569 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
571 #define FUNC_NAME s_scm_even_p
575 scm_t_inum val
= SCM_I_INUM (n
);
576 return scm_from_bool ((val
& 1L) == 0);
578 else if (SCM_BIGP (n
))
580 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
581 scm_remember_upto_here_1 (n
);
582 return scm_from_bool (even_p
);
584 else if (scm_is_true (scm_inf_p (n
)))
585 SCM_WRONG_TYPE_ARG (1, n
);
586 else if (SCM_REALP (n
))
588 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
594 SCM_WRONG_TYPE_ARG (1, n
);
597 SCM_WRONG_TYPE_ARG (1, n
);
601 SCM_DEFINE (scm_finite_p
, "finite?", 1, 0, 0,
603 "Return @code{#t} if @var{x} is neither infinite\n"
604 "nor a NaN, @code{#f} otherwise.")
605 #define FUNC_NAME s_scm_finite_p
608 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
609 else if (SCM_COMPLEXP (x
))
610 return scm_from_bool (DOUBLE_IS_FINITE (SCM_COMPLEX_REAL (x
))
611 && DOUBLE_IS_FINITE (SCM_COMPLEX_IMAG (x
)));
612 else if (SCM_NUMBERP (x
))
615 SCM_WRONG_TYPE_ARG (1, x
);
619 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
621 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
622 "or @samp{-inf.0}, @code{#f} otherwise.")
623 #define FUNC_NAME s_scm_inf_p
626 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
627 else if (SCM_COMPLEXP (x
))
628 return scm_from_bool (isinf (SCM_COMPLEX_REAL (x
))
629 || isinf (SCM_COMPLEX_IMAG (x
)));
635 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
637 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
639 #define FUNC_NAME s_scm_nan_p
642 return scm_from_bool (isnan (SCM_REAL_VALUE (n
)));
643 else if (SCM_COMPLEXP (n
))
644 return scm_from_bool (isnan (SCM_COMPLEX_REAL (n
))
645 || isnan (SCM_COMPLEX_IMAG (n
)));
651 /* Guile's idea of infinity. */
652 static double guile_Inf
;
654 /* Guile's idea of not a number. */
655 static double guile_NaN
;
658 guile_ieee_init (void)
660 /* Some version of gcc on some old version of Linux used to crash when
661 trying to make Inf and NaN. */
664 /* C99 INFINITY, when available.
665 FIXME: The standard allows for INFINITY to be something that overflows
666 at compile time. We ought to have a configure test to check for that
667 before trying to use it. (But in practice we believe this is not a
668 problem on any system guile is likely to target.) */
669 guile_Inf
= INFINITY
;
670 #elif defined HAVE_DINFINITY
672 extern unsigned int DINFINITY
[2];
673 guile_Inf
= (*((double *) (DINFINITY
)));
680 if (guile_Inf
== tmp
)
687 /* C99 NAN, when available */
689 #elif defined HAVE_DQNAN
692 extern unsigned int DQNAN
[2];
693 guile_NaN
= (*((double *)(DQNAN
)));
696 guile_NaN
= guile_Inf
/ guile_Inf
;
700 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
703 #define FUNC_NAME s_scm_inf
705 static int initialized
= 0;
711 return scm_from_double (guile_Inf
);
715 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
718 #define FUNC_NAME s_scm_nan
720 static int initialized
= 0;
726 return scm_from_double (guile_NaN
);
731 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
733 "Return the absolute value of @var{x}.")
738 scm_t_inum xx
= SCM_I_INUM (x
);
741 else if (SCM_POSFIXABLE (-xx
))
742 return SCM_I_MAKINUM (-xx
);
744 return scm_i_inum2big (-xx
);
746 else if (SCM_BIGP (x
))
748 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
750 return scm_i_clonebig (x
, 0);
754 else if (SCM_REALP (x
))
756 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
757 double xx
= SCM_REAL_VALUE (x
);
759 return scm_from_double (-xx
);
763 else if (SCM_FRACTIONP (x
))
765 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
767 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
768 SCM_FRACTION_DENOMINATOR (x
));
771 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
776 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
777 /* "Return the quotient of the numbers @var{x} and @var{y}."
780 scm_quotient (SCM x
, SCM y
)
784 scm_t_inum xx
= SCM_I_INUM (x
);
787 scm_t_inum yy
= SCM_I_INUM (y
);
789 scm_num_overflow (s_quotient
);
792 scm_t_inum z
= xx
/ yy
;
794 return SCM_I_MAKINUM (z
);
796 return scm_i_inum2big (z
);
799 else if (SCM_BIGP (y
))
801 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
802 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
803 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
805 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
806 scm_remember_upto_here_1 (y
);
807 return SCM_I_MAKINUM (-1);
813 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
815 else if (SCM_BIGP (x
))
819 scm_t_inum yy
= SCM_I_INUM (y
);
821 scm_num_overflow (s_quotient
);
826 SCM result
= scm_i_mkbig ();
829 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
832 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
835 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
836 scm_remember_upto_here_1 (x
);
837 return scm_i_normbig (result
);
840 else if (SCM_BIGP (y
))
842 SCM result
= scm_i_mkbig ();
843 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
846 scm_remember_upto_here_2 (x
, y
);
847 return scm_i_normbig (result
);
850 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
853 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
856 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
857 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
859 * "(remainder 13 4) @result{} 1\n"
860 * "(remainder -13 4) @result{} -1\n"
864 scm_remainder (SCM x
, SCM y
)
870 scm_t_inum yy
= SCM_I_INUM (y
);
872 scm_num_overflow (s_remainder
);
875 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
876 return SCM_I_MAKINUM (z
);
879 else if (SCM_BIGP (y
))
881 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
882 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
883 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
885 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
886 scm_remember_upto_here_1 (y
);
893 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
895 else if (SCM_BIGP (x
))
899 scm_t_inum yy
= SCM_I_INUM (y
);
901 scm_num_overflow (s_remainder
);
904 SCM result
= scm_i_mkbig ();
907 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
908 scm_remember_upto_here_1 (x
);
909 return scm_i_normbig (result
);
912 else if (SCM_BIGP (y
))
914 SCM result
= scm_i_mkbig ();
915 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
918 scm_remember_upto_here_2 (x
, y
);
919 return scm_i_normbig (result
);
922 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
925 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
929 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
930 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
932 * "(modulo 13 4) @result{} 1\n"
933 * "(modulo -13 4) @result{} 3\n"
937 scm_modulo (SCM x
, SCM y
)
941 scm_t_inum xx
= SCM_I_INUM (x
);
944 scm_t_inum yy
= SCM_I_INUM (y
);
946 scm_num_overflow (s_modulo
);
949 /* C99 specifies that "%" is the remainder corresponding to a
950 quotient rounded towards zero, and that's also traditional
951 for machine division, so z here should be well defined. */
952 scm_t_inum z
= xx
% yy
;
969 return SCM_I_MAKINUM (result
);
972 else if (SCM_BIGP (y
))
974 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
981 SCM pos_y
= scm_i_clonebig (y
, 0);
982 /* do this after the last scm_op */
983 mpz_init_set_si (z_x
, xx
);
984 result
= pos_y
; /* re-use this bignum */
985 mpz_mod (SCM_I_BIG_MPZ (result
),
987 SCM_I_BIG_MPZ (pos_y
));
988 scm_remember_upto_here_1 (pos_y
);
992 result
= scm_i_mkbig ();
993 /* do this after the last scm_op */
994 mpz_init_set_si (z_x
, xx
);
995 mpz_mod (SCM_I_BIG_MPZ (result
),
998 scm_remember_upto_here_1 (y
);
1001 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1002 mpz_add (SCM_I_BIG_MPZ (result
),
1004 SCM_I_BIG_MPZ (result
));
1005 scm_remember_upto_here_1 (y
);
1006 /* and do this before the next one */
1008 return scm_i_normbig (result
);
1012 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1014 else if (SCM_BIGP (x
))
1016 if (SCM_I_INUMP (y
))
1018 scm_t_inum yy
= SCM_I_INUM (y
);
1020 scm_num_overflow (s_modulo
);
1023 SCM result
= scm_i_mkbig ();
1024 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1026 (yy
< 0) ? - yy
: yy
);
1027 scm_remember_upto_here_1 (x
);
1028 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1029 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1030 SCM_I_BIG_MPZ (result
),
1032 return scm_i_normbig (result
);
1035 else if (SCM_BIGP (y
))
1038 SCM result
= scm_i_mkbig ();
1039 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1040 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1041 mpz_mod (SCM_I_BIG_MPZ (result
),
1043 SCM_I_BIG_MPZ (pos_y
));
1045 scm_remember_upto_here_1 (x
);
1046 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1047 mpz_add (SCM_I_BIG_MPZ (result
),
1049 SCM_I_BIG_MPZ (result
));
1050 scm_remember_upto_here_2 (y
, pos_y
);
1051 return scm_i_normbig (result
);
1055 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1058 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1061 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1062 (SCM x
, SCM y
, SCM rest
),
1063 "Return the greatest common divisor of all parameter values.\n"
1064 "If called without arguments, 0 is returned.")
1065 #define FUNC_NAME s_scm_i_gcd
1067 while (!scm_is_null (rest
))
1068 { x
= scm_gcd (x
, y
);
1070 rest
= scm_cdr (rest
);
1072 return scm_gcd (x
, y
);
1076 #define s_gcd s_scm_i_gcd
1077 #define g_gcd g_scm_i_gcd
1080 scm_gcd (SCM x
, SCM y
)
1083 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1085 if (SCM_I_INUMP (x
))
1087 if (SCM_I_INUMP (y
))
1089 scm_t_inum xx
= SCM_I_INUM (x
);
1090 scm_t_inum yy
= SCM_I_INUM (y
);
1091 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
1092 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
1102 /* Determine a common factor 2^k */
1103 while (!(1 & (u
| v
)))
1109 /* Now, any factor 2^n can be eliminated */
1129 return (SCM_POSFIXABLE (result
)
1130 ? SCM_I_MAKINUM (result
)
1131 : scm_i_inum2big (result
));
1133 else if (SCM_BIGP (y
))
1139 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1141 else if (SCM_BIGP (x
))
1143 if (SCM_I_INUMP (y
))
1148 yy
= SCM_I_INUM (y
);
1153 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1154 scm_remember_upto_here_1 (x
);
1155 return (SCM_POSFIXABLE (result
)
1156 ? SCM_I_MAKINUM (result
)
1157 : scm_from_unsigned_integer (result
));
1159 else if (SCM_BIGP (y
))
1161 SCM result
= scm_i_mkbig ();
1162 mpz_gcd (SCM_I_BIG_MPZ (result
),
1165 scm_remember_upto_here_2 (x
, y
);
1166 return scm_i_normbig (result
);
1169 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1172 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1175 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1176 (SCM x
, SCM y
, SCM rest
),
1177 "Return the least common multiple of the arguments.\n"
1178 "If called without arguments, 1 is returned.")
1179 #define FUNC_NAME s_scm_i_lcm
1181 while (!scm_is_null (rest
))
1182 { x
= scm_lcm (x
, y
);
1184 rest
= scm_cdr (rest
);
1186 return scm_lcm (x
, y
);
1190 #define s_lcm s_scm_i_lcm
1191 #define g_lcm g_scm_i_lcm
1194 scm_lcm (SCM n1
, SCM n2
)
1196 if (SCM_UNBNDP (n2
))
1198 if (SCM_UNBNDP (n1
))
1199 return SCM_I_MAKINUM (1L);
1200 n2
= SCM_I_MAKINUM (1L);
1203 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1204 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1205 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1206 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1208 if (SCM_I_INUMP (n1
))
1210 if (SCM_I_INUMP (n2
))
1212 SCM d
= scm_gcd (n1
, n2
);
1213 if (scm_is_eq (d
, SCM_INUM0
))
1216 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1220 /* inum n1, big n2 */
1223 SCM result
= scm_i_mkbig ();
1224 scm_t_inum nn1
= SCM_I_INUM (n1
);
1225 if (nn1
== 0) return SCM_INUM0
;
1226 if (nn1
< 0) nn1
= - nn1
;
1227 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1228 scm_remember_upto_here_1 (n2
);
1236 if (SCM_I_INUMP (n2
))
1243 SCM result
= scm_i_mkbig ();
1244 mpz_lcm(SCM_I_BIG_MPZ (result
),
1246 SCM_I_BIG_MPZ (n2
));
1247 scm_remember_upto_here_2(n1
, n2
);
1248 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1254 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1259 + + + x (map digit:logand X Y)
1260 + - + x (map digit:logand X (lognot (+ -1 Y)))
1261 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1262 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1267 + + + (map digit:logior X Y)
1268 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1269 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1270 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1275 + + + (map digit:logxor X Y)
1276 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1277 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1278 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1283 + + (any digit:logand X Y)
1284 + - (any digit:logand X (lognot (+ -1 Y)))
1285 - + (any digit:logand (lognot (+ -1 X)) Y)
1290 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1291 (SCM x
, SCM y
, SCM rest
),
1292 "Return the bitwise AND of the integer arguments.\n\n"
1294 "(logand) @result{} -1\n"
1295 "(logand 7) @result{} 7\n"
1296 "(logand #b111 #b011 #b001) @result{} 1\n"
1298 #define FUNC_NAME s_scm_i_logand
1300 while (!scm_is_null (rest
))
1301 { x
= scm_logand (x
, y
);
1303 rest
= scm_cdr (rest
);
1305 return scm_logand (x
, y
);
1309 #define s_scm_logand s_scm_i_logand
1311 SCM
scm_logand (SCM n1
, SCM n2
)
1312 #define FUNC_NAME s_scm_logand
1316 if (SCM_UNBNDP (n2
))
1318 if (SCM_UNBNDP (n1
))
1319 return SCM_I_MAKINUM (-1);
1320 else if (!SCM_NUMBERP (n1
))
1321 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1322 else if (SCM_NUMBERP (n1
))
1325 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1328 if (SCM_I_INUMP (n1
))
1330 nn1
= SCM_I_INUM (n1
);
1331 if (SCM_I_INUMP (n2
))
1333 scm_t_inum nn2
= SCM_I_INUM (n2
);
1334 return SCM_I_MAKINUM (nn1
& nn2
);
1336 else if SCM_BIGP (n2
)
1342 SCM result_z
= scm_i_mkbig ();
1344 mpz_init_set_si (nn1_z
, nn1
);
1345 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1346 scm_remember_upto_here_1 (n2
);
1348 return scm_i_normbig (result_z
);
1352 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1354 else if (SCM_BIGP (n1
))
1356 if (SCM_I_INUMP (n2
))
1359 nn1
= SCM_I_INUM (n1
);
1362 else if (SCM_BIGP (n2
))
1364 SCM result_z
= scm_i_mkbig ();
1365 mpz_and (SCM_I_BIG_MPZ (result_z
),
1367 SCM_I_BIG_MPZ (n2
));
1368 scm_remember_upto_here_2 (n1
, n2
);
1369 return scm_i_normbig (result_z
);
1372 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1375 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1380 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1381 (SCM x
, SCM y
, SCM rest
),
1382 "Return the bitwise OR of the integer arguments.\n\n"
1384 "(logior) @result{} 0\n"
1385 "(logior 7) @result{} 7\n"
1386 "(logior #b000 #b001 #b011) @result{} 3\n"
1388 #define FUNC_NAME s_scm_i_logior
1390 while (!scm_is_null (rest
))
1391 { x
= scm_logior (x
, y
);
1393 rest
= scm_cdr (rest
);
1395 return scm_logior (x
, y
);
1399 #define s_scm_logior s_scm_i_logior
1401 SCM
scm_logior (SCM n1
, SCM n2
)
1402 #define FUNC_NAME s_scm_logior
1406 if (SCM_UNBNDP (n2
))
1408 if (SCM_UNBNDP (n1
))
1410 else if (SCM_NUMBERP (n1
))
1413 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1416 if (SCM_I_INUMP (n1
))
1418 nn1
= SCM_I_INUM (n1
);
1419 if (SCM_I_INUMP (n2
))
1421 long nn2
= SCM_I_INUM (n2
);
1422 return SCM_I_MAKINUM (nn1
| nn2
);
1424 else if (SCM_BIGP (n2
))
1430 SCM result_z
= scm_i_mkbig ();
1432 mpz_init_set_si (nn1_z
, nn1
);
1433 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1434 scm_remember_upto_here_1 (n2
);
1436 return scm_i_normbig (result_z
);
1440 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1442 else if (SCM_BIGP (n1
))
1444 if (SCM_I_INUMP (n2
))
1447 nn1
= SCM_I_INUM (n1
);
1450 else if (SCM_BIGP (n2
))
1452 SCM result_z
= scm_i_mkbig ();
1453 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1455 SCM_I_BIG_MPZ (n2
));
1456 scm_remember_upto_here_2 (n1
, n2
);
1457 return scm_i_normbig (result_z
);
1460 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1463 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1468 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1469 (SCM x
, SCM y
, SCM rest
),
1470 "Return the bitwise XOR of the integer arguments. A bit is\n"
1471 "set in the result if it is set in an odd number of arguments.\n"
1473 "(logxor) @result{} 0\n"
1474 "(logxor 7) @result{} 7\n"
1475 "(logxor #b000 #b001 #b011) @result{} 2\n"
1476 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1478 #define FUNC_NAME s_scm_i_logxor
1480 while (!scm_is_null (rest
))
1481 { x
= scm_logxor (x
, y
);
1483 rest
= scm_cdr (rest
);
1485 return scm_logxor (x
, y
);
1489 #define s_scm_logxor s_scm_i_logxor
1491 SCM
scm_logxor (SCM n1
, SCM n2
)
1492 #define FUNC_NAME s_scm_logxor
1496 if (SCM_UNBNDP (n2
))
1498 if (SCM_UNBNDP (n1
))
1500 else if (SCM_NUMBERP (n1
))
1503 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1506 if (SCM_I_INUMP (n1
))
1508 nn1
= SCM_I_INUM (n1
);
1509 if (SCM_I_INUMP (n2
))
1511 scm_t_inum nn2
= SCM_I_INUM (n2
);
1512 return SCM_I_MAKINUM (nn1
^ nn2
);
1514 else if (SCM_BIGP (n2
))
1518 SCM result_z
= scm_i_mkbig ();
1520 mpz_init_set_si (nn1_z
, nn1
);
1521 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1522 scm_remember_upto_here_1 (n2
);
1524 return scm_i_normbig (result_z
);
1528 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1530 else if (SCM_BIGP (n1
))
1532 if (SCM_I_INUMP (n2
))
1535 nn1
= SCM_I_INUM (n1
);
1538 else if (SCM_BIGP (n2
))
1540 SCM result_z
= scm_i_mkbig ();
1541 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1543 SCM_I_BIG_MPZ (n2
));
1544 scm_remember_upto_here_2 (n1
, n2
);
1545 return scm_i_normbig (result_z
);
1548 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1551 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1556 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1558 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1559 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1560 "without actually calculating the @code{logand}, just testing\n"
1564 "(logtest #b0100 #b1011) @result{} #f\n"
1565 "(logtest #b0100 #b0111) @result{} #t\n"
1567 #define FUNC_NAME s_scm_logtest
1571 if (SCM_I_INUMP (j
))
1573 nj
= SCM_I_INUM (j
);
1574 if (SCM_I_INUMP (k
))
1576 scm_t_inum nk
= SCM_I_INUM (k
);
1577 return scm_from_bool (nj
& nk
);
1579 else if (SCM_BIGP (k
))
1587 mpz_init_set_si (nj_z
, nj
);
1588 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1589 scm_remember_upto_here_1 (k
);
1590 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1596 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1598 else if (SCM_BIGP (j
))
1600 if (SCM_I_INUMP (k
))
1603 nj
= SCM_I_INUM (j
);
1606 else if (SCM_BIGP (k
))
1610 mpz_init (result_z
);
1614 scm_remember_upto_here_2 (j
, k
);
1615 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1616 mpz_clear (result_z
);
1620 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1623 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1628 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1630 "Test whether bit number @var{index} in @var{j} is set.\n"
1631 "@var{index} starts from 0 for the least significant bit.\n"
1634 "(logbit? 0 #b1101) @result{} #t\n"
1635 "(logbit? 1 #b1101) @result{} #f\n"
1636 "(logbit? 2 #b1101) @result{} #t\n"
1637 "(logbit? 3 #b1101) @result{} #t\n"
1638 "(logbit? 4 #b1101) @result{} #f\n"
1640 #define FUNC_NAME s_scm_logbit_p
1642 unsigned long int iindex
;
1643 iindex
= scm_to_ulong (index
);
1645 if (SCM_I_INUMP (j
))
1647 /* bits above what's in an inum follow the sign bit */
1648 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1649 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1651 else if (SCM_BIGP (j
))
1653 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1654 scm_remember_upto_here_1 (j
);
1655 return scm_from_bool (val
);
1658 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1663 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1665 "Return the integer which is the ones-complement of the integer\n"
1669 "(number->string (lognot #b10000000) 2)\n"
1670 " @result{} \"-10000001\"\n"
1671 "(number->string (lognot #b0) 2)\n"
1672 " @result{} \"-1\"\n"
1674 #define FUNC_NAME s_scm_lognot
1676 if (SCM_I_INUMP (n
)) {
1677 /* No overflow here, just need to toggle all the bits making up the inum.
1678 Enhancement: No need to strip the tag and add it back, could just xor
1679 a block of 1 bits, if that worked with the various debug versions of
1681 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1683 } else if (SCM_BIGP (n
)) {
1684 SCM result
= scm_i_mkbig ();
1685 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1686 scm_remember_upto_here_1 (n
);
1690 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1695 /* returns 0 if IN is not an integer. OUT must already be
1698 coerce_to_big (SCM in
, mpz_t out
)
1701 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1702 else if (SCM_I_INUMP (in
))
1703 mpz_set_si (out
, SCM_I_INUM (in
));
1710 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1711 (SCM n
, SCM k
, SCM m
),
1712 "Return @var{n} raised to the integer exponent\n"
1713 "@var{k}, modulo @var{m}.\n"
1716 "(modulo-expt 2 3 5)\n"
1719 #define FUNC_NAME s_scm_modulo_expt
1725 /* There are two classes of error we might encounter --
1726 1) Math errors, which we'll report by calling scm_num_overflow,
1728 2) wrong-type errors, which of course we'll report by calling
1730 We don't report those errors immediately, however; instead we do
1731 some cleanup first. These variables tell us which error (if
1732 any) we should report after cleaning up.
1734 int report_overflow
= 0;
1736 int position_of_wrong_type
= 0;
1737 SCM value_of_wrong_type
= SCM_INUM0
;
1739 SCM result
= SCM_UNDEFINED
;
1745 if (scm_is_eq (m
, SCM_INUM0
))
1747 report_overflow
= 1;
1751 if (!coerce_to_big (n
, n_tmp
))
1753 value_of_wrong_type
= n
;
1754 position_of_wrong_type
= 1;
1758 if (!coerce_to_big (k
, k_tmp
))
1760 value_of_wrong_type
= k
;
1761 position_of_wrong_type
= 2;
1765 if (!coerce_to_big (m
, m_tmp
))
1767 value_of_wrong_type
= m
;
1768 position_of_wrong_type
= 3;
1772 /* if the exponent K is negative, and we simply call mpz_powm, we
1773 will get a divide-by-zero exception when an inverse 1/n mod m
1774 doesn't exist (or is not unique). Since exceptions are hard to
1775 handle, we'll attempt the inversion "by hand" -- that way, we get
1776 a simple failure code, which is easy to handle. */
1778 if (-1 == mpz_sgn (k_tmp
))
1780 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1782 report_overflow
= 1;
1785 mpz_neg (k_tmp
, k_tmp
);
1788 result
= scm_i_mkbig ();
1789 mpz_powm (SCM_I_BIG_MPZ (result
),
1794 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1795 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1802 if (report_overflow
)
1803 scm_num_overflow (FUNC_NAME
);
1805 if (position_of_wrong_type
)
1806 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1807 value_of_wrong_type
);
1809 return scm_i_normbig (result
);
1813 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1815 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1816 "exact integer, @var{n} can be any number.\n"
1818 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1819 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1820 "includes @math{0^0} is 1.\n"
1823 "(integer-expt 2 5) @result{} 32\n"
1824 "(integer-expt -3 3) @result{} -27\n"
1825 "(integer-expt 5 -3) @result{} 1/125\n"
1826 "(integer-expt 0 0) @result{} 1\n"
1828 #define FUNC_NAME s_scm_integer_expt
1831 SCM z_i2
= SCM_BOOL_F
;
1833 SCM acc
= SCM_I_MAKINUM (1L);
1835 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1836 if (!SCM_I_INUMP (k
) && !SCM_BIGP (k
))
1837 SCM_WRONG_TYPE_ARG (2, k
);
1839 if (scm_is_true (scm_zero_p (n
)))
1841 if (scm_is_true (scm_zero_p (k
))) /* 0^0 == 1 per R5RS */
1842 return acc
; /* return exact 1, regardless of n */
1843 else if (scm_is_true (scm_positive_p (k
)))
1845 else /* return NaN for (0 ^ k) for negative k per R6RS */
1848 else if (scm_is_eq (n
, acc
))
1850 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1851 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1853 if (SCM_I_INUMP (k
))
1854 i2
= SCM_I_INUM (k
);
1855 else if (SCM_BIGP (k
))
1857 z_i2
= scm_i_clonebig (k
, 1);
1858 scm_remember_upto_here_1 (k
);
1862 SCM_WRONG_TYPE_ARG (2, k
);
1866 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1868 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1869 n
= scm_divide (n
, SCM_UNDEFINED
);
1873 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1877 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1879 return scm_product (acc
, n
);
1881 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1882 acc
= scm_product (acc
, n
);
1883 n
= scm_product (n
, n
);
1884 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1892 n
= scm_divide (n
, SCM_UNDEFINED
);
1899 return scm_product (acc
, n
);
1901 acc
= scm_product (acc
, n
);
1902 n
= scm_product (n
, n
);
1909 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1911 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1912 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1914 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1915 "@var{cnt} is negative it's a division, rounded towards negative\n"
1916 "infinity. (Note that this is not the same rounding as\n"
1917 "@code{quotient} does.)\n"
1919 "With @var{n} viewed as an infinite precision twos complement,\n"
1920 "@code{ash} means a left shift introducing zero bits, or a right\n"
1921 "shift dropping bits.\n"
1924 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1925 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1927 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1928 "(ash -23 -2) @result{} -6\n"
1930 #define FUNC_NAME s_scm_ash
1933 bits_to_shift
= scm_to_long (cnt
);
1935 if (SCM_I_INUMP (n
))
1937 scm_t_inum nn
= SCM_I_INUM (n
);
1939 if (bits_to_shift
> 0)
1941 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1942 overflow a non-zero fixnum. For smaller shifts we check the
1943 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1944 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1945 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1951 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1953 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1956 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1960 SCM result
= scm_i_inum2big (nn
);
1961 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1968 bits_to_shift
= -bits_to_shift
;
1969 if (bits_to_shift
>= SCM_LONG_BIT
)
1970 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
1972 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1976 else if (SCM_BIGP (n
))
1980 if (bits_to_shift
== 0)
1983 result
= scm_i_mkbig ();
1984 if (bits_to_shift
>= 0)
1986 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1992 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1993 we have to allocate a bignum even if the result is going to be a
1995 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1997 return scm_i_normbig (result
);
2003 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2009 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
2010 (SCM n
, SCM start
, SCM end
),
2011 "Return the integer composed of the @var{start} (inclusive)\n"
2012 "through @var{end} (exclusive) bits of @var{n}. The\n"
2013 "@var{start}th bit becomes the 0-th bit in the result.\n"
2016 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
2017 " @result{} \"1010\"\n"
2018 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
2019 " @result{} \"10110\"\n"
2021 #define FUNC_NAME s_scm_bit_extract
2023 unsigned long int istart
, iend
, bits
;
2024 istart
= scm_to_ulong (start
);
2025 iend
= scm_to_ulong (end
);
2026 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
2028 /* how many bits to keep */
2029 bits
= iend
- istart
;
2031 if (SCM_I_INUMP (n
))
2033 scm_t_inum in
= SCM_I_INUM (n
);
2035 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
2036 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
2037 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
2039 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
2041 /* Since we emulate two's complement encoded numbers, this
2042 * special case requires us to produce a result that has
2043 * more bits than can be stored in a fixnum.
2045 SCM result
= scm_i_inum2big (in
);
2046 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2051 /* mask down to requisite bits */
2052 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2053 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2055 else if (SCM_BIGP (n
))
2060 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2064 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2065 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2066 such bits into a ulong. */
2067 result
= scm_i_mkbig ();
2068 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2069 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2070 result
= scm_i_normbig (result
);
2072 scm_remember_upto_here_1 (n
);
2076 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2081 static const char scm_logtab
[] = {
2082 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2085 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2087 "Return the number of bits in integer @var{n}. If integer is\n"
2088 "positive, the 1-bits in its binary representation are counted.\n"
2089 "If negative, the 0-bits in its two's-complement binary\n"
2090 "representation are counted. If 0, 0 is returned.\n"
2093 "(logcount #b10101010)\n"
2100 #define FUNC_NAME s_scm_logcount
2102 if (SCM_I_INUMP (n
))
2104 unsigned long c
= 0;
2105 scm_t_inum nn
= SCM_I_INUM (n
);
2110 c
+= scm_logtab
[15 & nn
];
2113 return SCM_I_MAKINUM (c
);
2115 else if (SCM_BIGP (n
))
2117 unsigned long count
;
2118 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2119 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2121 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2122 scm_remember_upto_here_1 (n
);
2123 return SCM_I_MAKINUM (count
);
2126 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2131 static const char scm_ilentab
[] = {
2132 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2136 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2138 "Return the number of bits necessary to represent @var{n}.\n"
2141 "(integer-length #b10101010)\n"
2143 "(integer-length 0)\n"
2145 "(integer-length #b1111)\n"
2148 #define FUNC_NAME s_scm_integer_length
2150 if (SCM_I_INUMP (n
))
2152 unsigned long c
= 0;
2154 scm_t_inum nn
= SCM_I_INUM (n
);
2160 l
= scm_ilentab
[15 & nn
];
2163 return SCM_I_MAKINUM (c
- 4 + l
);
2165 else if (SCM_BIGP (n
))
2167 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2168 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2169 1 too big, so check for that and adjust. */
2170 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2171 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2172 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2173 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2175 scm_remember_upto_here_1 (n
);
2176 return SCM_I_MAKINUM (size
);
2179 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2183 /*** NUMBERS -> STRINGS ***/
2184 #define SCM_MAX_DBL_PREC 60
2185 #define SCM_MAX_DBL_RADIX 36
2187 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2188 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2189 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2192 void init_dblprec(int *prec
, int radix
) {
2193 /* determine floating point precision by adding successively
2194 smaller increments to 1.0 until it is considered == 1.0 */
2195 double f
= ((double)1.0)/radix
;
2196 double fsum
= 1.0 + f
;
2201 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2213 void init_fx_radix(double *fx_list
, int radix
)
2215 /* initialize a per-radix list of tolerances. When added
2216 to a number < 1.0, we can determine if we should raund
2217 up and quit converting a number to a string. */
2221 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2222 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2225 /* use this array as a way to generate a single digit */
2226 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
2229 idbl2str (double f
, char *a
, int radix
)
2231 int efmt
, dpt
, d
, i
, wp
;
2233 #ifdef DBL_MIN_10_EXP
2236 #endif /* DBL_MIN_10_EXP */
2241 radix
> SCM_MAX_DBL_RADIX
)
2243 /* revert to existing behavior */
2247 wp
= scm_dblprec
[radix
-2];
2248 fx
= fx_per_radix
[radix
-2];
2252 #ifdef HAVE_COPYSIGN
2253 double sgn
= copysign (1.0, f
);
2258 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2264 strcpy (a
, "-inf.0");
2266 strcpy (a
, "+inf.0");
2271 strcpy (a
, "+nan.0");
2281 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2282 make-uniform-vector, from causing infinite loops. */
2283 /* just do the checking...if it passes, we do the conversion for our
2284 radix again below */
2291 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2299 while (f_cpy
> 10.0)
2302 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2323 if (f
+ fx
[wp
] >= radix
)
2330 /* adding 9999 makes this equivalent to abs(x) % 3 */
2331 dpt
= (exp
+ 9999) % 3;
2335 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2357 a
[ch
++] = number_chars
[d
];
2360 if (f
+ fx
[wp
] >= 1.0)
2362 a
[ch
- 1] = number_chars
[d
+1];
2374 if ((dpt
> 4) && (exp
> 6))
2376 d
= (a
[0] == '-' ? 2 : 1);
2377 for (i
= ch
++; i
> d
; i
--)
2390 if (a
[ch
- 1] == '.')
2391 a
[ch
++] = '0'; /* trailing zero */
2400 for (i
= radix
; i
<= exp
; i
*= radix
);
2401 for (i
/= radix
; i
; i
/= radix
)
2403 a
[ch
++] = number_chars
[exp
/ i
];
2412 icmplx2str (double real
, double imag
, char *str
, int radix
)
2416 i
= idbl2str (real
, str
, radix
);
2419 /* Don't output a '+' for negative numbers or for Inf and
2420 NaN. They will provide their own sign. */
2421 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
2423 i
+= idbl2str (imag
, &str
[i
], radix
);
2430 iflo2str (SCM flt
, char *str
, int radix
)
2433 if (SCM_REALP (flt
))
2434 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2436 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2441 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2442 characters in the result.
2444 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2446 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2451 return scm_iuint2str (-num
, rad
, p
) + 1;
2454 return scm_iuint2str (num
, rad
, p
);
2457 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2458 characters in the result.
2460 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2462 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2466 scm_t_uintmax n
= num
;
2468 if (rad
< 2 || rad
> 36)
2469 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
2471 for (n
/= rad
; n
> 0; n
/= rad
)
2481 p
[i
] = number_chars
[d
];
2486 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2488 "Return a string holding the external representation of the\n"
2489 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2490 "inexact, a radix of 10 will be used.")
2491 #define FUNC_NAME s_scm_number_to_string
2495 if (SCM_UNBNDP (radix
))
2498 base
= scm_to_signed_integer (radix
, 2, 36);
2500 if (SCM_I_INUMP (n
))
2502 char num_buf
[SCM_INTBUFLEN
];
2503 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2504 return scm_from_locale_stringn (num_buf
, length
);
2506 else if (SCM_BIGP (n
))
2508 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2509 scm_remember_upto_here_1 (n
);
2510 return scm_take_locale_string (str
);
2512 else if (SCM_FRACTIONP (n
))
2514 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2515 scm_from_locale_string ("/"),
2516 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2518 else if (SCM_INEXACTP (n
))
2520 char num_buf
[FLOBUFLEN
];
2521 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2524 SCM_WRONG_TYPE_ARG (1, n
);
2529 /* These print routines used to be stubbed here so that scm_repl.c
2530 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2533 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2535 char num_buf
[FLOBUFLEN
];
2536 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2541 scm_i_print_double (double val
, SCM port
)
2543 char num_buf
[FLOBUFLEN
];
2544 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2548 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2551 char num_buf
[FLOBUFLEN
];
2552 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2557 scm_i_print_complex (double real
, double imag
, SCM port
)
2559 char num_buf
[FLOBUFLEN
];
2560 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2564 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2567 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2568 scm_display (str
, port
);
2569 scm_remember_upto_here_1 (str
);
2574 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2576 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2577 scm_remember_upto_here_1 (exp
);
2578 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2582 /*** END nums->strs ***/
2585 /*** STRINGS -> NUMBERS ***/
2587 /* The following functions implement the conversion from strings to numbers.
2588 * The implementation somehow follows the grammar for numbers as it is given
2589 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2590 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2591 * points should be noted about the implementation:
2592 * * Each function keeps a local index variable 'idx' that points at the
2593 * current position within the parsed string. The global index is only
2594 * updated if the function could parse the corresponding syntactic unit
2596 * * Similarly, the functions keep track of indicators of inexactness ('#',
2597 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2598 * global exactness information is only updated after each part has been
2599 * successfully parsed.
2600 * * Sequences of digits are parsed into temporary variables holding fixnums.
2601 * Only if these fixnums would overflow, the result variables are updated
2602 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2603 * the temporary variables holding the fixnums are cleared, and the process
2604 * starts over again. If for example fixnums were able to store five decimal
2605 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2606 * and the result was computed as 12345 * 100000 + 67890. In other words,
2607 * only every five digits two bignum operations were performed.
2610 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2612 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2614 /* Caller is responsible for checking that the return value is in range
2615 for the given radix, which should be <= 36. */
2617 char_decimal_value (scm_t_uint32 c
)
2619 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
2620 that's certainly above any valid decimal, so we take advantage of
2621 that to elide some tests. */
2622 unsigned int d
= (unsigned int) uc_decimal_value (c
);
2624 /* If that failed, try extended hexadecimals, then. Only accept ascii
2629 if (c
>= (scm_t_uint32
) 'a')
2630 d
= c
- (scm_t_uint32
)'a' + 10U;
2636 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2637 unsigned int radix
, enum t_exactness
*p_exactness
)
2639 unsigned int idx
= *p_idx
;
2640 unsigned int hash_seen
= 0;
2641 scm_t_bits shift
= 1;
2643 unsigned int digit_value
;
2646 size_t len
= scm_i_string_length (mem
);
2651 c
= scm_i_string_ref (mem
, idx
);
2652 digit_value
= char_decimal_value (c
);
2653 if (digit_value
>= radix
)
2657 result
= SCM_I_MAKINUM (digit_value
);
2660 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2670 digit_value
= char_decimal_value (c
);
2671 /* This check catches non-decimals in addition to out-of-range
2673 if (digit_value
>= radix
)
2678 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2680 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2682 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2689 shift
= shift
* radix
;
2690 add
= add
* radix
+ digit_value
;
2695 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2697 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2701 *p_exactness
= INEXACT
;
2707 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2708 * covers the parts of the rules that start at a potential point. The value
2709 * of the digits up to the point have been parsed by the caller and are given
2710 * in variable result. The content of *p_exactness indicates, whether a hash
2711 * has already been seen in the digits before the point.
2714 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2717 mem2decimal_from_point (SCM result
, SCM mem
,
2718 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2720 unsigned int idx
= *p_idx
;
2721 enum t_exactness x
= *p_exactness
;
2722 size_t len
= scm_i_string_length (mem
);
2727 if (scm_i_string_ref (mem
, idx
) == '.')
2729 scm_t_bits shift
= 1;
2731 unsigned int digit_value
;
2732 SCM big_shift
= SCM_INUM1
;
2737 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2738 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2743 digit_value
= DIGIT2UINT (c
);
2754 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2756 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2757 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2759 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2767 add
= add
* 10 + digit_value
;
2773 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2774 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2775 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2778 result
= scm_divide (result
, big_shift
);
2780 /* We've seen a decimal point, thus the value is implicitly inexact. */
2792 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2794 switch (scm_i_string_ref (mem
, idx
))
2806 c
= scm_i_string_ref (mem
, idx
);
2814 c
= scm_i_string_ref (mem
, idx
);
2823 c
= scm_i_string_ref (mem
, idx
);
2828 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2832 exponent
= DIGIT2UINT (c
);
2835 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2836 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2839 if (exponent
<= SCM_MAXEXP
)
2840 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2846 if (exponent
> SCM_MAXEXP
)
2848 size_t exp_len
= idx
- start
;
2849 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2850 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2851 scm_out_of_range ("string->number", exp_num
);
2854 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2856 result
= scm_product (result
, e
);
2858 result
= scm_divide2real (result
, e
);
2860 /* We've seen an exponent, thus the value is implicitly inexact. */
2878 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2881 mem2ureal (SCM mem
, unsigned int *p_idx
,
2882 unsigned int radix
, enum t_exactness
*p_exactness
)
2884 unsigned int idx
= *p_idx
;
2886 size_t len
= scm_i_string_length (mem
);
2888 /* Start off believing that the number will be exact. This changes
2889 to INEXACT if we see a decimal point or a hash. */
2890 enum t_exactness x
= EXACT
;
2895 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2901 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2903 /* Cobble up the fractional part. We might want to set the
2904 NaN's mantissa from it. */
2906 mem2uinteger (mem
, &idx
, 10, &x
);
2911 if (scm_i_string_ref (mem
, idx
) == '.')
2915 else if (idx
+ 1 == len
)
2917 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2920 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
2927 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2928 if (scm_is_false (uinteger
))
2933 else if (scm_i_string_ref (mem
, idx
) == '/')
2941 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2942 if (scm_is_false (divisor
))
2945 /* both are int/big here, I assume */
2946 result
= scm_i_make_ratio (uinteger
, divisor
);
2948 else if (radix
== 10)
2950 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2951 if (scm_is_false (result
))
2960 /* Update *p_exactness if the number just read was inexact. This is
2961 important for complex numbers, so that a complex number is
2962 treated as inexact overall if either its real or imaginary part
2968 /* When returning an inexact zero, make sure it is represented as a
2969 floating point value so that we can change its sign.
2971 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
2972 result
= scm_from_double (0.0);
2978 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2981 mem2complex (SCM mem
, unsigned int idx
,
2982 unsigned int radix
, enum t_exactness
*p_exactness
)
2987 size_t len
= scm_i_string_length (mem
);
2992 c
= scm_i_string_ref (mem
, idx
);
3007 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3008 if (scm_is_false (ureal
))
3010 /* input must be either +i or -i */
3015 if (scm_i_string_ref (mem
, idx
) == 'i'
3016 || scm_i_string_ref (mem
, idx
) == 'I')
3022 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
3029 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3030 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
3035 c
= scm_i_string_ref (mem
, idx
);
3039 /* either +<ureal>i or -<ureal>i */
3046 return scm_make_rectangular (SCM_INUM0
, ureal
);
3049 /* polar input: <real>@<real>. */
3060 c
= scm_i_string_ref (mem
, idx
);
3078 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3079 if (scm_is_false (angle
))
3084 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3085 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3087 result
= scm_make_polar (ureal
, angle
);
3092 /* expecting input matching <real>[+-]<ureal>?i */
3099 int sign
= (c
== '+') ? 1 : -1;
3100 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3102 if (scm_is_false (imag
))
3103 imag
= SCM_I_MAKINUM (sign
);
3104 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
3105 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3109 if (scm_i_string_ref (mem
, idx
) != 'i'
3110 && scm_i_string_ref (mem
, idx
) != 'I')
3117 return scm_make_rectangular (ureal
, imag
);
3126 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3128 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3131 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3133 unsigned int idx
= 0;
3134 unsigned int radix
= NO_RADIX
;
3135 enum t_exactness forced_x
= NO_EXACTNESS
;
3136 enum t_exactness implicit_x
= EXACT
;
3138 size_t len
= scm_i_string_length (mem
);
3140 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3141 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3143 switch (scm_i_string_ref (mem
, idx
+ 1))
3146 if (radix
!= NO_RADIX
)
3151 if (radix
!= NO_RADIX
)
3156 if (forced_x
!= NO_EXACTNESS
)
3161 if (forced_x
!= NO_EXACTNESS
)
3166 if (radix
!= NO_RADIX
)
3171 if (radix
!= NO_RADIX
)
3181 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3182 if (radix
== NO_RADIX
)
3183 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3185 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3187 if (scm_is_false (result
))
3193 if (SCM_INEXACTP (result
))
3194 return scm_inexact_to_exact (result
);
3198 if (SCM_INEXACTP (result
))
3201 return scm_exact_to_inexact (result
);
3204 if (implicit_x
== INEXACT
)
3206 if (SCM_INEXACTP (result
))
3209 return scm_exact_to_inexact (result
);
3217 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3218 unsigned int default_radix
)
3220 SCM str
= scm_from_locale_stringn (mem
, len
);
3222 return scm_i_string_to_number (str
, default_radix
);
3226 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3227 (SCM string
, SCM radix
),
3228 "Return a number of the maximally precise representation\n"
3229 "expressed by the given @var{string}. @var{radix} must be an\n"
3230 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3231 "is a default radix that may be overridden by an explicit radix\n"
3232 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3233 "supplied, then the default radix is 10. If string is not a\n"
3234 "syntactically valid notation for a number, then\n"
3235 "@code{string->number} returns @code{#f}.")
3236 #define FUNC_NAME s_scm_string_to_number
3240 SCM_VALIDATE_STRING (1, string
);
3242 if (SCM_UNBNDP (radix
))
3245 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3247 answer
= scm_i_string_to_number (string
, base
);
3248 scm_remember_upto_here_1 (string
);
3254 /*** END strs->nums ***/
3258 scm_bigequal (SCM x
, SCM y
)
3260 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3261 scm_remember_upto_here_2 (x
, y
);
3262 return scm_from_bool (0 == result
);
3266 scm_real_equalp (SCM x
, SCM y
)
3268 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3272 scm_complex_equalp (SCM x
, SCM y
)
3274 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3275 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3279 scm_i_fraction_equalp (SCM x
, SCM y
)
3281 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3282 SCM_FRACTION_NUMERATOR (y
)))
3283 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3284 SCM_FRACTION_DENOMINATOR (y
))))
3291 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3293 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3295 #define FUNC_NAME s_scm_number_p
3297 return scm_from_bool (SCM_NUMBERP (x
));
3301 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3303 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3304 "otherwise. Note that the sets of real, rational and integer\n"
3305 "values form subsets of the set of complex numbers, i. e. the\n"
3306 "predicate will also be fulfilled if @var{x} is a real,\n"
3307 "rational or integer number.")
3308 #define FUNC_NAME s_scm_complex_p
3310 /* all numbers are complex. */
3311 return scm_number_p (x
);
3315 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3317 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3318 "otherwise. Note that the set of integer values forms a subset of\n"
3319 "the set of real numbers, i. e. the predicate will also be\n"
3320 "fulfilled if @var{x} is an integer number.")
3321 #define FUNC_NAME s_scm_real_p
3323 /* we can't represent irrational numbers. */
3324 return scm_rational_p (x
);
3328 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3330 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3331 "otherwise. Note that the set of integer values forms a subset of\n"
3332 "the set of rational numbers, i. e. the predicate will also be\n"
3333 "fulfilled if @var{x} is an integer number.")
3334 #define FUNC_NAME s_scm_rational_p
3336 if (SCM_I_INUMP (x
))
3338 else if (SCM_IMP (x
))
3340 else if (SCM_BIGP (x
))
3342 else if (SCM_FRACTIONP (x
))
3344 else if (SCM_REALP (x
))
3345 /* due to their limited precision, all floating point numbers are
3346 rational as well. */
3353 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3355 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3357 #define FUNC_NAME s_scm_integer_p
3360 if (SCM_I_INUMP (x
))
3366 if (!SCM_INEXACTP (x
))
3368 if (SCM_COMPLEXP (x
))
3370 r
= SCM_REAL_VALUE (x
);
3380 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3381 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3382 (SCM x
, SCM y
, SCM rest
),
3383 "Return @code{#t} if all parameters are numerically equal.")
3384 #define FUNC_NAME s_scm_i_num_eq_p
3386 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3388 while (!scm_is_null (rest
))
3390 if (scm_is_false (scm_num_eq_p (x
, y
)))
3394 rest
= scm_cdr (rest
);
3396 return scm_num_eq_p (x
, y
);
3400 scm_num_eq_p (SCM x
, SCM y
)
3403 if (SCM_I_INUMP (x
))
3405 scm_t_signed_bits xx
= SCM_I_INUM (x
);
3406 if (SCM_I_INUMP (y
))
3408 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3409 return scm_from_bool (xx
== yy
);
3411 else if (SCM_BIGP (y
))
3413 else if (SCM_REALP (y
))
3415 /* On a 32-bit system an inum fits a double, we can cast the inum
3416 to a double and compare.
3418 But on a 64-bit system an inum is bigger than a double and
3419 casting it to a double (call that dxx) will round. dxx is at
3420 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3421 an integer and fits a long. So we cast yy to a long and
3422 compare with plain xx.
3424 An alternative (for any size system actually) would be to check
3425 yy is an integer (with floor) and is in range of an inum
3426 (compare against appropriate powers of 2) then test
3427 xx==(scm_t_signed_bits)yy. It's just a matter of which
3428 casts/comparisons might be fastest or easiest for the cpu. */
3430 double yy
= SCM_REAL_VALUE (y
);
3431 return scm_from_bool ((double) xx
== yy
3432 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3433 || xx
== (scm_t_signed_bits
) yy
));
3435 else if (SCM_COMPLEXP (y
))
3436 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3437 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3438 else if (SCM_FRACTIONP (y
))
3441 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3443 else if (SCM_BIGP (x
))
3445 if (SCM_I_INUMP (y
))
3447 else if (SCM_BIGP (y
))
3449 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3450 scm_remember_upto_here_2 (x
, y
);
3451 return scm_from_bool (0 == cmp
);
3453 else if (SCM_REALP (y
))
3456 if (isnan (SCM_REAL_VALUE (y
)))
3458 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3459 scm_remember_upto_here_1 (x
);
3460 return scm_from_bool (0 == cmp
);
3462 else if (SCM_COMPLEXP (y
))
3465 if (0.0 != SCM_COMPLEX_IMAG (y
))
3467 if (isnan (SCM_COMPLEX_REAL (y
)))
3469 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3470 scm_remember_upto_here_1 (x
);
3471 return scm_from_bool (0 == cmp
);
3473 else if (SCM_FRACTIONP (y
))
3476 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3478 else if (SCM_REALP (x
))
3480 double xx
= SCM_REAL_VALUE (x
);
3481 if (SCM_I_INUMP (y
))
3483 /* see comments with inum/real above */
3484 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3485 return scm_from_bool (xx
== (double) yy
3486 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3487 || (scm_t_signed_bits
) xx
== yy
));
3489 else if (SCM_BIGP (y
))
3492 if (isnan (SCM_REAL_VALUE (x
)))
3494 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3495 scm_remember_upto_here_1 (y
);
3496 return scm_from_bool (0 == cmp
);
3498 else if (SCM_REALP (y
))
3499 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3500 else if (SCM_COMPLEXP (y
))
3501 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3502 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3503 else if (SCM_FRACTIONP (y
))
3505 double xx
= SCM_REAL_VALUE (x
);
3509 return scm_from_bool (xx
< 0.0);
3510 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3514 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3516 else if (SCM_COMPLEXP (x
))
3518 if (SCM_I_INUMP (y
))
3519 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3520 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3521 else if (SCM_BIGP (y
))
3524 if (0.0 != SCM_COMPLEX_IMAG (x
))
3526 if (isnan (SCM_COMPLEX_REAL (x
)))
3528 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3529 scm_remember_upto_here_1 (y
);
3530 return scm_from_bool (0 == cmp
);
3532 else if (SCM_REALP (y
))
3533 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3534 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3535 else if (SCM_COMPLEXP (y
))
3536 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3537 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3538 else if (SCM_FRACTIONP (y
))
3541 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3543 xx
= SCM_COMPLEX_REAL (x
);
3547 return scm_from_bool (xx
< 0.0);
3548 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3552 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3554 else if (SCM_FRACTIONP (x
))
3556 if (SCM_I_INUMP (y
))
3558 else if (SCM_BIGP (y
))
3560 else if (SCM_REALP (y
))
3562 double yy
= SCM_REAL_VALUE (y
);
3566 return scm_from_bool (0.0 < yy
);
3567 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3570 else if (SCM_COMPLEXP (y
))
3573 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3575 yy
= SCM_COMPLEX_REAL (y
);
3579 return scm_from_bool (0.0 < yy
);
3580 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3583 else if (SCM_FRACTIONP (y
))
3584 return scm_i_fraction_equalp (x
, y
);
3586 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3589 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3593 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3594 done are good for inums, but for bignums an answer can almost always be
3595 had by just examining a few high bits of the operands, as done by GMP in
3596 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3597 of the float exponent to take into account. */
3599 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3600 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3601 (SCM x
, SCM y
, SCM rest
),
3602 "Return @code{#t} if the list of parameters is monotonically\n"
3604 #define FUNC_NAME s_scm_i_num_less_p
3606 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3608 while (!scm_is_null (rest
))
3610 if (scm_is_false (scm_less_p (x
, y
)))
3614 rest
= scm_cdr (rest
);
3616 return scm_less_p (x
, y
);
3620 scm_less_p (SCM x
, SCM y
)
3623 if (SCM_I_INUMP (x
))
3625 scm_t_inum xx
= SCM_I_INUM (x
);
3626 if (SCM_I_INUMP (y
))
3628 scm_t_inum yy
= SCM_I_INUM (y
);
3629 return scm_from_bool (xx
< yy
);
3631 else if (SCM_BIGP (y
))
3633 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3634 scm_remember_upto_here_1 (y
);
3635 return scm_from_bool (sgn
> 0);
3637 else if (SCM_REALP (y
))
3638 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3639 else if (SCM_FRACTIONP (y
))
3641 /* "x < a/b" becomes "x*b < a" */
3643 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3644 y
= SCM_FRACTION_NUMERATOR (y
);
3648 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3650 else if (SCM_BIGP (x
))
3652 if (SCM_I_INUMP (y
))
3654 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3655 scm_remember_upto_here_1 (x
);
3656 return scm_from_bool (sgn
< 0);
3658 else if (SCM_BIGP (y
))
3660 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3661 scm_remember_upto_here_2 (x
, y
);
3662 return scm_from_bool (cmp
< 0);
3664 else if (SCM_REALP (y
))
3667 if (isnan (SCM_REAL_VALUE (y
)))
3669 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3670 scm_remember_upto_here_1 (x
);
3671 return scm_from_bool (cmp
< 0);
3673 else if (SCM_FRACTIONP (y
))
3676 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3678 else if (SCM_REALP (x
))
3680 if (SCM_I_INUMP (y
))
3681 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3682 else if (SCM_BIGP (y
))
3685 if (isnan (SCM_REAL_VALUE (x
)))
3687 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3688 scm_remember_upto_here_1 (y
);
3689 return scm_from_bool (cmp
> 0);
3691 else if (SCM_REALP (y
))
3692 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3693 else if (SCM_FRACTIONP (y
))
3695 double xx
= SCM_REAL_VALUE (x
);
3699 return scm_from_bool (xx
< 0.0);
3700 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3704 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3706 else if (SCM_FRACTIONP (x
))
3708 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3710 /* "a/b < y" becomes "a < y*b" */
3711 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3712 x
= SCM_FRACTION_NUMERATOR (x
);
3715 else if (SCM_REALP (y
))
3717 double yy
= SCM_REAL_VALUE (y
);
3721 return scm_from_bool (0.0 < yy
);
3722 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3725 else if (SCM_FRACTIONP (y
))
3727 /* "a/b < c/d" becomes "a*d < c*b" */
3728 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3729 SCM_FRACTION_DENOMINATOR (y
));
3730 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3731 SCM_FRACTION_DENOMINATOR (x
));
3737 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3740 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3744 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3745 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3746 (SCM x
, SCM y
, SCM rest
),
3747 "Return @code{#t} if the list of parameters is monotonically\n"
3749 #define FUNC_NAME s_scm_i_num_gr_p
3751 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3753 while (!scm_is_null (rest
))
3755 if (scm_is_false (scm_gr_p (x
, y
)))
3759 rest
= scm_cdr (rest
);
3761 return scm_gr_p (x
, y
);
3764 #define FUNC_NAME s_scm_i_num_gr_p
3766 scm_gr_p (SCM x
, SCM y
)
3768 if (!SCM_NUMBERP (x
))
3769 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3770 else if (!SCM_NUMBERP (y
))
3771 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3773 return scm_less_p (y
, x
);
3778 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3779 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3780 (SCM x
, SCM y
, SCM rest
),
3781 "Return @code{#t} if the list of parameters is monotonically\n"
3783 #define FUNC_NAME s_scm_i_num_leq_p
3785 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3787 while (!scm_is_null (rest
))
3789 if (scm_is_false (scm_leq_p (x
, y
)))
3793 rest
= scm_cdr (rest
);
3795 return scm_leq_p (x
, y
);
3798 #define FUNC_NAME s_scm_i_num_leq_p
3800 scm_leq_p (SCM x
, SCM y
)
3802 if (!SCM_NUMBERP (x
))
3803 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3804 else if (!SCM_NUMBERP (y
))
3805 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3806 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3809 return scm_not (scm_less_p (y
, x
));
3814 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3815 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3816 (SCM x
, SCM y
, SCM rest
),
3817 "Return @code{#t} if the list of parameters is monotonically\n"
3819 #define FUNC_NAME s_scm_i_num_geq_p
3821 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3823 while (!scm_is_null (rest
))
3825 if (scm_is_false (scm_geq_p (x
, y
)))
3829 rest
= scm_cdr (rest
);
3831 return scm_geq_p (x
, y
);
3834 #define FUNC_NAME s_scm_i_num_geq_p
3836 scm_geq_p (SCM x
, SCM y
)
3838 if (!SCM_NUMBERP (x
))
3839 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3840 else if (!SCM_NUMBERP (y
))
3841 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3842 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3845 return scm_not (scm_less_p (x
, y
));
3850 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3851 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3857 if (SCM_I_INUMP (z
))
3858 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3859 else if (SCM_BIGP (z
))
3861 else if (SCM_REALP (z
))
3862 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3863 else if (SCM_COMPLEXP (z
))
3864 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3865 && SCM_COMPLEX_IMAG (z
) == 0.0);
3866 else if (SCM_FRACTIONP (z
))
3869 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3873 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3874 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3878 scm_positive_p (SCM x
)
3880 if (SCM_I_INUMP (x
))
3881 return scm_from_bool (SCM_I_INUM (x
) > 0);
3882 else if (SCM_BIGP (x
))
3884 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3885 scm_remember_upto_here_1 (x
);
3886 return scm_from_bool (sgn
> 0);
3888 else if (SCM_REALP (x
))
3889 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3890 else if (SCM_FRACTIONP (x
))
3891 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3893 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3897 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3898 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3902 scm_negative_p (SCM x
)
3904 if (SCM_I_INUMP (x
))
3905 return scm_from_bool (SCM_I_INUM (x
) < 0);
3906 else if (SCM_BIGP (x
))
3908 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3909 scm_remember_upto_here_1 (x
);
3910 return scm_from_bool (sgn
< 0);
3912 else if (SCM_REALP (x
))
3913 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3914 else if (SCM_FRACTIONP (x
))
3915 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3917 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3921 /* scm_min and scm_max return an inexact when either argument is inexact, as
3922 required by r5rs. On that basis, for exact/inexact combinations the
3923 exact is converted to inexact to compare and possibly return. This is
3924 unlike scm_less_p above which takes some trouble to preserve all bits in
3925 its test, such trouble is not required for min and max. */
3927 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3928 (SCM x
, SCM y
, SCM rest
),
3929 "Return the maximum of all parameter values.")
3930 #define FUNC_NAME s_scm_i_max
3932 while (!scm_is_null (rest
))
3933 { x
= scm_max (x
, y
);
3935 rest
= scm_cdr (rest
);
3937 return scm_max (x
, y
);
3941 #define s_max s_scm_i_max
3942 #define g_max g_scm_i_max
3945 scm_max (SCM x
, SCM y
)
3950 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3951 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3954 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3957 if (SCM_I_INUMP (x
))
3959 scm_t_inum xx
= SCM_I_INUM (x
);
3960 if (SCM_I_INUMP (y
))
3962 scm_t_inum yy
= SCM_I_INUM (y
);
3963 return (xx
< yy
) ? y
: x
;
3965 else if (SCM_BIGP (y
))
3967 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3968 scm_remember_upto_here_1 (y
);
3969 return (sgn
< 0) ? x
: y
;
3971 else if (SCM_REALP (y
))
3974 /* if y==NaN then ">" is false and we return NaN */
3975 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3977 else if (SCM_FRACTIONP (y
))
3980 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3983 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3985 else if (SCM_BIGP (x
))
3987 if (SCM_I_INUMP (y
))
3989 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3990 scm_remember_upto_here_1 (x
);
3991 return (sgn
< 0) ? y
: x
;
3993 else if (SCM_BIGP (y
))
3995 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3996 scm_remember_upto_here_2 (x
, y
);
3997 return (cmp
> 0) ? x
: y
;
3999 else if (SCM_REALP (y
))
4001 /* if y==NaN then xx>yy is false, so we return the NaN y */
4004 xx
= scm_i_big2dbl (x
);
4005 yy
= SCM_REAL_VALUE (y
);
4006 return (xx
> yy
? scm_from_double (xx
) : y
);
4008 else if (SCM_FRACTIONP (y
))
4013 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4015 else if (SCM_REALP (x
))
4017 if (SCM_I_INUMP (y
))
4019 double z
= SCM_I_INUM (y
);
4020 /* if x==NaN then "<" is false and we return NaN */
4021 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
4023 else if (SCM_BIGP (y
))
4028 else if (SCM_REALP (y
))
4030 /* if x==NaN then our explicit check means we return NaN
4031 if y==NaN then ">" is false and we return NaN
4032 calling isnan is unavoidable, since it's the only way to know
4033 which of x or y causes any compares to be false */
4034 double xx
= SCM_REAL_VALUE (x
);
4035 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
4037 else if (SCM_FRACTIONP (y
))
4039 double yy
= scm_i_fraction2double (y
);
4040 double xx
= SCM_REAL_VALUE (x
);
4041 return (xx
< yy
) ? scm_from_double (yy
) : x
;
4044 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4046 else if (SCM_FRACTIONP (x
))
4048 if (SCM_I_INUMP (y
))
4052 else if (SCM_BIGP (y
))
4056 else if (SCM_REALP (y
))
4058 double xx
= scm_i_fraction2double (x
);
4059 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4061 else if (SCM_FRACTIONP (y
))
4066 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4069 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4073 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4074 (SCM x
, SCM y
, SCM rest
),
4075 "Return the minimum of all parameter values.")
4076 #define FUNC_NAME s_scm_i_min
4078 while (!scm_is_null (rest
))
4079 { x
= scm_min (x
, y
);
4081 rest
= scm_cdr (rest
);
4083 return scm_min (x
, y
);
4087 #define s_min s_scm_i_min
4088 #define g_min g_scm_i_min
4091 scm_min (SCM x
, SCM y
)
4096 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4097 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4100 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4103 if (SCM_I_INUMP (x
))
4105 scm_t_inum xx
= SCM_I_INUM (x
);
4106 if (SCM_I_INUMP (y
))
4108 scm_t_inum yy
= SCM_I_INUM (y
);
4109 return (xx
< yy
) ? x
: y
;
4111 else if (SCM_BIGP (y
))
4113 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4114 scm_remember_upto_here_1 (y
);
4115 return (sgn
< 0) ? y
: x
;
4117 else if (SCM_REALP (y
))
4120 /* if y==NaN then "<" is false and we return NaN */
4121 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4123 else if (SCM_FRACTIONP (y
))
4126 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4129 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4131 else if (SCM_BIGP (x
))
4133 if (SCM_I_INUMP (y
))
4135 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4136 scm_remember_upto_here_1 (x
);
4137 return (sgn
< 0) ? x
: y
;
4139 else if (SCM_BIGP (y
))
4141 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4142 scm_remember_upto_here_2 (x
, y
);
4143 return (cmp
> 0) ? y
: x
;
4145 else if (SCM_REALP (y
))
4147 /* if y==NaN then xx<yy is false, so we return the NaN y */
4150 xx
= scm_i_big2dbl (x
);
4151 yy
= SCM_REAL_VALUE (y
);
4152 return (xx
< yy
? scm_from_double (xx
) : y
);
4154 else if (SCM_FRACTIONP (y
))
4159 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4161 else if (SCM_REALP (x
))
4163 if (SCM_I_INUMP (y
))
4165 double z
= SCM_I_INUM (y
);
4166 /* if x==NaN then "<" is false and we return NaN */
4167 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4169 else if (SCM_BIGP (y
))
4174 else if (SCM_REALP (y
))
4176 /* if x==NaN then our explicit check means we return NaN
4177 if y==NaN then "<" is false and we return NaN
4178 calling isnan is unavoidable, since it's the only way to know
4179 which of x or y causes any compares to be false */
4180 double xx
= SCM_REAL_VALUE (x
);
4181 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4183 else if (SCM_FRACTIONP (y
))
4185 double yy
= scm_i_fraction2double (y
);
4186 double xx
= SCM_REAL_VALUE (x
);
4187 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4190 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4192 else if (SCM_FRACTIONP (x
))
4194 if (SCM_I_INUMP (y
))
4198 else if (SCM_BIGP (y
))
4202 else if (SCM_REALP (y
))
4204 double xx
= scm_i_fraction2double (x
);
4205 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4207 else if (SCM_FRACTIONP (y
))
4212 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4215 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4219 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4220 (SCM x
, SCM y
, SCM rest
),
4221 "Return the sum of all parameter values. Return 0 if called without\n"
4223 #define FUNC_NAME s_scm_i_sum
4225 while (!scm_is_null (rest
))
4226 { x
= scm_sum (x
, y
);
4228 rest
= scm_cdr (rest
);
4230 return scm_sum (x
, y
);
4234 #define s_sum s_scm_i_sum
4235 #define g_sum g_scm_i_sum
4238 scm_sum (SCM x
, SCM y
)
4240 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4242 if (SCM_NUMBERP (x
)) return x
;
4243 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4244 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4247 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4249 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4251 scm_t_inum xx
= SCM_I_INUM (x
);
4252 scm_t_inum yy
= SCM_I_INUM (y
);
4253 scm_t_inum z
= xx
+ yy
;
4254 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
4256 else if (SCM_BIGP (y
))
4261 else if (SCM_REALP (y
))
4263 scm_t_inum xx
= SCM_I_INUM (x
);
4264 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4266 else if (SCM_COMPLEXP (y
))
4268 scm_t_inum xx
= SCM_I_INUM (x
);
4269 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4270 SCM_COMPLEX_IMAG (y
));
4272 else if (SCM_FRACTIONP (y
))
4273 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4274 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4275 SCM_FRACTION_DENOMINATOR (y
));
4277 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4278 } else if (SCM_BIGP (x
))
4280 if (SCM_I_INUMP (y
))
4285 inum
= SCM_I_INUM (y
);
4288 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4291 SCM result
= scm_i_mkbig ();
4292 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4293 scm_remember_upto_here_1 (x
);
4294 /* we know the result will have to be a bignum */
4297 return scm_i_normbig (result
);
4301 SCM result
= scm_i_mkbig ();
4302 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4303 scm_remember_upto_here_1 (x
);
4304 /* we know the result will have to be a bignum */
4307 return scm_i_normbig (result
);
4310 else if (SCM_BIGP (y
))
4312 SCM result
= scm_i_mkbig ();
4313 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4314 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4315 mpz_add (SCM_I_BIG_MPZ (result
),
4318 scm_remember_upto_here_2 (x
, y
);
4319 /* we know the result will have to be a bignum */
4322 return scm_i_normbig (result
);
4324 else if (SCM_REALP (y
))
4326 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4327 scm_remember_upto_here_1 (x
);
4328 return scm_from_double (result
);
4330 else if (SCM_COMPLEXP (y
))
4332 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4333 + SCM_COMPLEX_REAL (y
));
4334 scm_remember_upto_here_1 (x
);
4335 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4337 else if (SCM_FRACTIONP (y
))
4338 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4339 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4340 SCM_FRACTION_DENOMINATOR (y
));
4342 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4344 else if (SCM_REALP (x
))
4346 if (SCM_I_INUMP (y
))
4347 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4348 else if (SCM_BIGP (y
))
4350 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4351 scm_remember_upto_here_1 (y
);
4352 return scm_from_double (result
);
4354 else if (SCM_REALP (y
))
4355 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4356 else if (SCM_COMPLEXP (y
))
4357 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4358 SCM_COMPLEX_IMAG (y
));
4359 else if (SCM_FRACTIONP (y
))
4360 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4362 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4364 else if (SCM_COMPLEXP (x
))
4366 if (SCM_I_INUMP (y
))
4367 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4368 SCM_COMPLEX_IMAG (x
));
4369 else if (SCM_BIGP (y
))
4371 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4372 + SCM_COMPLEX_REAL (x
));
4373 scm_remember_upto_here_1 (y
);
4374 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4376 else if (SCM_REALP (y
))
4377 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4378 SCM_COMPLEX_IMAG (x
));
4379 else if (SCM_COMPLEXP (y
))
4380 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4381 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4382 else if (SCM_FRACTIONP (y
))
4383 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4384 SCM_COMPLEX_IMAG (x
));
4386 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4388 else if (SCM_FRACTIONP (x
))
4390 if (SCM_I_INUMP (y
))
4391 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4392 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4393 SCM_FRACTION_DENOMINATOR (x
));
4394 else if (SCM_BIGP (y
))
4395 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4396 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4397 SCM_FRACTION_DENOMINATOR (x
));
4398 else if (SCM_REALP (y
))
4399 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4400 else if (SCM_COMPLEXP (y
))
4401 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4402 SCM_COMPLEX_IMAG (y
));
4403 else if (SCM_FRACTIONP (y
))
4404 /* a/b + c/d = (ad + bc) / bd */
4405 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4406 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4407 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4409 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4412 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4416 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4418 "Return @math{@var{x}+1}.")
4419 #define FUNC_NAME s_scm_oneplus
4421 return scm_sum (x
, SCM_INUM1
);
4426 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4427 (SCM x
, SCM y
, SCM rest
),
4428 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4429 "the sum of all but the first argument are subtracted from the first\n"
4431 #define FUNC_NAME s_scm_i_difference
4433 while (!scm_is_null (rest
))
4434 { x
= scm_difference (x
, y
);
4436 rest
= scm_cdr (rest
);
4438 return scm_difference (x
, y
);
4442 #define s_difference s_scm_i_difference
4443 #define g_difference g_scm_i_difference
4446 scm_difference (SCM x
, SCM y
)
4447 #define FUNC_NAME s_difference
4449 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4452 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4454 if (SCM_I_INUMP (x
))
4456 scm_t_inum xx
= -SCM_I_INUM (x
);
4457 if (SCM_FIXABLE (xx
))
4458 return SCM_I_MAKINUM (xx
);
4460 return scm_i_inum2big (xx
);
4462 else if (SCM_BIGP (x
))
4463 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4464 bignum, but negating that gives a fixnum. */
4465 return scm_i_normbig (scm_i_clonebig (x
, 0));
4466 else if (SCM_REALP (x
))
4467 return scm_from_double (-SCM_REAL_VALUE (x
));
4468 else if (SCM_COMPLEXP (x
))
4469 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4470 -SCM_COMPLEX_IMAG (x
));
4471 else if (SCM_FRACTIONP (x
))
4472 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4473 SCM_FRACTION_DENOMINATOR (x
));
4475 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4478 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4480 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4482 scm_t_inum xx
= SCM_I_INUM (x
);
4483 scm_t_inum yy
= SCM_I_INUM (y
);
4484 scm_t_inum z
= xx
- yy
;
4485 if (SCM_FIXABLE (z
))
4486 return SCM_I_MAKINUM (z
);
4488 return scm_i_inum2big (z
);
4490 else if (SCM_BIGP (y
))
4492 /* inum-x - big-y */
4493 scm_t_inum xx
= SCM_I_INUM (x
);
4496 return scm_i_clonebig (y
, 0);
4499 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4500 SCM result
= scm_i_mkbig ();
4503 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4506 /* x - y == -(y + -x) */
4507 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4508 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4510 scm_remember_upto_here_1 (y
);
4512 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4513 /* we know the result will have to be a bignum */
4516 return scm_i_normbig (result
);
4519 else if (SCM_REALP (y
))
4521 scm_t_inum xx
= SCM_I_INUM (x
);
4522 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4524 else if (SCM_COMPLEXP (y
))
4526 scm_t_inum xx
= SCM_I_INUM (x
);
4527 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4528 - SCM_COMPLEX_IMAG (y
));
4530 else if (SCM_FRACTIONP (y
))
4531 /* a - b/c = (ac - b) / c */
4532 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4533 SCM_FRACTION_NUMERATOR (y
)),
4534 SCM_FRACTION_DENOMINATOR (y
));
4536 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4538 else if (SCM_BIGP (x
))
4540 if (SCM_I_INUMP (y
))
4542 /* big-x - inum-y */
4543 scm_t_inum yy
= SCM_I_INUM (y
);
4544 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4546 scm_remember_upto_here_1 (x
);
4548 return (SCM_FIXABLE (-yy
) ?
4549 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
4552 SCM result
= scm_i_mkbig ();
4555 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4557 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4558 scm_remember_upto_here_1 (x
);
4560 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4561 /* we know the result will have to be a bignum */
4564 return scm_i_normbig (result
);
4567 else if (SCM_BIGP (y
))
4569 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4570 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4571 SCM result
= scm_i_mkbig ();
4572 mpz_sub (SCM_I_BIG_MPZ (result
),
4575 scm_remember_upto_here_2 (x
, y
);
4576 /* we know the result will have to be a bignum */
4577 if ((sgn_x
== 1) && (sgn_y
== -1))
4579 if ((sgn_x
== -1) && (sgn_y
== 1))
4581 return scm_i_normbig (result
);
4583 else if (SCM_REALP (y
))
4585 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4586 scm_remember_upto_here_1 (x
);
4587 return scm_from_double (result
);
4589 else if (SCM_COMPLEXP (y
))
4591 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4592 - SCM_COMPLEX_REAL (y
));
4593 scm_remember_upto_here_1 (x
);
4594 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4596 else if (SCM_FRACTIONP (y
))
4597 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4598 SCM_FRACTION_NUMERATOR (y
)),
4599 SCM_FRACTION_DENOMINATOR (y
));
4600 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4602 else if (SCM_REALP (x
))
4604 if (SCM_I_INUMP (y
))
4605 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4606 else if (SCM_BIGP (y
))
4608 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4609 scm_remember_upto_here_1 (x
);
4610 return scm_from_double (result
);
4612 else if (SCM_REALP (y
))
4613 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4614 else if (SCM_COMPLEXP (y
))
4615 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4616 -SCM_COMPLEX_IMAG (y
));
4617 else if (SCM_FRACTIONP (y
))
4618 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4620 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4622 else if (SCM_COMPLEXP (x
))
4624 if (SCM_I_INUMP (y
))
4625 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4626 SCM_COMPLEX_IMAG (x
));
4627 else if (SCM_BIGP (y
))
4629 double real_part
= (SCM_COMPLEX_REAL (x
)
4630 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4631 scm_remember_upto_here_1 (x
);
4632 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4634 else if (SCM_REALP (y
))
4635 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4636 SCM_COMPLEX_IMAG (x
));
4637 else if (SCM_COMPLEXP (y
))
4638 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4639 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4640 else if (SCM_FRACTIONP (y
))
4641 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4642 SCM_COMPLEX_IMAG (x
));
4644 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4646 else if (SCM_FRACTIONP (x
))
4648 if (SCM_I_INUMP (y
))
4649 /* a/b - c = (a - cb) / b */
4650 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4651 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4652 SCM_FRACTION_DENOMINATOR (x
));
4653 else if (SCM_BIGP (y
))
4654 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4655 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4656 SCM_FRACTION_DENOMINATOR (x
));
4657 else if (SCM_REALP (y
))
4658 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4659 else if (SCM_COMPLEXP (y
))
4660 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4661 -SCM_COMPLEX_IMAG (y
));
4662 else if (SCM_FRACTIONP (y
))
4663 /* a/b - c/d = (ad - bc) / bd */
4664 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4665 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4666 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4668 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4671 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4676 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4678 "Return @math{@var{x}-1}.")
4679 #define FUNC_NAME s_scm_oneminus
4681 return scm_difference (x
, SCM_INUM1
);
4686 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4687 (SCM x
, SCM y
, SCM rest
),
4688 "Return the product of all arguments. If called without arguments,\n"
4690 #define FUNC_NAME s_scm_i_product
4692 while (!scm_is_null (rest
))
4693 { x
= scm_product (x
, y
);
4695 rest
= scm_cdr (rest
);
4697 return scm_product (x
, y
);
4701 #define s_product s_scm_i_product
4702 #define g_product g_scm_i_product
4705 scm_product (SCM x
, SCM y
)
4707 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4710 return SCM_I_MAKINUM (1L);
4711 else if (SCM_NUMBERP (x
))
4714 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4717 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4722 xx
= SCM_I_INUM (x
);
4726 case 0: return x
; break;
4727 case 1: return y
; break;
4730 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4732 scm_t_inum yy
= SCM_I_INUM (y
);
4733 scm_t_inum kk
= xx
* yy
;
4734 SCM k
= SCM_I_MAKINUM (kk
);
4735 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4739 SCM result
= scm_i_inum2big (xx
);
4740 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4741 return scm_i_normbig (result
);
4744 else if (SCM_BIGP (y
))
4746 SCM result
= scm_i_mkbig ();
4747 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4748 scm_remember_upto_here_1 (y
);
4751 else if (SCM_REALP (y
))
4752 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4753 else if (SCM_COMPLEXP (y
))
4754 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4755 xx
* SCM_COMPLEX_IMAG (y
));
4756 else if (SCM_FRACTIONP (y
))
4757 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4758 SCM_FRACTION_DENOMINATOR (y
));
4760 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4762 else if (SCM_BIGP (x
))
4764 if (SCM_I_INUMP (y
))
4769 else if (SCM_BIGP (y
))
4771 SCM result
= scm_i_mkbig ();
4772 mpz_mul (SCM_I_BIG_MPZ (result
),
4775 scm_remember_upto_here_2 (x
, y
);
4778 else if (SCM_REALP (y
))
4780 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4781 scm_remember_upto_here_1 (x
);
4782 return scm_from_double (result
);
4784 else if (SCM_COMPLEXP (y
))
4786 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4787 scm_remember_upto_here_1 (x
);
4788 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4789 z
* SCM_COMPLEX_IMAG (y
));
4791 else if (SCM_FRACTIONP (y
))
4792 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4793 SCM_FRACTION_DENOMINATOR (y
));
4795 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4797 else if (SCM_REALP (x
))
4799 if (SCM_I_INUMP (y
))
4801 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4802 if (scm_is_eq (y
, SCM_INUM0
))
4804 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4806 else if (SCM_BIGP (y
))
4808 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4809 scm_remember_upto_here_1 (y
);
4810 return scm_from_double (result
);
4812 else if (SCM_REALP (y
))
4813 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4814 else if (SCM_COMPLEXP (y
))
4815 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4816 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4817 else if (SCM_FRACTIONP (y
))
4818 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4820 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4822 else if (SCM_COMPLEXP (x
))
4824 if (SCM_I_INUMP (y
))
4826 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4827 if (scm_is_eq (y
, SCM_INUM0
))
4829 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4830 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4832 else if (SCM_BIGP (y
))
4834 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4835 scm_remember_upto_here_1 (y
);
4836 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4837 z
* SCM_COMPLEX_IMAG (x
));
4839 else if (SCM_REALP (y
))
4840 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4841 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4842 else if (SCM_COMPLEXP (y
))
4844 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4845 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4846 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4847 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4849 else if (SCM_FRACTIONP (y
))
4851 double yy
= scm_i_fraction2double (y
);
4852 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4853 yy
* SCM_COMPLEX_IMAG (x
));
4856 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4858 else if (SCM_FRACTIONP (x
))
4860 if (SCM_I_INUMP (y
))
4861 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4862 SCM_FRACTION_DENOMINATOR (x
));
4863 else if (SCM_BIGP (y
))
4864 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4865 SCM_FRACTION_DENOMINATOR (x
));
4866 else if (SCM_REALP (y
))
4867 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4868 else if (SCM_COMPLEXP (y
))
4870 double xx
= scm_i_fraction2double (x
);
4871 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4872 xx
* SCM_COMPLEX_IMAG (y
));
4874 else if (SCM_FRACTIONP (y
))
4875 /* a/b * c/d = ac / bd */
4876 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4877 SCM_FRACTION_NUMERATOR (y
)),
4878 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4879 SCM_FRACTION_DENOMINATOR (y
)));
4881 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4884 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4887 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4888 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4889 #define ALLOW_DIVIDE_BY_ZERO
4890 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4893 /* The code below for complex division is adapted from the GNU
4894 libstdc++, which adapted it from f2c's libF77, and is subject to
4897 /****************************************************************
4898 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4900 Permission to use, copy, modify, and distribute this software
4901 and its documentation for any purpose and without fee is hereby
4902 granted, provided that the above copyright notice appear in all
4903 copies and that both that the copyright notice and this
4904 permission notice and warranty disclaimer appear in supporting
4905 documentation, and that the names of AT&T Bell Laboratories or
4906 Bellcore or any of their entities not be used in advertising or
4907 publicity pertaining to distribution of the software without
4908 specific, written prior permission.
4910 AT&T and Bellcore disclaim all warranties with regard to this
4911 software, including all implied warranties of merchantability
4912 and fitness. In no event shall AT&T or Bellcore be liable for
4913 any special, indirect or consequential damages or any damages
4914 whatsoever resulting from loss of use, data or profits, whether
4915 in an action of contract, negligence or other tortious action,
4916 arising out of or in connection with the use or performance of
4918 ****************************************************************/
4920 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4921 (SCM x
, SCM y
, SCM rest
),
4922 "Divide the first argument by the product of the remaining\n"
4923 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4925 #define FUNC_NAME s_scm_i_divide
4927 while (!scm_is_null (rest
))
4928 { x
= scm_divide (x
, y
);
4930 rest
= scm_cdr (rest
);
4932 return scm_divide (x
, y
);
4936 #define s_divide s_scm_i_divide
4937 #define g_divide g_scm_i_divide
4940 do_divide (SCM x
, SCM y
, int inexact
)
4941 #define FUNC_NAME s_divide
4945 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4948 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4949 else if (SCM_I_INUMP (x
))
4951 scm_t_inum xx
= SCM_I_INUM (x
);
4952 if (xx
== 1 || xx
== -1)
4954 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4956 scm_num_overflow (s_divide
);
4961 return scm_from_double (1.0 / (double) xx
);
4962 else return scm_i_make_ratio (SCM_INUM1
, x
);
4965 else if (SCM_BIGP (x
))
4968 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4969 else return scm_i_make_ratio (SCM_INUM1
, x
);
4971 else if (SCM_REALP (x
))
4973 double xx
= SCM_REAL_VALUE (x
);
4974 #ifndef ALLOW_DIVIDE_BY_ZERO
4976 scm_num_overflow (s_divide
);
4979 return scm_from_double (1.0 / xx
);
4981 else if (SCM_COMPLEXP (x
))
4983 double r
= SCM_COMPLEX_REAL (x
);
4984 double i
= SCM_COMPLEX_IMAG (x
);
4985 if (fabs(r
) <= fabs(i
))
4988 double d
= i
* (1.0 + t
* t
);
4989 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4994 double d
= r
* (1.0 + t
* t
);
4995 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4998 else if (SCM_FRACTIONP (x
))
4999 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
5000 SCM_FRACTION_NUMERATOR (x
));
5002 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
5005 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5007 scm_t_inum xx
= SCM_I_INUM (x
);
5008 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5010 scm_t_inum yy
= SCM_I_INUM (y
);
5013 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5014 scm_num_overflow (s_divide
);
5016 return scm_from_double ((double) xx
/ (double) yy
);
5019 else if (xx
% yy
!= 0)
5022 return scm_from_double ((double) xx
/ (double) yy
);
5023 else return scm_i_make_ratio (x
, y
);
5027 scm_t_inum z
= xx
/ yy
;
5028 if (SCM_FIXABLE (z
))
5029 return SCM_I_MAKINUM (z
);
5031 return scm_i_inum2big (z
);
5034 else if (SCM_BIGP (y
))
5037 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
5038 else return scm_i_make_ratio (x
, y
);
5040 else if (SCM_REALP (y
))
5042 double yy
= SCM_REAL_VALUE (y
);
5043 #ifndef ALLOW_DIVIDE_BY_ZERO
5045 scm_num_overflow (s_divide
);
5048 return scm_from_double ((double) xx
/ yy
);
5050 else if (SCM_COMPLEXP (y
))
5053 complex_div
: /* y _must_ be a complex number */
5055 double r
= SCM_COMPLEX_REAL (y
);
5056 double i
= SCM_COMPLEX_IMAG (y
);
5057 if (fabs(r
) <= fabs(i
))
5060 double d
= i
* (1.0 + t
* t
);
5061 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5066 double d
= r
* (1.0 + t
* t
);
5067 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5071 else if (SCM_FRACTIONP (y
))
5072 /* a / b/c = ac / b */
5073 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5074 SCM_FRACTION_NUMERATOR (y
));
5076 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5078 else if (SCM_BIGP (x
))
5080 if (SCM_I_INUMP (y
))
5082 scm_t_inum yy
= SCM_I_INUM (y
);
5085 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5086 scm_num_overflow (s_divide
);
5088 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5089 scm_remember_upto_here_1 (x
);
5090 return (sgn
== 0) ? scm_nan () : scm_inf ();
5097 /* FIXME: HMM, what are the relative performance issues here?
5098 We need to test. Is it faster on average to test
5099 divisible_p, then perform whichever operation, or is it
5100 faster to perform the integer div opportunistically and
5101 switch to real if there's a remainder? For now we take the
5102 middle ground: test, then if divisible, use the faster div
5105 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
5106 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5110 SCM result
= scm_i_mkbig ();
5111 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5112 scm_remember_upto_here_1 (x
);
5114 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5115 return scm_i_normbig (result
);
5120 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5121 else return scm_i_make_ratio (x
, y
);
5125 else if (SCM_BIGP (y
))
5130 /* It's easily possible for the ratio x/y to fit a double
5131 but one or both x and y be too big to fit a double,
5132 hence the use of mpq_get_d rather than converting and
5135 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5136 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5137 return scm_from_double (mpq_get_d (q
));
5141 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5145 SCM result
= scm_i_mkbig ();
5146 mpz_divexact (SCM_I_BIG_MPZ (result
),
5149 scm_remember_upto_here_2 (x
, y
);
5150 return scm_i_normbig (result
);
5153 return scm_i_make_ratio (x
, y
);
5156 else if (SCM_REALP (y
))
5158 double yy
= SCM_REAL_VALUE (y
);
5159 #ifndef ALLOW_DIVIDE_BY_ZERO
5161 scm_num_overflow (s_divide
);
5164 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5166 else if (SCM_COMPLEXP (y
))
5168 a
= scm_i_big2dbl (x
);
5171 else if (SCM_FRACTIONP (y
))
5172 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5173 SCM_FRACTION_NUMERATOR (y
));
5175 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5177 else if (SCM_REALP (x
))
5179 double rx
= SCM_REAL_VALUE (x
);
5180 if (SCM_I_INUMP (y
))
5182 scm_t_inum yy
= SCM_I_INUM (y
);
5183 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5185 scm_num_overflow (s_divide
);
5188 return scm_from_double (rx
/ (double) yy
);
5190 else if (SCM_BIGP (y
))
5192 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5193 scm_remember_upto_here_1 (y
);
5194 return scm_from_double (rx
/ dby
);
5196 else if (SCM_REALP (y
))
5198 double yy
= SCM_REAL_VALUE (y
);
5199 #ifndef ALLOW_DIVIDE_BY_ZERO
5201 scm_num_overflow (s_divide
);
5204 return scm_from_double (rx
/ yy
);
5206 else if (SCM_COMPLEXP (y
))
5211 else if (SCM_FRACTIONP (y
))
5212 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5214 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5216 else if (SCM_COMPLEXP (x
))
5218 double rx
= SCM_COMPLEX_REAL (x
);
5219 double ix
= SCM_COMPLEX_IMAG (x
);
5220 if (SCM_I_INUMP (y
))
5222 scm_t_inum yy
= SCM_I_INUM (y
);
5223 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5225 scm_num_overflow (s_divide
);
5230 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5233 else if (SCM_BIGP (y
))
5235 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5236 scm_remember_upto_here_1 (y
);
5237 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5239 else if (SCM_REALP (y
))
5241 double yy
= SCM_REAL_VALUE (y
);
5242 #ifndef ALLOW_DIVIDE_BY_ZERO
5244 scm_num_overflow (s_divide
);
5247 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5249 else if (SCM_COMPLEXP (y
))
5251 double ry
= SCM_COMPLEX_REAL (y
);
5252 double iy
= SCM_COMPLEX_IMAG (y
);
5253 if (fabs(ry
) <= fabs(iy
))
5256 double d
= iy
* (1.0 + t
* t
);
5257 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5262 double d
= ry
* (1.0 + t
* t
);
5263 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5266 else if (SCM_FRACTIONP (y
))
5268 double yy
= scm_i_fraction2double (y
);
5269 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5272 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5274 else if (SCM_FRACTIONP (x
))
5276 if (SCM_I_INUMP (y
))
5278 scm_t_inum yy
= SCM_I_INUM (y
);
5279 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5281 scm_num_overflow (s_divide
);
5284 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5285 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5287 else if (SCM_BIGP (y
))
5289 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5290 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5292 else if (SCM_REALP (y
))
5294 double yy
= SCM_REAL_VALUE (y
);
5295 #ifndef ALLOW_DIVIDE_BY_ZERO
5297 scm_num_overflow (s_divide
);
5300 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5302 else if (SCM_COMPLEXP (y
))
5304 a
= scm_i_fraction2double (x
);
5307 else if (SCM_FRACTIONP (y
))
5308 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5309 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5311 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5314 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5318 scm_divide (SCM x
, SCM y
)
5320 return do_divide (x
, y
, 0);
5323 static SCM
scm_divide2real (SCM x
, SCM y
)
5325 return do_divide (x
, y
, 1);
5331 scm_c_truncate (double x
)
5342 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5343 half-way case (ie. when x is an integer plus 0.5) going upwards.
5344 Then half-way cases are identified and adjusted down if the
5345 round-upwards didn't give the desired even integer.
5347 "plus_half == result" identifies a half-way case. If plus_half, which is
5348 x + 0.5, is an integer then x must be an integer plus 0.5.
5350 An odd "result" value is identified with result/2 != floor(result/2).
5351 This is done with plus_half, since that value is ready for use sooner in
5352 a pipelined cpu, and we're already requiring plus_half == result.
5354 Note however that we need to be careful when x is big and already an
5355 integer. In that case "x+0.5" may round to an adjacent integer, causing
5356 us to return such a value, incorrectly. For instance if the hardware is
5357 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5358 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5359 returned. Or if the hardware is in round-upwards mode, then other bigger
5360 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5361 representable value, 2^128+2^76 (or whatever), again incorrect.
5363 These bad roundings of x+0.5 are avoided by testing at the start whether
5364 x is already an integer. If it is then clearly that's the desired result
5365 already. And if it's not then the exponent must be small enough to allow
5366 an 0.5 to be represented, and hence added without a bad rounding. */
5369 scm_c_round (double x
)
5371 double plus_half
, result
;
5376 plus_half
= x
+ 0.5;
5377 result
= floor (plus_half
);
5378 /* Adjust so that the rounding is towards even. */
5379 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5384 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5386 "Round the number @var{x} towards zero.")
5387 #define FUNC_NAME s_scm_truncate_number
5389 if (scm_is_false (scm_negative_p (x
)))
5390 return scm_floor (x
);
5392 return scm_ceiling (x
);
5396 static SCM exactly_one_half
;
5398 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5400 "Round the number @var{x} towards the nearest integer. "
5401 "When it is exactly halfway between two integers, "
5402 "round towards the even one.")
5403 #define FUNC_NAME s_scm_round_number
5405 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5407 else if (SCM_REALP (x
))
5408 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5411 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5412 single quotient+remainder division then examining to see which way
5413 the rounding should go. */
5414 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5415 SCM result
= scm_floor (plus_half
);
5416 /* Adjust so that the rounding is towards even. */
5417 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5418 && scm_is_true (scm_odd_p (result
)))
5419 return scm_difference (result
, SCM_INUM1
);
5426 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5428 "Round the number @var{x} towards minus infinity.")
5429 #define FUNC_NAME s_scm_floor
5431 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5433 else if (SCM_REALP (x
))
5434 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5435 else if (SCM_FRACTIONP (x
))
5437 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5438 SCM_FRACTION_DENOMINATOR (x
));
5439 if (scm_is_false (scm_negative_p (x
)))
5441 /* For positive x, rounding towards zero is correct. */
5446 /* For negative x, we need to return q-1 unless x is an
5447 integer. But fractions are never integer, per our
5449 return scm_difference (q
, SCM_INUM1
);
5453 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5457 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5459 "Round the number @var{x} towards infinity.")
5460 #define FUNC_NAME s_scm_ceiling
5462 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5464 else if (SCM_REALP (x
))
5465 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5466 else if (SCM_FRACTIONP (x
))
5468 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5469 SCM_FRACTION_DENOMINATOR (x
));
5470 if (scm_is_false (scm_positive_p (x
)))
5472 /* For negative x, rounding towards zero is correct. */
5477 /* For positive x, we need to return q+1 unless x is an
5478 integer. But fractions are never integer, per our
5480 return scm_sum (q
, SCM_INUM1
);
5484 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5488 /* sin/cos/tan/asin/acos/atan
5489 sinh/cosh/tanh/asinh/acosh/atanh
5490 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5491 Written by Jerry D. Hedden, (C) FSF.
5492 See the file `COPYING' for terms applying to this program. */
5494 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5496 "Return @var{x} raised to the power of @var{y}.")
5497 #define FUNC_NAME s_scm_expt
5499 if (scm_is_integer (y
))
5501 if (scm_is_true (scm_exact_p (y
)))
5502 return scm_integer_expt (x
, y
);
5505 /* Here we handle the case where the exponent is an inexact
5506 integer. We make the exponent exact in order to use
5507 scm_integer_expt, and thus avoid the spurious imaginary
5508 parts that may result from round-off errors in the general
5509 e^(y log x) method below (for example when squaring a large
5510 negative number). In this case, we must return an inexact
5511 result for correctness. We also make the base inexact so
5512 that scm_integer_expt will use fast inexact arithmetic
5513 internally. Note that making the base inexact is not
5514 sufficient to guarantee an inexact result, because
5515 scm_integer_expt will return an exact 1 when the exponent
5516 is 0, even if the base is inexact. */
5517 return scm_exact_to_inexact
5518 (scm_integer_expt (scm_exact_to_inexact (x
),
5519 scm_inexact_to_exact (y
)));
5522 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5524 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5527 return scm_exp (scm_product (scm_log (x
), y
));
5531 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5533 "Compute the sine of @var{z}.")
5534 #define FUNC_NAME s_scm_sin
5536 if (scm_is_real (z
))
5537 return scm_from_double (sin (scm_to_double (z
)));
5538 else if (SCM_COMPLEXP (z
))
5540 x
= SCM_COMPLEX_REAL (z
);
5541 y
= SCM_COMPLEX_IMAG (z
);
5542 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5543 cos (x
) * sinh (y
));
5546 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5550 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5552 "Compute the cosine of @var{z}.")
5553 #define FUNC_NAME s_scm_cos
5555 if (scm_is_real (z
))
5556 return scm_from_double (cos (scm_to_double (z
)));
5557 else if (SCM_COMPLEXP (z
))
5559 x
= SCM_COMPLEX_REAL (z
);
5560 y
= SCM_COMPLEX_IMAG (z
);
5561 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5562 -sin (x
) * sinh (y
));
5565 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5569 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5571 "Compute the tangent of @var{z}.")
5572 #define FUNC_NAME s_scm_tan
5574 if (scm_is_real (z
))
5575 return scm_from_double (tan (scm_to_double (z
)));
5576 else if (SCM_COMPLEXP (z
))
5578 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5579 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5580 w
= cos (x
) + cosh (y
);
5581 #ifndef ALLOW_DIVIDE_BY_ZERO
5583 scm_num_overflow (s_scm_tan
);
5585 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5588 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5592 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5594 "Compute the hyperbolic sine of @var{z}.")
5595 #define FUNC_NAME s_scm_sinh
5597 if (scm_is_real (z
))
5598 return scm_from_double (sinh (scm_to_double (z
)));
5599 else if (SCM_COMPLEXP (z
))
5601 x
= SCM_COMPLEX_REAL (z
);
5602 y
= SCM_COMPLEX_IMAG (z
);
5603 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5604 cosh (x
) * sin (y
));
5607 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5611 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5613 "Compute the hyperbolic cosine of @var{z}.")
5614 #define FUNC_NAME s_scm_cosh
5616 if (scm_is_real (z
))
5617 return scm_from_double (cosh (scm_to_double (z
)));
5618 else if (SCM_COMPLEXP (z
))
5620 x
= SCM_COMPLEX_REAL (z
);
5621 y
= SCM_COMPLEX_IMAG (z
);
5622 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5623 sinh (x
) * sin (y
));
5626 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5630 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5632 "Compute the hyperbolic tangent of @var{z}.")
5633 #define FUNC_NAME s_scm_tanh
5635 if (scm_is_real (z
))
5636 return scm_from_double (tanh (scm_to_double (z
)));
5637 else if (SCM_COMPLEXP (z
))
5639 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5640 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5641 w
= cosh (x
) + cos (y
);
5642 #ifndef ALLOW_DIVIDE_BY_ZERO
5644 scm_num_overflow (s_scm_tanh
);
5646 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5649 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5653 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5655 "Compute the arc sine of @var{z}.")
5656 #define FUNC_NAME s_scm_asin
5658 if (scm_is_real (z
))
5660 double w
= scm_to_double (z
);
5661 if (w
>= -1.0 && w
<= 1.0)
5662 return scm_from_double (asin (w
));
5664 return scm_product (scm_c_make_rectangular (0, -1),
5665 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5667 else if (SCM_COMPLEXP (z
))
5669 x
= SCM_COMPLEX_REAL (z
);
5670 y
= SCM_COMPLEX_IMAG (z
);
5671 return scm_product (scm_c_make_rectangular (0, -1),
5672 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5675 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5679 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5681 "Compute the arc cosine of @var{z}.")
5682 #define FUNC_NAME s_scm_acos
5684 if (scm_is_real (z
))
5686 double w
= scm_to_double (z
);
5687 if (w
>= -1.0 && w
<= 1.0)
5688 return scm_from_double (acos (w
));
5690 return scm_sum (scm_from_double (acos (0.0)),
5691 scm_product (scm_c_make_rectangular (0, 1),
5692 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5694 else if (SCM_COMPLEXP (z
))
5696 x
= SCM_COMPLEX_REAL (z
);
5697 y
= SCM_COMPLEX_IMAG (z
);
5698 return scm_sum (scm_from_double (acos (0.0)),
5699 scm_product (scm_c_make_rectangular (0, 1),
5700 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5703 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5707 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5709 "With one argument, compute the arc tangent of @var{z}.\n"
5710 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5711 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5712 #define FUNC_NAME s_scm_atan
5716 if (scm_is_real (z
))
5717 return scm_from_double (atan (scm_to_double (z
)));
5718 else if (SCM_COMPLEXP (z
))
5721 v
= SCM_COMPLEX_REAL (z
);
5722 w
= SCM_COMPLEX_IMAG (z
);
5723 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5724 scm_c_make_rectangular (v
, w
+ 1.0))),
5725 scm_c_make_rectangular (0, 2));
5728 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5730 else if (scm_is_real (z
))
5732 if (scm_is_real (y
))
5733 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5735 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5738 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5742 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5744 "Compute the inverse hyperbolic sine of @var{z}.")
5745 #define FUNC_NAME s_scm_sys_asinh
5747 if (scm_is_real (z
))
5748 return scm_from_double (asinh (scm_to_double (z
)));
5749 else if (scm_is_number (z
))
5750 return scm_log (scm_sum (z
,
5751 scm_sqrt (scm_sum (scm_product (z
, z
),
5754 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5758 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5760 "Compute the inverse hyperbolic cosine of @var{z}.")
5761 #define FUNC_NAME s_scm_sys_acosh
5763 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5764 return scm_from_double (acosh (scm_to_double (z
)));
5765 else if (scm_is_number (z
))
5766 return scm_log (scm_sum (z
,
5767 scm_sqrt (scm_difference (scm_product (z
, z
),
5770 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5774 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5776 "Compute the inverse hyperbolic tangent of @var{z}.")
5777 #define FUNC_NAME s_scm_sys_atanh
5779 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5780 return scm_from_double (atanh (scm_to_double (z
)));
5781 else if (scm_is_number (z
))
5782 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
5783 scm_difference (SCM_INUM1
, z
))),
5786 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5791 scm_c_make_rectangular (double re
, double im
)
5794 return scm_from_double (re
);
5799 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5801 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
5802 SCM_COMPLEX_REAL (z
) = re
;
5803 SCM_COMPLEX_IMAG (z
) = im
;
5808 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5809 (SCM real_part
, SCM imaginary_part
),
5810 "Return a complex number constructed of the given @var{real-part} "
5811 "and @var{imaginary-part} parts.")
5812 #define FUNC_NAME s_scm_make_rectangular
5814 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5815 SCM_ARG1
, FUNC_NAME
, "real");
5816 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5817 SCM_ARG2
, FUNC_NAME
, "real");
5818 return scm_c_make_rectangular (scm_to_double (real_part
),
5819 scm_to_double (imaginary_part
));
5824 scm_c_make_polar (double mag
, double ang
)
5828 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5829 use it on Glibc-based systems that have it (it's a GNU extension). See
5830 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5832 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5833 sincos (ang
, &s
, &c
);
5838 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5841 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5843 "Return the complex number @var{x} * e^(i * @var{y}).")
5844 #define FUNC_NAME s_scm_make_polar
5846 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5847 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5848 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5853 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5854 /* "Return the real part of the number @var{z}."
5857 scm_real_part (SCM z
)
5859 if (SCM_I_INUMP (z
))
5861 else if (SCM_BIGP (z
))
5863 else if (SCM_REALP (z
))
5865 else if (SCM_COMPLEXP (z
))
5866 return scm_from_double (SCM_COMPLEX_REAL (z
));
5867 else if (SCM_FRACTIONP (z
))
5870 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5874 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5875 /* "Return the imaginary part of the number @var{z}."
5878 scm_imag_part (SCM z
)
5880 if (SCM_I_INUMP (z
))
5882 else if (SCM_BIGP (z
))
5884 else if (SCM_REALP (z
))
5886 else if (SCM_COMPLEXP (z
))
5887 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5888 else if (SCM_FRACTIONP (z
))
5891 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5894 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5895 /* "Return the numerator of the number @var{z}."
5898 scm_numerator (SCM z
)
5900 if (SCM_I_INUMP (z
))
5902 else if (SCM_BIGP (z
))
5904 else if (SCM_FRACTIONP (z
))
5905 return SCM_FRACTION_NUMERATOR (z
);
5906 else if (SCM_REALP (z
))
5907 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5909 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5913 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5914 /* "Return the denominator of the number @var{z}."
5917 scm_denominator (SCM z
)
5919 if (SCM_I_INUMP (z
))
5921 else if (SCM_BIGP (z
))
5923 else if (SCM_FRACTIONP (z
))
5924 return SCM_FRACTION_DENOMINATOR (z
);
5925 else if (SCM_REALP (z
))
5926 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5928 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5931 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5932 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5933 * "@code{abs} for real arguments, but also allows complex numbers."
5936 scm_magnitude (SCM z
)
5938 if (SCM_I_INUMP (z
))
5940 scm_t_inum zz
= SCM_I_INUM (z
);
5943 else if (SCM_POSFIXABLE (-zz
))
5944 return SCM_I_MAKINUM (-zz
);
5946 return scm_i_inum2big (-zz
);
5948 else if (SCM_BIGP (z
))
5950 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5951 scm_remember_upto_here_1 (z
);
5953 return scm_i_clonebig (z
, 0);
5957 else if (SCM_REALP (z
))
5958 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5959 else if (SCM_COMPLEXP (z
))
5960 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5961 else if (SCM_FRACTIONP (z
))
5963 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5965 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5966 SCM_FRACTION_DENOMINATOR (z
));
5969 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5973 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5974 /* "Return the angle of the complex number @var{z}."
5979 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5980 flo0 to save allocating a new flonum with scm_from_double each time.
5981 But if atan2 follows the floating point rounding mode, then the value
5982 is not a constant. Maybe it'd be close enough though. */
5983 if (SCM_I_INUMP (z
))
5985 if (SCM_I_INUM (z
) >= 0)
5988 return scm_from_double (atan2 (0.0, -1.0));
5990 else if (SCM_BIGP (z
))
5992 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5993 scm_remember_upto_here_1 (z
);
5995 return scm_from_double (atan2 (0.0, -1.0));
5999 else if (SCM_REALP (z
))
6001 if (SCM_REAL_VALUE (z
) >= 0)
6004 return scm_from_double (atan2 (0.0, -1.0));
6006 else if (SCM_COMPLEXP (z
))
6007 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
6008 else if (SCM_FRACTIONP (z
))
6010 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
6012 else return scm_from_double (atan2 (0.0, -1.0));
6015 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
6019 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
6020 /* Convert the number @var{x} to its inexact representation.\n"
6023 scm_exact_to_inexact (SCM z
)
6025 if (SCM_I_INUMP (z
))
6026 return scm_from_double ((double) SCM_I_INUM (z
));
6027 else if (SCM_BIGP (z
))
6028 return scm_from_double (scm_i_big2dbl (z
));
6029 else if (SCM_FRACTIONP (z
))
6030 return scm_from_double (scm_i_fraction2double (z
));
6031 else if (SCM_INEXACTP (z
))
6034 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
6038 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
6040 "Return an exact number that is numerically closest to @var{z}.")
6041 #define FUNC_NAME s_scm_inexact_to_exact
6043 if (SCM_I_INUMP (z
))
6045 else if (SCM_BIGP (z
))
6047 else if (SCM_REALP (z
))
6049 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
6050 SCM_OUT_OF_RANGE (1, z
);
6057 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6058 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6059 scm_i_mpz2num (mpq_denref (frac
)));
6061 /* When scm_i_make_ratio throws, we leak the memory allocated
6068 else if (SCM_FRACTIONP (z
))
6071 SCM_WRONG_TYPE_ARG (1, z
);
6075 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6077 "Returns the @emph{simplest} rational number differing\n"
6078 "from @var{x} by no more than @var{eps}.\n"
6080 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6081 "exact result when both its arguments are exact. Thus, you might need\n"
6082 "to use @code{inexact->exact} on the arguments.\n"
6085 "(rationalize (inexact->exact 1.2) 1/100)\n"
6088 #define FUNC_NAME s_scm_rationalize
6090 if (SCM_I_INUMP (x
))
6092 else if (SCM_BIGP (x
))
6094 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6096 /* Use continued fractions to find closest ratio. All
6097 arithmetic is done with exact numbers.
6100 SCM ex
= scm_inexact_to_exact (x
);
6101 SCM int_part
= scm_floor (ex
);
6103 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
6104 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
6108 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6111 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6112 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6114 /* We stop after a million iterations just to be absolutely sure
6115 that we don't go into an infinite loop. The process normally
6116 converges after less than a dozen iterations.
6119 eps
= scm_abs (eps
);
6120 while (++i
< 1000000)
6122 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6123 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6124 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6126 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6127 eps
))) /* abs(x-a/b) <= eps */
6129 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6130 if (scm_is_false (scm_exact_p (x
))
6131 || scm_is_false (scm_exact_p (eps
)))
6132 return scm_exact_to_inexact (res
);
6136 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6138 tt
= scm_floor (rx
); /* tt = floor (rx) */
6144 scm_num_overflow (s_scm_rationalize
);
6147 SCM_WRONG_TYPE_ARG (1, x
);
6151 /* conversion functions */
6154 scm_is_integer (SCM val
)
6156 return scm_is_true (scm_integer_p (val
));
6160 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6162 if (SCM_I_INUMP (val
))
6164 scm_t_signed_bits n
= SCM_I_INUM (val
);
6165 return n
>= min
&& n
<= max
;
6167 else if (SCM_BIGP (val
))
6169 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6171 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6173 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6175 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6176 return n
>= min
&& n
<= max
;
6186 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6187 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6190 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6191 SCM_I_BIG_MPZ (val
));
6193 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6205 return n
>= min
&& n
<= max
;
6213 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6215 if (SCM_I_INUMP (val
))
6217 scm_t_signed_bits n
= SCM_I_INUM (val
);
6218 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6220 else if (SCM_BIGP (val
))
6222 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6224 else if (max
<= ULONG_MAX
)
6226 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6228 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6229 return n
>= min
&& n
<= max
;
6239 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6242 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6243 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6246 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6247 SCM_I_BIG_MPZ (val
));
6249 return n
>= min
&& n
<= max
;
6257 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6259 scm_error (scm_out_of_range_key
,
6261 "Value out of range ~S to ~S: ~S",
6262 scm_list_3 (min
, max
, bad_val
),
6263 scm_list_1 (bad_val
));
6266 #define TYPE scm_t_intmax
6267 #define TYPE_MIN min
6268 #define TYPE_MAX max
6269 #define SIZEOF_TYPE 0
6270 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6271 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6272 #include "libguile/conv-integer.i.c"
6274 #define TYPE scm_t_uintmax
6275 #define TYPE_MIN min
6276 #define TYPE_MAX max
6277 #define SIZEOF_TYPE 0
6278 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6279 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6280 #include "libguile/conv-uinteger.i.c"
6282 #define TYPE scm_t_int8
6283 #define TYPE_MIN SCM_T_INT8_MIN
6284 #define TYPE_MAX SCM_T_INT8_MAX
6285 #define SIZEOF_TYPE 1
6286 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6287 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6288 #include "libguile/conv-integer.i.c"
6290 #define TYPE scm_t_uint8
6292 #define TYPE_MAX SCM_T_UINT8_MAX
6293 #define SIZEOF_TYPE 1
6294 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6295 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6296 #include "libguile/conv-uinteger.i.c"
6298 #define TYPE scm_t_int16
6299 #define TYPE_MIN SCM_T_INT16_MIN
6300 #define TYPE_MAX SCM_T_INT16_MAX
6301 #define SIZEOF_TYPE 2
6302 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6303 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6304 #include "libguile/conv-integer.i.c"
6306 #define TYPE scm_t_uint16
6308 #define TYPE_MAX SCM_T_UINT16_MAX
6309 #define SIZEOF_TYPE 2
6310 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6311 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6312 #include "libguile/conv-uinteger.i.c"
6314 #define TYPE scm_t_int32
6315 #define TYPE_MIN SCM_T_INT32_MIN
6316 #define TYPE_MAX SCM_T_INT32_MAX
6317 #define SIZEOF_TYPE 4
6318 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6319 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6320 #include "libguile/conv-integer.i.c"
6322 #define TYPE scm_t_uint32
6324 #define TYPE_MAX SCM_T_UINT32_MAX
6325 #define SIZEOF_TYPE 4
6326 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6327 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6328 #include "libguile/conv-uinteger.i.c"
6330 #define TYPE scm_t_wchar
6331 #define TYPE_MIN (scm_t_int32)-1
6332 #define TYPE_MAX (scm_t_int32)0x10ffff
6333 #define SIZEOF_TYPE 4
6334 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6335 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6336 #include "libguile/conv-integer.i.c"
6338 #define TYPE scm_t_int64
6339 #define TYPE_MIN SCM_T_INT64_MIN
6340 #define TYPE_MAX SCM_T_INT64_MAX
6341 #define SIZEOF_TYPE 8
6342 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6343 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6344 #include "libguile/conv-integer.i.c"
6346 #define TYPE scm_t_uint64
6348 #define TYPE_MAX SCM_T_UINT64_MAX
6349 #define SIZEOF_TYPE 8
6350 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6351 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6352 #include "libguile/conv-uinteger.i.c"
6355 scm_to_mpz (SCM val
, mpz_t rop
)
6357 if (SCM_I_INUMP (val
))
6358 mpz_set_si (rop
, SCM_I_INUM (val
));
6359 else if (SCM_BIGP (val
))
6360 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6362 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6366 scm_from_mpz (mpz_t val
)
6368 return scm_i_mpz2num (val
);
6372 scm_is_real (SCM val
)
6374 return scm_is_true (scm_real_p (val
));
6378 scm_is_rational (SCM val
)
6380 return scm_is_true (scm_rational_p (val
));
6384 scm_to_double (SCM val
)
6386 if (SCM_I_INUMP (val
))
6387 return SCM_I_INUM (val
);
6388 else if (SCM_BIGP (val
))
6389 return scm_i_big2dbl (val
);
6390 else if (SCM_FRACTIONP (val
))
6391 return scm_i_fraction2double (val
);
6392 else if (SCM_REALP (val
))
6393 return SCM_REAL_VALUE (val
);
6395 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6399 scm_from_double (double val
)
6403 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
6405 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
6406 SCM_REAL_VALUE (z
) = val
;
6411 #if SCM_ENABLE_DEPRECATED == 1
6414 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
6416 scm_c_issue_deprecation_warning
6417 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6421 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6425 scm_out_of_range (NULL
, num
);
6428 return scm_to_double (num
);
6432 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
6434 scm_c_issue_deprecation_warning
6435 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6439 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6443 scm_out_of_range (NULL
, num
);
6446 return scm_to_double (num
);
6452 scm_is_complex (SCM val
)
6454 return scm_is_true (scm_complex_p (val
));
6458 scm_c_real_part (SCM z
)
6460 if (SCM_COMPLEXP (z
))
6461 return SCM_COMPLEX_REAL (z
);
6464 /* Use the scm_real_part to get proper error checking and
6467 return scm_to_double (scm_real_part (z
));
6472 scm_c_imag_part (SCM z
)
6474 if (SCM_COMPLEXP (z
))
6475 return SCM_COMPLEX_IMAG (z
);
6478 /* Use the scm_imag_part to get proper error checking and
6479 dispatching. The result will almost always be 0.0, but not
6482 return scm_to_double (scm_imag_part (z
));
6487 scm_c_magnitude (SCM z
)
6489 return scm_to_double (scm_magnitude (z
));
6495 return scm_to_double (scm_angle (z
));
6499 scm_is_number (SCM z
)
6501 return scm_is_true (scm_number_p (z
));
6505 /* In the following functions we dispatch to the real-arg funcs like log()
6506 when we know the arg is real, instead of just handing everything to
6507 clog() for instance. This is in case clog() doesn't optimize for a
6508 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6509 well use it to go straight to the applicable C func. */
6511 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6513 "Return the natural logarithm of @var{z}.")
6514 #define FUNC_NAME s_scm_log
6516 if (SCM_COMPLEXP (z
))
6518 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6519 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6521 double re
= SCM_COMPLEX_REAL (z
);
6522 double im
= SCM_COMPLEX_IMAG (z
);
6523 return scm_c_make_rectangular (log (hypot (re
, im
)),
6529 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6530 although the value itself overflows. */
6531 double re
= scm_to_double (z
);
6532 double l
= log (fabs (re
));
6534 return scm_from_double (l
);
6536 return scm_c_make_rectangular (l
, M_PI
);
6542 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6544 "Return the base 10 logarithm of @var{z}.")
6545 #define FUNC_NAME s_scm_log10
6547 if (SCM_COMPLEXP (z
))
6549 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6550 clog() and a multiply by M_LOG10E, rather than the fallback
6551 log10+hypot+atan2.) */
6552 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6553 && defined SCM_COMPLEX_VALUE
6554 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6556 double re
= SCM_COMPLEX_REAL (z
);
6557 double im
= SCM_COMPLEX_IMAG (z
);
6558 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6559 M_LOG10E
* atan2 (im
, re
));
6564 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6565 although the value itself overflows. */
6566 double re
= scm_to_double (z
);
6567 double l
= log10 (fabs (re
));
6569 return scm_from_double (l
);
6571 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6577 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6579 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6580 "base of natural logarithms (2.71828@dots{}).")
6581 #define FUNC_NAME s_scm_exp
6583 if (SCM_COMPLEXP (z
))
6585 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6586 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6588 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6589 SCM_COMPLEX_IMAG (z
));
6594 /* When z is a negative bignum the conversion to double overflows,
6595 giving -infinity, but that's ok, the exp is still 0.0. */
6596 return scm_from_double (exp (scm_to_double (z
)));
6602 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6604 "Return the square root of @var{z}. Of the two possible roots\n"
6605 "(positive and negative), the one with the a positive real part\n"
6606 "is returned, or if that's zero then a positive imaginary part.\n"
6610 "(sqrt 9.0) @result{} 3.0\n"
6611 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6612 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6613 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6615 #define FUNC_NAME s_scm_sqrt
6617 if (SCM_COMPLEXP (x
))
6619 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6620 && defined SCM_COMPLEX_VALUE
6621 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6623 double re
= SCM_COMPLEX_REAL (x
);
6624 double im
= SCM_COMPLEX_IMAG (x
);
6625 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6626 0.5 * atan2 (im
, re
));
6631 double xx
= scm_to_double (x
);
6633 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6635 return scm_from_double (sqrt (xx
));
6647 mpz_init_set_si (z_negative_one
, -1);
6649 /* It may be possible to tune the performance of some algorithms by using
6650 * the following constants to avoid the creation of bignums. Please, before
6651 * using these values, remember the two rules of program optimization:
6652 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6653 scm_c_define ("most-positive-fixnum",
6654 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6655 scm_c_define ("most-negative-fixnum",
6656 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6658 scm_add_feature ("complex");
6659 scm_add_feature ("inexact");
6660 flo0
= scm_from_double (0.0);
6662 /* determine floating point precision */
6663 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6665 init_dblprec(&scm_dblprec
[i
-2],i
);
6666 init_fx_radix(fx_per_radix
[i
-2],i
);
6669 /* hard code precision for base 10 if the preprocessor tells us to... */
6670 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6673 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
6674 #include "libguile/numbers.x"