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 the real number @var{x} is neither\n"
604 "infinite 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_is_real (x
))
612 SCM_WRONG_TYPE_ARG (1, x
);
616 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
618 "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
619 "@samp{-inf.0}. Otherwise return @code{#f}.")
620 #define FUNC_NAME s_scm_inf_p
623 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
624 else if (scm_is_real (x
))
627 SCM_WRONG_TYPE_ARG (1, x
);
631 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
633 "Return @code{#t} if the real number @var{x} is a NaN,\n"
634 "or @code{#f} otherwise.")
635 #define FUNC_NAME s_scm_nan_p
638 return scm_from_bool (isnan (SCM_REAL_VALUE (x
)));
639 else if (scm_is_real (x
))
642 SCM_WRONG_TYPE_ARG (1, x
);
646 /* Guile's idea of infinity. */
647 static double guile_Inf
;
649 /* Guile's idea of not a number. */
650 static double guile_NaN
;
653 guile_ieee_init (void)
655 /* Some version of gcc on some old version of Linux used to crash when
656 trying to make Inf and NaN. */
659 /* C99 INFINITY, when available.
660 FIXME: The standard allows for INFINITY to be something that overflows
661 at compile time. We ought to have a configure test to check for that
662 before trying to use it. (But in practice we believe this is not a
663 problem on any system guile is likely to target.) */
664 guile_Inf
= INFINITY
;
665 #elif defined HAVE_DINFINITY
667 extern unsigned int DINFINITY
[2];
668 guile_Inf
= (*((double *) (DINFINITY
)));
675 if (guile_Inf
== tmp
)
682 /* C99 NAN, when available */
684 #elif defined HAVE_DQNAN
687 extern unsigned int DQNAN
[2];
688 guile_NaN
= (*((double *)(DQNAN
)));
691 guile_NaN
= guile_Inf
/ guile_Inf
;
695 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
698 #define FUNC_NAME s_scm_inf
700 static int initialized
= 0;
706 return scm_from_double (guile_Inf
);
710 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
713 #define FUNC_NAME s_scm_nan
715 static int initialized
= 0;
721 return scm_from_double (guile_NaN
);
726 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
728 "Return the absolute value of @var{x}.")
733 scm_t_inum xx
= SCM_I_INUM (x
);
736 else if (SCM_POSFIXABLE (-xx
))
737 return SCM_I_MAKINUM (-xx
);
739 return scm_i_inum2big (-xx
);
741 else if (SCM_BIGP (x
))
743 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
745 return scm_i_clonebig (x
, 0);
749 else if (SCM_REALP (x
))
751 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
752 double xx
= SCM_REAL_VALUE (x
);
754 return scm_from_double (-xx
);
758 else if (SCM_FRACTIONP (x
))
760 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
762 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
763 SCM_FRACTION_DENOMINATOR (x
));
766 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
771 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
772 /* "Return the quotient of the numbers @var{x} and @var{y}."
775 scm_quotient (SCM x
, SCM y
)
779 scm_t_inum xx
= SCM_I_INUM (x
);
782 scm_t_inum yy
= SCM_I_INUM (y
);
784 scm_num_overflow (s_quotient
);
787 scm_t_inum z
= xx
/ yy
;
789 return SCM_I_MAKINUM (z
);
791 return scm_i_inum2big (z
);
794 else if (SCM_BIGP (y
))
796 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
797 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
798 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
800 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
801 scm_remember_upto_here_1 (y
);
802 return SCM_I_MAKINUM (-1);
808 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
810 else if (SCM_BIGP (x
))
814 scm_t_inum yy
= SCM_I_INUM (y
);
816 scm_num_overflow (s_quotient
);
821 SCM result
= scm_i_mkbig ();
824 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
827 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
830 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
831 scm_remember_upto_here_1 (x
);
832 return scm_i_normbig (result
);
835 else if (SCM_BIGP (y
))
837 SCM result
= scm_i_mkbig ();
838 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
841 scm_remember_upto_here_2 (x
, y
);
842 return scm_i_normbig (result
);
845 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
848 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
851 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
852 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
854 * "(remainder 13 4) @result{} 1\n"
855 * "(remainder -13 4) @result{} -1\n"
859 scm_remainder (SCM x
, SCM y
)
865 scm_t_inum yy
= SCM_I_INUM (y
);
867 scm_num_overflow (s_remainder
);
870 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
871 return SCM_I_MAKINUM (z
);
874 else if (SCM_BIGP (y
))
876 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
877 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
878 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
880 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
881 scm_remember_upto_here_1 (y
);
888 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
890 else if (SCM_BIGP (x
))
894 scm_t_inum yy
= SCM_I_INUM (y
);
896 scm_num_overflow (s_remainder
);
899 SCM result
= scm_i_mkbig ();
902 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
903 scm_remember_upto_here_1 (x
);
904 return scm_i_normbig (result
);
907 else if (SCM_BIGP (y
))
909 SCM result
= scm_i_mkbig ();
910 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
913 scm_remember_upto_here_2 (x
, y
);
914 return scm_i_normbig (result
);
917 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
920 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
924 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
925 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
927 * "(modulo 13 4) @result{} 1\n"
928 * "(modulo -13 4) @result{} 3\n"
932 scm_modulo (SCM x
, SCM y
)
936 scm_t_inum xx
= SCM_I_INUM (x
);
939 scm_t_inum yy
= SCM_I_INUM (y
);
941 scm_num_overflow (s_modulo
);
944 /* C99 specifies that "%" is the remainder corresponding to a
945 quotient rounded towards zero, and that's also traditional
946 for machine division, so z here should be well defined. */
947 scm_t_inum z
= xx
% yy
;
964 return SCM_I_MAKINUM (result
);
967 else if (SCM_BIGP (y
))
969 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
976 SCM pos_y
= scm_i_clonebig (y
, 0);
977 /* do this after the last scm_op */
978 mpz_init_set_si (z_x
, xx
);
979 result
= pos_y
; /* re-use this bignum */
980 mpz_mod (SCM_I_BIG_MPZ (result
),
982 SCM_I_BIG_MPZ (pos_y
));
983 scm_remember_upto_here_1 (pos_y
);
987 result
= scm_i_mkbig ();
988 /* do this after the last scm_op */
989 mpz_init_set_si (z_x
, xx
);
990 mpz_mod (SCM_I_BIG_MPZ (result
),
993 scm_remember_upto_here_1 (y
);
996 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
997 mpz_add (SCM_I_BIG_MPZ (result
),
999 SCM_I_BIG_MPZ (result
));
1000 scm_remember_upto_here_1 (y
);
1001 /* and do this before the next one */
1003 return scm_i_normbig (result
);
1007 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1009 else if (SCM_BIGP (x
))
1011 if (SCM_I_INUMP (y
))
1013 scm_t_inum yy
= SCM_I_INUM (y
);
1015 scm_num_overflow (s_modulo
);
1018 SCM result
= scm_i_mkbig ();
1019 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1021 (yy
< 0) ? - yy
: yy
);
1022 scm_remember_upto_here_1 (x
);
1023 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1024 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1025 SCM_I_BIG_MPZ (result
),
1027 return scm_i_normbig (result
);
1030 else if (SCM_BIGP (y
))
1033 SCM result
= scm_i_mkbig ();
1034 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1035 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1036 mpz_mod (SCM_I_BIG_MPZ (result
),
1038 SCM_I_BIG_MPZ (pos_y
));
1040 scm_remember_upto_here_1 (x
);
1041 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1042 mpz_add (SCM_I_BIG_MPZ (result
),
1044 SCM_I_BIG_MPZ (result
));
1045 scm_remember_upto_here_2 (y
, pos_y
);
1046 return scm_i_normbig (result
);
1050 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1053 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1056 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1057 (SCM x
, SCM y
, SCM rest
),
1058 "Return the greatest common divisor of all parameter values.\n"
1059 "If called without arguments, 0 is returned.")
1060 #define FUNC_NAME s_scm_i_gcd
1062 while (!scm_is_null (rest
))
1063 { x
= scm_gcd (x
, y
);
1065 rest
= scm_cdr (rest
);
1067 return scm_gcd (x
, y
);
1071 #define s_gcd s_scm_i_gcd
1072 #define g_gcd g_scm_i_gcd
1075 scm_gcd (SCM x
, SCM y
)
1078 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1080 if (SCM_I_INUMP (x
))
1082 if (SCM_I_INUMP (y
))
1084 scm_t_inum xx
= SCM_I_INUM (x
);
1085 scm_t_inum yy
= SCM_I_INUM (y
);
1086 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
1087 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
1097 /* Determine a common factor 2^k */
1098 while (!(1 & (u
| v
)))
1104 /* Now, any factor 2^n can be eliminated */
1124 return (SCM_POSFIXABLE (result
)
1125 ? SCM_I_MAKINUM (result
)
1126 : scm_i_inum2big (result
));
1128 else if (SCM_BIGP (y
))
1134 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1136 else if (SCM_BIGP (x
))
1138 if (SCM_I_INUMP (y
))
1143 yy
= SCM_I_INUM (y
);
1148 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1149 scm_remember_upto_here_1 (x
);
1150 return (SCM_POSFIXABLE (result
)
1151 ? SCM_I_MAKINUM (result
)
1152 : scm_from_unsigned_integer (result
));
1154 else if (SCM_BIGP (y
))
1156 SCM result
= scm_i_mkbig ();
1157 mpz_gcd (SCM_I_BIG_MPZ (result
),
1160 scm_remember_upto_here_2 (x
, y
);
1161 return scm_i_normbig (result
);
1164 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1167 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1170 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1171 (SCM x
, SCM y
, SCM rest
),
1172 "Return the least common multiple of the arguments.\n"
1173 "If called without arguments, 1 is returned.")
1174 #define FUNC_NAME s_scm_i_lcm
1176 while (!scm_is_null (rest
))
1177 { x
= scm_lcm (x
, y
);
1179 rest
= scm_cdr (rest
);
1181 return scm_lcm (x
, y
);
1185 #define s_lcm s_scm_i_lcm
1186 #define g_lcm g_scm_i_lcm
1189 scm_lcm (SCM n1
, SCM n2
)
1191 if (SCM_UNBNDP (n2
))
1193 if (SCM_UNBNDP (n1
))
1194 return SCM_I_MAKINUM (1L);
1195 n2
= SCM_I_MAKINUM (1L);
1198 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1199 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1200 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1201 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1203 if (SCM_I_INUMP (n1
))
1205 if (SCM_I_INUMP (n2
))
1207 SCM d
= scm_gcd (n1
, n2
);
1208 if (scm_is_eq (d
, SCM_INUM0
))
1211 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1215 /* inum n1, big n2 */
1218 SCM result
= scm_i_mkbig ();
1219 scm_t_inum nn1
= SCM_I_INUM (n1
);
1220 if (nn1
== 0) return SCM_INUM0
;
1221 if (nn1
< 0) nn1
= - nn1
;
1222 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1223 scm_remember_upto_here_1 (n2
);
1231 if (SCM_I_INUMP (n2
))
1238 SCM result
= scm_i_mkbig ();
1239 mpz_lcm(SCM_I_BIG_MPZ (result
),
1241 SCM_I_BIG_MPZ (n2
));
1242 scm_remember_upto_here_2(n1
, n2
);
1243 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1249 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1254 + + + x (map digit:logand X Y)
1255 + - + x (map digit:logand X (lognot (+ -1 Y)))
1256 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1257 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1262 + + + (map digit:logior X Y)
1263 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1264 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1265 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1270 + + + (map digit:logxor X Y)
1271 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1272 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1273 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1278 + + (any digit:logand X Y)
1279 + - (any digit:logand X (lognot (+ -1 Y)))
1280 - + (any digit:logand (lognot (+ -1 X)) Y)
1285 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1286 (SCM x
, SCM y
, SCM rest
),
1287 "Return the bitwise AND of the integer arguments.\n\n"
1289 "(logand) @result{} -1\n"
1290 "(logand 7) @result{} 7\n"
1291 "(logand #b111 #b011 #b001) @result{} 1\n"
1293 #define FUNC_NAME s_scm_i_logand
1295 while (!scm_is_null (rest
))
1296 { x
= scm_logand (x
, y
);
1298 rest
= scm_cdr (rest
);
1300 return scm_logand (x
, y
);
1304 #define s_scm_logand s_scm_i_logand
1306 SCM
scm_logand (SCM n1
, SCM n2
)
1307 #define FUNC_NAME s_scm_logand
1311 if (SCM_UNBNDP (n2
))
1313 if (SCM_UNBNDP (n1
))
1314 return SCM_I_MAKINUM (-1);
1315 else if (!SCM_NUMBERP (n1
))
1316 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1317 else if (SCM_NUMBERP (n1
))
1320 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1323 if (SCM_I_INUMP (n1
))
1325 nn1
= SCM_I_INUM (n1
);
1326 if (SCM_I_INUMP (n2
))
1328 scm_t_inum nn2
= SCM_I_INUM (n2
);
1329 return SCM_I_MAKINUM (nn1
& nn2
);
1331 else if SCM_BIGP (n2
)
1337 SCM result_z
= scm_i_mkbig ();
1339 mpz_init_set_si (nn1_z
, nn1
);
1340 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1341 scm_remember_upto_here_1 (n2
);
1343 return scm_i_normbig (result_z
);
1347 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1349 else if (SCM_BIGP (n1
))
1351 if (SCM_I_INUMP (n2
))
1354 nn1
= SCM_I_INUM (n1
);
1357 else if (SCM_BIGP (n2
))
1359 SCM result_z
= scm_i_mkbig ();
1360 mpz_and (SCM_I_BIG_MPZ (result_z
),
1362 SCM_I_BIG_MPZ (n2
));
1363 scm_remember_upto_here_2 (n1
, n2
);
1364 return scm_i_normbig (result_z
);
1367 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1370 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1375 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1376 (SCM x
, SCM y
, SCM rest
),
1377 "Return the bitwise OR of the integer arguments.\n\n"
1379 "(logior) @result{} 0\n"
1380 "(logior 7) @result{} 7\n"
1381 "(logior #b000 #b001 #b011) @result{} 3\n"
1383 #define FUNC_NAME s_scm_i_logior
1385 while (!scm_is_null (rest
))
1386 { x
= scm_logior (x
, y
);
1388 rest
= scm_cdr (rest
);
1390 return scm_logior (x
, y
);
1394 #define s_scm_logior s_scm_i_logior
1396 SCM
scm_logior (SCM n1
, SCM n2
)
1397 #define FUNC_NAME s_scm_logior
1401 if (SCM_UNBNDP (n2
))
1403 if (SCM_UNBNDP (n1
))
1405 else if (SCM_NUMBERP (n1
))
1408 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1411 if (SCM_I_INUMP (n1
))
1413 nn1
= SCM_I_INUM (n1
);
1414 if (SCM_I_INUMP (n2
))
1416 long nn2
= SCM_I_INUM (n2
);
1417 return SCM_I_MAKINUM (nn1
| nn2
);
1419 else if (SCM_BIGP (n2
))
1425 SCM result_z
= scm_i_mkbig ();
1427 mpz_init_set_si (nn1_z
, nn1
);
1428 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1429 scm_remember_upto_here_1 (n2
);
1431 return scm_i_normbig (result_z
);
1435 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1437 else if (SCM_BIGP (n1
))
1439 if (SCM_I_INUMP (n2
))
1442 nn1
= SCM_I_INUM (n1
);
1445 else if (SCM_BIGP (n2
))
1447 SCM result_z
= scm_i_mkbig ();
1448 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1450 SCM_I_BIG_MPZ (n2
));
1451 scm_remember_upto_here_2 (n1
, n2
);
1452 return scm_i_normbig (result_z
);
1455 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1458 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1463 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1464 (SCM x
, SCM y
, SCM rest
),
1465 "Return the bitwise XOR of the integer arguments. A bit is\n"
1466 "set in the result if it is set in an odd number of arguments.\n"
1468 "(logxor) @result{} 0\n"
1469 "(logxor 7) @result{} 7\n"
1470 "(logxor #b000 #b001 #b011) @result{} 2\n"
1471 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1473 #define FUNC_NAME s_scm_i_logxor
1475 while (!scm_is_null (rest
))
1476 { x
= scm_logxor (x
, y
);
1478 rest
= scm_cdr (rest
);
1480 return scm_logxor (x
, y
);
1484 #define s_scm_logxor s_scm_i_logxor
1486 SCM
scm_logxor (SCM n1
, SCM n2
)
1487 #define FUNC_NAME s_scm_logxor
1491 if (SCM_UNBNDP (n2
))
1493 if (SCM_UNBNDP (n1
))
1495 else if (SCM_NUMBERP (n1
))
1498 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1501 if (SCM_I_INUMP (n1
))
1503 nn1
= SCM_I_INUM (n1
);
1504 if (SCM_I_INUMP (n2
))
1506 scm_t_inum nn2
= SCM_I_INUM (n2
);
1507 return SCM_I_MAKINUM (nn1
^ nn2
);
1509 else if (SCM_BIGP (n2
))
1513 SCM result_z
= scm_i_mkbig ();
1515 mpz_init_set_si (nn1_z
, nn1
);
1516 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1517 scm_remember_upto_here_1 (n2
);
1519 return scm_i_normbig (result_z
);
1523 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1525 else if (SCM_BIGP (n1
))
1527 if (SCM_I_INUMP (n2
))
1530 nn1
= SCM_I_INUM (n1
);
1533 else if (SCM_BIGP (n2
))
1535 SCM result_z
= scm_i_mkbig ();
1536 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1538 SCM_I_BIG_MPZ (n2
));
1539 scm_remember_upto_here_2 (n1
, n2
);
1540 return scm_i_normbig (result_z
);
1543 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1546 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1551 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1553 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1554 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1555 "without actually calculating the @code{logand}, just testing\n"
1559 "(logtest #b0100 #b1011) @result{} #f\n"
1560 "(logtest #b0100 #b0111) @result{} #t\n"
1562 #define FUNC_NAME s_scm_logtest
1566 if (SCM_I_INUMP (j
))
1568 nj
= SCM_I_INUM (j
);
1569 if (SCM_I_INUMP (k
))
1571 scm_t_inum nk
= SCM_I_INUM (k
);
1572 return scm_from_bool (nj
& nk
);
1574 else if (SCM_BIGP (k
))
1582 mpz_init_set_si (nj_z
, nj
);
1583 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1584 scm_remember_upto_here_1 (k
);
1585 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1591 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1593 else if (SCM_BIGP (j
))
1595 if (SCM_I_INUMP (k
))
1598 nj
= SCM_I_INUM (j
);
1601 else if (SCM_BIGP (k
))
1605 mpz_init (result_z
);
1609 scm_remember_upto_here_2 (j
, k
);
1610 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1611 mpz_clear (result_z
);
1615 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1618 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1623 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1625 "Test whether bit number @var{index} in @var{j} is set.\n"
1626 "@var{index} starts from 0 for the least significant bit.\n"
1629 "(logbit? 0 #b1101) @result{} #t\n"
1630 "(logbit? 1 #b1101) @result{} #f\n"
1631 "(logbit? 2 #b1101) @result{} #t\n"
1632 "(logbit? 3 #b1101) @result{} #t\n"
1633 "(logbit? 4 #b1101) @result{} #f\n"
1635 #define FUNC_NAME s_scm_logbit_p
1637 unsigned long int iindex
;
1638 iindex
= scm_to_ulong (index
);
1640 if (SCM_I_INUMP (j
))
1642 /* bits above what's in an inum follow the sign bit */
1643 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1644 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1646 else if (SCM_BIGP (j
))
1648 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1649 scm_remember_upto_here_1 (j
);
1650 return scm_from_bool (val
);
1653 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1658 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1660 "Return the integer which is the ones-complement of the integer\n"
1664 "(number->string (lognot #b10000000) 2)\n"
1665 " @result{} \"-10000001\"\n"
1666 "(number->string (lognot #b0) 2)\n"
1667 " @result{} \"-1\"\n"
1669 #define FUNC_NAME s_scm_lognot
1671 if (SCM_I_INUMP (n
)) {
1672 /* No overflow here, just need to toggle all the bits making up the inum.
1673 Enhancement: No need to strip the tag and add it back, could just xor
1674 a block of 1 bits, if that worked with the various debug versions of
1676 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1678 } else if (SCM_BIGP (n
)) {
1679 SCM result
= scm_i_mkbig ();
1680 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1681 scm_remember_upto_here_1 (n
);
1685 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1690 /* returns 0 if IN is not an integer. OUT must already be
1693 coerce_to_big (SCM in
, mpz_t out
)
1696 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1697 else if (SCM_I_INUMP (in
))
1698 mpz_set_si (out
, SCM_I_INUM (in
));
1705 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1706 (SCM n
, SCM k
, SCM m
),
1707 "Return @var{n} raised to the integer exponent\n"
1708 "@var{k}, modulo @var{m}.\n"
1711 "(modulo-expt 2 3 5)\n"
1714 #define FUNC_NAME s_scm_modulo_expt
1720 /* There are two classes of error we might encounter --
1721 1) Math errors, which we'll report by calling scm_num_overflow,
1723 2) wrong-type errors, which of course we'll report by calling
1725 We don't report those errors immediately, however; instead we do
1726 some cleanup first. These variables tell us which error (if
1727 any) we should report after cleaning up.
1729 int report_overflow
= 0;
1731 int position_of_wrong_type
= 0;
1732 SCM value_of_wrong_type
= SCM_INUM0
;
1734 SCM result
= SCM_UNDEFINED
;
1740 if (scm_is_eq (m
, SCM_INUM0
))
1742 report_overflow
= 1;
1746 if (!coerce_to_big (n
, n_tmp
))
1748 value_of_wrong_type
= n
;
1749 position_of_wrong_type
= 1;
1753 if (!coerce_to_big (k
, k_tmp
))
1755 value_of_wrong_type
= k
;
1756 position_of_wrong_type
= 2;
1760 if (!coerce_to_big (m
, m_tmp
))
1762 value_of_wrong_type
= m
;
1763 position_of_wrong_type
= 3;
1767 /* if the exponent K is negative, and we simply call mpz_powm, we
1768 will get a divide-by-zero exception when an inverse 1/n mod m
1769 doesn't exist (or is not unique). Since exceptions are hard to
1770 handle, we'll attempt the inversion "by hand" -- that way, we get
1771 a simple failure code, which is easy to handle. */
1773 if (-1 == mpz_sgn (k_tmp
))
1775 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1777 report_overflow
= 1;
1780 mpz_neg (k_tmp
, k_tmp
);
1783 result
= scm_i_mkbig ();
1784 mpz_powm (SCM_I_BIG_MPZ (result
),
1789 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1790 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1797 if (report_overflow
)
1798 scm_num_overflow (FUNC_NAME
);
1800 if (position_of_wrong_type
)
1801 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1802 value_of_wrong_type
);
1804 return scm_i_normbig (result
);
1808 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1810 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1811 "exact integer, @var{n} can be any number.\n"
1813 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1814 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1815 "includes @math{0^0} is 1.\n"
1818 "(integer-expt 2 5) @result{} 32\n"
1819 "(integer-expt -3 3) @result{} -27\n"
1820 "(integer-expt 5 -3) @result{} 1/125\n"
1821 "(integer-expt 0 0) @result{} 1\n"
1823 #define FUNC_NAME s_scm_integer_expt
1826 SCM z_i2
= SCM_BOOL_F
;
1828 SCM acc
= SCM_I_MAKINUM (1L);
1830 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1831 if (!SCM_I_INUMP (k
) && !SCM_BIGP (k
))
1832 SCM_WRONG_TYPE_ARG (2, k
);
1834 if (scm_is_true (scm_zero_p (n
)))
1836 if (scm_is_true (scm_zero_p (k
))) /* 0^0 == 1 per R5RS */
1837 return acc
; /* return exact 1, regardless of n */
1838 else if (scm_is_true (scm_positive_p (k
)))
1840 else /* return NaN for (0 ^ k) for negative k per R6RS */
1843 else if (scm_is_eq (n
, acc
))
1845 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1846 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1848 if (SCM_I_INUMP (k
))
1849 i2
= SCM_I_INUM (k
);
1850 else if (SCM_BIGP (k
))
1852 z_i2
= scm_i_clonebig (k
, 1);
1853 scm_remember_upto_here_1 (k
);
1857 SCM_WRONG_TYPE_ARG (2, k
);
1861 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1863 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1864 n
= scm_divide (n
, SCM_UNDEFINED
);
1868 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1872 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1874 return scm_product (acc
, n
);
1876 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1877 acc
= scm_product (acc
, n
);
1878 n
= scm_product (n
, n
);
1879 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1887 n
= scm_divide (n
, SCM_UNDEFINED
);
1894 return scm_product (acc
, n
);
1896 acc
= scm_product (acc
, n
);
1897 n
= scm_product (n
, n
);
1904 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1906 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1907 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1909 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1910 "@var{cnt} is negative it's a division, rounded towards negative\n"
1911 "infinity. (Note that this is not the same rounding as\n"
1912 "@code{quotient} does.)\n"
1914 "With @var{n} viewed as an infinite precision twos complement,\n"
1915 "@code{ash} means a left shift introducing zero bits, or a right\n"
1916 "shift dropping bits.\n"
1919 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1920 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1922 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1923 "(ash -23 -2) @result{} -6\n"
1925 #define FUNC_NAME s_scm_ash
1928 bits_to_shift
= scm_to_long (cnt
);
1930 if (SCM_I_INUMP (n
))
1932 scm_t_inum nn
= SCM_I_INUM (n
);
1934 if (bits_to_shift
> 0)
1936 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1937 overflow a non-zero fixnum. For smaller shifts we check the
1938 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1939 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1940 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1946 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1948 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1951 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1955 SCM result
= scm_i_inum2big (nn
);
1956 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1963 bits_to_shift
= -bits_to_shift
;
1964 if (bits_to_shift
>= SCM_LONG_BIT
)
1965 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
1967 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1971 else if (SCM_BIGP (n
))
1975 if (bits_to_shift
== 0)
1978 result
= scm_i_mkbig ();
1979 if (bits_to_shift
>= 0)
1981 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1987 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1988 we have to allocate a bignum even if the result is going to be a
1990 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1992 return scm_i_normbig (result
);
1998 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2004 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
2005 (SCM n
, SCM start
, SCM end
),
2006 "Return the integer composed of the @var{start} (inclusive)\n"
2007 "through @var{end} (exclusive) bits of @var{n}. The\n"
2008 "@var{start}th bit becomes the 0-th bit in the result.\n"
2011 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
2012 " @result{} \"1010\"\n"
2013 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
2014 " @result{} \"10110\"\n"
2016 #define FUNC_NAME s_scm_bit_extract
2018 unsigned long int istart
, iend
, bits
;
2019 istart
= scm_to_ulong (start
);
2020 iend
= scm_to_ulong (end
);
2021 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
2023 /* how many bits to keep */
2024 bits
= iend
- istart
;
2026 if (SCM_I_INUMP (n
))
2028 scm_t_inum in
= SCM_I_INUM (n
);
2030 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
2031 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
2032 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
2034 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
2036 /* Since we emulate two's complement encoded numbers, this
2037 * special case requires us to produce a result that has
2038 * more bits than can be stored in a fixnum.
2040 SCM result
= scm_i_inum2big (in
);
2041 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2046 /* mask down to requisite bits */
2047 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2048 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2050 else if (SCM_BIGP (n
))
2055 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2059 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2060 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2061 such bits into a ulong. */
2062 result
= scm_i_mkbig ();
2063 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2064 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2065 result
= scm_i_normbig (result
);
2067 scm_remember_upto_here_1 (n
);
2071 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2076 static const char scm_logtab
[] = {
2077 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2080 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2082 "Return the number of bits in integer @var{n}. If integer is\n"
2083 "positive, the 1-bits in its binary representation are counted.\n"
2084 "If negative, the 0-bits in its two's-complement binary\n"
2085 "representation are counted. If 0, 0 is returned.\n"
2088 "(logcount #b10101010)\n"
2095 #define FUNC_NAME s_scm_logcount
2097 if (SCM_I_INUMP (n
))
2099 unsigned long c
= 0;
2100 scm_t_inum nn
= SCM_I_INUM (n
);
2105 c
+= scm_logtab
[15 & nn
];
2108 return SCM_I_MAKINUM (c
);
2110 else if (SCM_BIGP (n
))
2112 unsigned long count
;
2113 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2114 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2116 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2117 scm_remember_upto_here_1 (n
);
2118 return SCM_I_MAKINUM (count
);
2121 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2126 static const char scm_ilentab
[] = {
2127 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2131 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2133 "Return the number of bits necessary to represent @var{n}.\n"
2136 "(integer-length #b10101010)\n"
2138 "(integer-length 0)\n"
2140 "(integer-length #b1111)\n"
2143 #define FUNC_NAME s_scm_integer_length
2145 if (SCM_I_INUMP (n
))
2147 unsigned long c
= 0;
2149 scm_t_inum nn
= SCM_I_INUM (n
);
2155 l
= scm_ilentab
[15 & nn
];
2158 return SCM_I_MAKINUM (c
- 4 + l
);
2160 else if (SCM_BIGP (n
))
2162 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2163 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2164 1 too big, so check for that and adjust. */
2165 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2166 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2167 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2168 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2170 scm_remember_upto_here_1 (n
);
2171 return SCM_I_MAKINUM (size
);
2174 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2178 /*** NUMBERS -> STRINGS ***/
2179 #define SCM_MAX_DBL_PREC 60
2180 #define SCM_MAX_DBL_RADIX 36
2182 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2183 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2184 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2187 void init_dblprec(int *prec
, int radix
) {
2188 /* determine floating point precision by adding successively
2189 smaller increments to 1.0 until it is considered == 1.0 */
2190 double f
= ((double)1.0)/radix
;
2191 double fsum
= 1.0 + f
;
2196 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2208 void init_fx_radix(double *fx_list
, int radix
)
2210 /* initialize a per-radix list of tolerances. When added
2211 to a number < 1.0, we can determine if we should raund
2212 up and quit converting a number to a string. */
2216 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2217 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2220 /* use this array as a way to generate a single digit */
2221 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
2224 idbl2str (double f
, char *a
, int radix
)
2226 int efmt
, dpt
, d
, i
, wp
;
2228 #ifdef DBL_MIN_10_EXP
2231 #endif /* DBL_MIN_10_EXP */
2236 radix
> SCM_MAX_DBL_RADIX
)
2238 /* revert to existing behavior */
2242 wp
= scm_dblprec
[radix
-2];
2243 fx
= fx_per_radix
[radix
-2];
2247 #ifdef HAVE_COPYSIGN
2248 double sgn
= copysign (1.0, f
);
2253 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2259 strcpy (a
, "-inf.0");
2261 strcpy (a
, "+inf.0");
2266 strcpy (a
, "+nan.0");
2276 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2277 make-uniform-vector, from causing infinite loops. */
2278 /* just do the checking...if it passes, we do the conversion for our
2279 radix again below */
2286 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2294 while (f_cpy
> 10.0)
2297 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2318 if (f
+ fx
[wp
] >= radix
)
2325 /* adding 9999 makes this equivalent to abs(x) % 3 */
2326 dpt
= (exp
+ 9999) % 3;
2330 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2352 a
[ch
++] = number_chars
[d
];
2355 if (f
+ fx
[wp
] >= 1.0)
2357 a
[ch
- 1] = number_chars
[d
+1];
2369 if ((dpt
> 4) && (exp
> 6))
2371 d
= (a
[0] == '-' ? 2 : 1);
2372 for (i
= ch
++; i
> d
; i
--)
2385 if (a
[ch
- 1] == '.')
2386 a
[ch
++] = '0'; /* trailing zero */
2395 for (i
= radix
; i
<= exp
; i
*= radix
);
2396 for (i
/= radix
; i
; i
/= radix
)
2398 a
[ch
++] = number_chars
[exp
/ i
];
2407 icmplx2str (double real
, double imag
, char *str
, int radix
)
2411 i
= idbl2str (real
, str
, radix
);
2414 /* Don't output a '+' for negative numbers or for Inf and
2415 NaN. They will provide their own sign. */
2416 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
2418 i
+= idbl2str (imag
, &str
[i
], radix
);
2425 iflo2str (SCM flt
, char *str
, int radix
)
2428 if (SCM_REALP (flt
))
2429 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2431 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2436 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2437 characters in the result.
2439 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2441 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2446 return scm_iuint2str (-num
, rad
, p
) + 1;
2449 return scm_iuint2str (num
, rad
, p
);
2452 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2453 characters in the result.
2455 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2457 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2461 scm_t_uintmax n
= num
;
2463 if (rad
< 2 || rad
> 36)
2464 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
2466 for (n
/= rad
; n
> 0; n
/= rad
)
2476 p
[i
] = number_chars
[d
];
2481 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2483 "Return a string holding the external representation of the\n"
2484 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2485 "inexact, a radix of 10 will be used.")
2486 #define FUNC_NAME s_scm_number_to_string
2490 if (SCM_UNBNDP (radix
))
2493 base
= scm_to_signed_integer (radix
, 2, 36);
2495 if (SCM_I_INUMP (n
))
2497 char num_buf
[SCM_INTBUFLEN
];
2498 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2499 return scm_from_locale_stringn (num_buf
, length
);
2501 else if (SCM_BIGP (n
))
2503 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2504 scm_remember_upto_here_1 (n
);
2505 return scm_take_locale_string (str
);
2507 else if (SCM_FRACTIONP (n
))
2509 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2510 scm_from_locale_string ("/"),
2511 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2513 else if (SCM_INEXACTP (n
))
2515 char num_buf
[FLOBUFLEN
];
2516 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2519 SCM_WRONG_TYPE_ARG (1, n
);
2524 /* These print routines used to be stubbed here so that scm_repl.c
2525 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2528 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2530 char num_buf
[FLOBUFLEN
];
2531 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2536 scm_i_print_double (double val
, SCM port
)
2538 char num_buf
[FLOBUFLEN
];
2539 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2543 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2546 char num_buf
[FLOBUFLEN
];
2547 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2552 scm_i_print_complex (double real
, double imag
, SCM port
)
2554 char num_buf
[FLOBUFLEN
];
2555 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2559 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2562 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2563 scm_display (str
, port
);
2564 scm_remember_upto_here_1 (str
);
2569 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2571 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2572 scm_remember_upto_here_1 (exp
);
2573 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2577 /*** END nums->strs ***/
2580 /*** STRINGS -> NUMBERS ***/
2582 /* The following functions implement the conversion from strings to numbers.
2583 * The implementation somehow follows the grammar for numbers as it is given
2584 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2585 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2586 * points should be noted about the implementation:
2587 * * Each function keeps a local index variable 'idx' that points at the
2588 * current position within the parsed string. The global index is only
2589 * updated if the function could parse the corresponding syntactic unit
2591 * * Similarly, the functions keep track of indicators of inexactness ('#',
2592 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2593 * global exactness information is only updated after each part has been
2594 * successfully parsed.
2595 * * Sequences of digits are parsed into temporary variables holding fixnums.
2596 * Only if these fixnums would overflow, the result variables are updated
2597 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2598 * the temporary variables holding the fixnums are cleared, and the process
2599 * starts over again. If for example fixnums were able to store five decimal
2600 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2601 * and the result was computed as 12345 * 100000 + 67890. In other words,
2602 * only every five digits two bignum operations were performed.
2605 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2607 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2609 /* Caller is responsible for checking that the return value is in range
2610 for the given radix, which should be <= 36. */
2612 char_decimal_value (scm_t_uint32 c
)
2614 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
2615 that's certainly above any valid decimal, so we take advantage of
2616 that to elide some tests. */
2617 unsigned int d
= (unsigned int) uc_decimal_value (c
);
2619 /* If that failed, try extended hexadecimals, then. Only accept ascii
2624 if (c
>= (scm_t_uint32
) 'a')
2625 d
= c
- (scm_t_uint32
)'a' + 10U;
2631 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2632 unsigned int radix
, enum t_exactness
*p_exactness
)
2634 unsigned int idx
= *p_idx
;
2635 unsigned int hash_seen
= 0;
2636 scm_t_bits shift
= 1;
2638 unsigned int digit_value
;
2641 size_t len
= scm_i_string_length (mem
);
2646 c
= scm_i_string_ref (mem
, idx
);
2647 digit_value
= char_decimal_value (c
);
2648 if (digit_value
>= radix
)
2652 result
= SCM_I_MAKINUM (digit_value
);
2655 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2665 digit_value
= char_decimal_value (c
);
2666 /* This check catches non-decimals in addition to out-of-range
2668 if (digit_value
>= radix
)
2673 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2675 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2677 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2684 shift
= shift
* radix
;
2685 add
= add
* radix
+ digit_value
;
2690 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2692 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2696 *p_exactness
= INEXACT
;
2702 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2703 * covers the parts of the rules that start at a potential point. The value
2704 * of the digits up to the point have been parsed by the caller and are given
2705 * in variable result. The content of *p_exactness indicates, whether a hash
2706 * has already been seen in the digits before the point.
2709 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2712 mem2decimal_from_point (SCM result
, SCM mem
,
2713 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2715 unsigned int idx
= *p_idx
;
2716 enum t_exactness x
= *p_exactness
;
2717 size_t len
= scm_i_string_length (mem
);
2722 if (scm_i_string_ref (mem
, idx
) == '.')
2724 scm_t_bits shift
= 1;
2726 unsigned int digit_value
;
2727 SCM big_shift
= SCM_INUM1
;
2732 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2733 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2738 digit_value
= DIGIT2UINT (c
);
2749 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2751 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2752 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2754 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2762 add
= add
* 10 + digit_value
;
2768 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2769 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2770 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2773 result
= scm_divide (result
, big_shift
);
2775 /* We've seen a decimal point, thus the value is implicitly inexact. */
2787 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2789 switch (scm_i_string_ref (mem
, idx
))
2801 c
= scm_i_string_ref (mem
, idx
);
2809 c
= scm_i_string_ref (mem
, idx
);
2818 c
= scm_i_string_ref (mem
, idx
);
2823 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2827 exponent
= DIGIT2UINT (c
);
2830 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2831 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2834 if (exponent
<= SCM_MAXEXP
)
2835 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2841 if (exponent
> SCM_MAXEXP
)
2843 size_t exp_len
= idx
- start
;
2844 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2845 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2846 scm_out_of_range ("string->number", exp_num
);
2849 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2851 result
= scm_product (result
, e
);
2853 result
= scm_divide2real (result
, e
);
2855 /* We've seen an exponent, thus the value is implicitly inexact. */
2873 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2876 mem2ureal (SCM mem
, unsigned int *p_idx
,
2877 unsigned int radix
, enum t_exactness
*p_exactness
)
2879 unsigned int idx
= *p_idx
;
2881 size_t len
= scm_i_string_length (mem
);
2883 /* Start off believing that the number will be exact. This changes
2884 to INEXACT if we see a decimal point or a hash. */
2885 enum t_exactness x
= EXACT
;
2890 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2896 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2898 /* Cobble up the fractional part. We might want to set the
2899 NaN's mantissa from it. */
2901 mem2uinteger (mem
, &idx
, 10, &x
);
2906 if (scm_i_string_ref (mem
, idx
) == '.')
2910 else if (idx
+ 1 == len
)
2912 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2915 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
2922 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2923 if (scm_is_false (uinteger
))
2928 else if (scm_i_string_ref (mem
, idx
) == '/')
2936 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2937 if (scm_is_false (divisor
))
2940 /* both are int/big here, I assume */
2941 result
= scm_i_make_ratio (uinteger
, divisor
);
2943 else if (radix
== 10)
2945 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2946 if (scm_is_false (result
))
2955 /* Update *p_exactness if the number just read was inexact. This is
2956 important for complex numbers, so that a complex number is
2957 treated as inexact overall if either its real or imaginary part
2963 /* When returning an inexact zero, make sure it is represented as a
2964 floating point value so that we can change its sign.
2966 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
2967 result
= scm_from_double (0.0);
2973 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2976 mem2complex (SCM mem
, unsigned int idx
,
2977 unsigned int radix
, enum t_exactness
*p_exactness
)
2982 size_t len
= scm_i_string_length (mem
);
2987 c
= scm_i_string_ref (mem
, idx
);
3002 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3003 if (scm_is_false (ureal
))
3005 /* input must be either +i or -i */
3010 if (scm_i_string_ref (mem
, idx
) == 'i'
3011 || scm_i_string_ref (mem
, idx
) == 'I')
3017 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
3024 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3025 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
3030 c
= scm_i_string_ref (mem
, idx
);
3034 /* either +<ureal>i or -<ureal>i */
3041 return scm_make_rectangular (SCM_INUM0
, ureal
);
3044 /* polar input: <real>@<real>. */
3055 c
= scm_i_string_ref (mem
, idx
);
3073 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3074 if (scm_is_false (angle
))
3079 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3080 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3082 result
= scm_make_polar (ureal
, angle
);
3087 /* expecting input matching <real>[+-]<ureal>?i */
3094 int sign
= (c
== '+') ? 1 : -1;
3095 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3097 if (scm_is_false (imag
))
3098 imag
= SCM_I_MAKINUM (sign
);
3099 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
3100 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3104 if (scm_i_string_ref (mem
, idx
) != 'i'
3105 && scm_i_string_ref (mem
, idx
) != 'I')
3112 return scm_make_rectangular (ureal
, imag
);
3121 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3123 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3126 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3128 unsigned int idx
= 0;
3129 unsigned int radix
= NO_RADIX
;
3130 enum t_exactness forced_x
= NO_EXACTNESS
;
3131 enum t_exactness implicit_x
= EXACT
;
3133 size_t len
= scm_i_string_length (mem
);
3135 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3136 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3138 switch (scm_i_string_ref (mem
, idx
+ 1))
3141 if (radix
!= NO_RADIX
)
3146 if (radix
!= NO_RADIX
)
3151 if (forced_x
!= NO_EXACTNESS
)
3156 if (forced_x
!= NO_EXACTNESS
)
3161 if (radix
!= NO_RADIX
)
3166 if (radix
!= NO_RADIX
)
3176 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3177 if (radix
== NO_RADIX
)
3178 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3180 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3182 if (scm_is_false (result
))
3188 if (SCM_INEXACTP (result
))
3189 return scm_inexact_to_exact (result
);
3193 if (SCM_INEXACTP (result
))
3196 return scm_exact_to_inexact (result
);
3199 if (implicit_x
== INEXACT
)
3201 if (SCM_INEXACTP (result
))
3204 return scm_exact_to_inexact (result
);
3212 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3213 unsigned int default_radix
)
3215 SCM str
= scm_from_locale_stringn (mem
, len
);
3217 return scm_i_string_to_number (str
, default_radix
);
3221 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3222 (SCM string
, SCM radix
),
3223 "Return a number of the maximally precise representation\n"
3224 "expressed by the given @var{string}. @var{radix} must be an\n"
3225 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3226 "is a default radix that may be overridden by an explicit radix\n"
3227 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3228 "supplied, then the default radix is 10. If string is not a\n"
3229 "syntactically valid notation for a number, then\n"
3230 "@code{string->number} returns @code{#f}.")
3231 #define FUNC_NAME s_scm_string_to_number
3235 SCM_VALIDATE_STRING (1, string
);
3237 if (SCM_UNBNDP (radix
))
3240 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3242 answer
= scm_i_string_to_number (string
, base
);
3243 scm_remember_upto_here_1 (string
);
3249 /*** END strs->nums ***/
3252 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3254 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3256 #define FUNC_NAME s_scm_number_p
3258 return scm_from_bool (SCM_NUMBERP (x
));
3262 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3264 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3265 "otherwise. Note that the sets of real, rational and integer\n"
3266 "values form subsets of the set of complex numbers, i. e. the\n"
3267 "predicate will also be fulfilled if @var{x} is a real,\n"
3268 "rational or integer number.")
3269 #define FUNC_NAME s_scm_complex_p
3271 /* all numbers are complex. */
3272 return scm_number_p (x
);
3276 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3278 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3279 "otherwise. Note that the set of integer values forms a subset of\n"
3280 "the set of real numbers, i. e. the predicate will also be\n"
3281 "fulfilled if @var{x} is an integer number.")
3282 #define FUNC_NAME s_scm_real_p
3284 return scm_from_bool
3285 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
3289 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3291 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3292 "otherwise. Note that the set of integer values forms a subset of\n"
3293 "the set of rational numbers, i. e. the predicate will also be\n"
3294 "fulfilled if @var{x} is an integer number.")
3295 #define FUNC_NAME s_scm_rational_p
3297 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
3299 else if (SCM_REALP (x
))
3300 /* due to their limited precision, finite floating point numbers are
3301 rational as well. (finite means neither infinity nor a NaN) */
3302 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
3308 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3310 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3312 #define FUNC_NAME s_scm_integer_p
3314 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
3316 else if (SCM_REALP (x
))
3318 double val
= SCM_REAL_VALUE (x
);
3319 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
3327 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3328 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3329 (SCM x
, SCM y
, SCM rest
),
3330 "Return @code{#t} if all parameters are numerically equal.")
3331 #define FUNC_NAME s_scm_i_num_eq_p
3333 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3335 while (!scm_is_null (rest
))
3337 if (scm_is_false (scm_num_eq_p (x
, y
)))
3341 rest
= scm_cdr (rest
);
3343 return scm_num_eq_p (x
, y
);
3347 scm_num_eq_p (SCM x
, SCM y
)
3350 if (SCM_I_INUMP (x
))
3352 scm_t_signed_bits xx
= SCM_I_INUM (x
);
3353 if (SCM_I_INUMP (y
))
3355 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3356 return scm_from_bool (xx
== yy
);
3358 else if (SCM_BIGP (y
))
3360 else if (SCM_REALP (y
))
3362 /* On a 32-bit system an inum fits a double, we can cast the inum
3363 to a double and compare.
3365 But on a 64-bit system an inum is bigger than a double and
3366 casting it to a double (call that dxx) will round. dxx is at
3367 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3368 an integer and fits a long. So we cast yy to a long and
3369 compare with plain xx.
3371 An alternative (for any size system actually) would be to check
3372 yy is an integer (with floor) and is in range of an inum
3373 (compare against appropriate powers of 2) then test
3374 xx==(scm_t_signed_bits)yy. It's just a matter of which
3375 casts/comparisons might be fastest or easiest for the cpu. */
3377 double yy
= SCM_REAL_VALUE (y
);
3378 return scm_from_bool ((double) xx
== yy
3379 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3380 || xx
== (scm_t_signed_bits
) yy
));
3382 else if (SCM_COMPLEXP (y
))
3383 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3384 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3385 else if (SCM_FRACTIONP (y
))
3388 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3390 else if (SCM_BIGP (x
))
3392 if (SCM_I_INUMP (y
))
3394 else if (SCM_BIGP (y
))
3396 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3397 scm_remember_upto_here_2 (x
, y
);
3398 return scm_from_bool (0 == cmp
);
3400 else if (SCM_REALP (y
))
3403 if (isnan (SCM_REAL_VALUE (y
)))
3405 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3406 scm_remember_upto_here_1 (x
);
3407 return scm_from_bool (0 == cmp
);
3409 else if (SCM_COMPLEXP (y
))
3412 if (0.0 != SCM_COMPLEX_IMAG (y
))
3414 if (isnan (SCM_COMPLEX_REAL (y
)))
3416 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3417 scm_remember_upto_here_1 (x
);
3418 return scm_from_bool (0 == cmp
);
3420 else if (SCM_FRACTIONP (y
))
3423 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3425 else if (SCM_REALP (x
))
3427 double xx
= SCM_REAL_VALUE (x
);
3428 if (SCM_I_INUMP (y
))
3430 /* see comments with inum/real above */
3431 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3432 return scm_from_bool (xx
== (double) yy
3433 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3434 || (scm_t_signed_bits
) xx
== yy
));
3436 else if (SCM_BIGP (y
))
3439 if (isnan (SCM_REAL_VALUE (x
)))
3441 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3442 scm_remember_upto_here_1 (y
);
3443 return scm_from_bool (0 == cmp
);
3445 else if (SCM_REALP (y
))
3446 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3447 else if (SCM_COMPLEXP (y
))
3448 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3449 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3450 else if (SCM_FRACTIONP (y
))
3452 double xx
= SCM_REAL_VALUE (x
);
3456 return scm_from_bool (xx
< 0.0);
3457 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3461 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3463 else if (SCM_COMPLEXP (x
))
3465 if (SCM_I_INUMP (y
))
3466 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3467 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3468 else if (SCM_BIGP (y
))
3471 if (0.0 != SCM_COMPLEX_IMAG (x
))
3473 if (isnan (SCM_COMPLEX_REAL (x
)))
3475 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3476 scm_remember_upto_here_1 (y
);
3477 return scm_from_bool (0 == cmp
);
3479 else if (SCM_REALP (y
))
3480 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3481 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3482 else if (SCM_COMPLEXP (y
))
3483 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3484 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3485 else if (SCM_FRACTIONP (y
))
3488 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3490 xx
= SCM_COMPLEX_REAL (x
);
3494 return scm_from_bool (xx
< 0.0);
3495 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3499 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3501 else if (SCM_FRACTIONP (x
))
3503 if (SCM_I_INUMP (y
))
3505 else if (SCM_BIGP (y
))
3507 else if (SCM_REALP (y
))
3509 double yy
= SCM_REAL_VALUE (y
);
3513 return scm_from_bool (0.0 < yy
);
3514 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3517 else if (SCM_COMPLEXP (y
))
3520 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3522 yy
= SCM_COMPLEX_REAL (y
);
3526 return scm_from_bool (0.0 < yy
);
3527 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3530 else if (SCM_FRACTIONP (y
))
3531 return scm_i_fraction_equalp (x
, y
);
3533 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3536 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3540 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3541 done are good for inums, but for bignums an answer can almost always be
3542 had by just examining a few high bits of the operands, as done by GMP in
3543 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3544 of the float exponent to take into account. */
3546 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3547 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3548 (SCM x
, SCM y
, SCM rest
),
3549 "Return @code{#t} if the list of parameters is monotonically\n"
3551 #define FUNC_NAME s_scm_i_num_less_p
3553 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3555 while (!scm_is_null (rest
))
3557 if (scm_is_false (scm_less_p (x
, y
)))
3561 rest
= scm_cdr (rest
);
3563 return scm_less_p (x
, y
);
3567 scm_less_p (SCM x
, SCM y
)
3570 if (SCM_I_INUMP (x
))
3572 scm_t_inum xx
= SCM_I_INUM (x
);
3573 if (SCM_I_INUMP (y
))
3575 scm_t_inum yy
= SCM_I_INUM (y
);
3576 return scm_from_bool (xx
< yy
);
3578 else if (SCM_BIGP (y
))
3580 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3581 scm_remember_upto_here_1 (y
);
3582 return scm_from_bool (sgn
> 0);
3584 else if (SCM_REALP (y
))
3585 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3586 else if (SCM_FRACTIONP (y
))
3588 /* "x < a/b" becomes "x*b < a" */
3590 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3591 y
= SCM_FRACTION_NUMERATOR (y
);
3595 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3597 else if (SCM_BIGP (x
))
3599 if (SCM_I_INUMP (y
))
3601 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3602 scm_remember_upto_here_1 (x
);
3603 return scm_from_bool (sgn
< 0);
3605 else if (SCM_BIGP (y
))
3607 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3608 scm_remember_upto_here_2 (x
, y
);
3609 return scm_from_bool (cmp
< 0);
3611 else if (SCM_REALP (y
))
3614 if (isnan (SCM_REAL_VALUE (y
)))
3616 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3617 scm_remember_upto_here_1 (x
);
3618 return scm_from_bool (cmp
< 0);
3620 else if (SCM_FRACTIONP (y
))
3623 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3625 else if (SCM_REALP (x
))
3627 if (SCM_I_INUMP (y
))
3628 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3629 else if (SCM_BIGP (y
))
3632 if (isnan (SCM_REAL_VALUE (x
)))
3634 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3635 scm_remember_upto_here_1 (y
);
3636 return scm_from_bool (cmp
> 0);
3638 else if (SCM_REALP (y
))
3639 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3640 else if (SCM_FRACTIONP (y
))
3642 double xx
= SCM_REAL_VALUE (x
);
3646 return scm_from_bool (xx
< 0.0);
3647 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3651 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3653 else if (SCM_FRACTIONP (x
))
3655 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3657 /* "a/b < y" becomes "a < y*b" */
3658 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3659 x
= SCM_FRACTION_NUMERATOR (x
);
3662 else if (SCM_REALP (y
))
3664 double yy
= SCM_REAL_VALUE (y
);
3668 return scm_from_bool (0.0 < yy
);
3669 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3672 else if (SCM_FRACTIONP (y
))
3674 /* "a/b < c/d" becomes "a*d < c*b" */
3675 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3676 SCM_FRACTION_DENOMINATOR (y
));
3677 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3678 SCM_FRACTION_DENOMINATOR (x
));
3684 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3687 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3691 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3692 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3693 (SCM x
, SCM y
, SCM rest
),
3694 "Return @code{#t} if the list of parameters is monotonically\n"
3696 #define FUNC_NAME s_scm_i_num_gr_p
3698 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3700 while (!scm_is_null (rest
))
3702 if (scm_is_false (scm_gr_p (x
, y
)))
3706 rest
= scm_cdr (rest
);
3708 return scm_gr_p (x
, y
);
3711 #define FUNC_NAME s_scm_i_num_gr_p
3713 scm_gr_p (SCM x
, SCM y
)
3715 if (!SCM_NUMBERP (x
))
3716 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3717 else if (!SCM_NUMBERP (y
))
3718 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3720 return scm_less_p (y
, x
);
3725 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3726 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3727 (SCM x
, SCM y
, SCM rest
),
3728 "Return @code{#t} if the list of parameters is monotonically\n"
3730 #define FUNC_NAME s_scm_i_num_leq_p
3732 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3734 while (!scm_is_null (rest
))
3736 if (scm_is_false (scm_leq_p (x
, y
)))
3740 rest
= scm_cdr (rest
);
3742 return scm_leq_p (x
, y
);
3745 #define FUNC_NAME s_scm_i_num_leq_p
3747 scm_leq_p (SCM x
, SCM y
)
3749 if (!SCM_NUMBERP (x
))
3750 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3751 else if (!SCM_NUMBERP (y
))
3752 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3753 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3756 return scm_not (scm_less_p (y
, x
));
3761 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3762 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3763 (SCM x
, SCM y
, SCM rest
),
3764 "Return @code{#t} if the list of parameters is monotonically\n"
3766 #define FUNC_NAME s_scm_i_num_geq_p
3768 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3770 while (!scm_is_null (rest
))
3772 if (scm_is_false (scm_geq_p (x
, y
)))
3776 rest
= scm_cdr (rest
);
3778 return scm_geq_p (x
, y
);
3781 #define FUNC_NAME s_scm_i_num_geq_p
3783 scm_geq_p (SCM x
, SCM y
)
3785 if (!SCM_NUMBERP (x
))
3786 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3787 else if (!SCM_NUMBERP (y
))
3788 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3789 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3792 return scm_not (scm_less_p (x
, y
));
3797 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3798 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3804 if (SCM_I_INUMP (z
))
3805 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3806 else if (SCM_BIGP (z
))
3808 else if (SCM_REALP (z
))
3809 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3810 else if (SCM_COMPLEXP (z
))
3811 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3812 && SCM_COMPLEX_IMAG (z
) == 0.0);
3813 else if (SCM_FRACTIONP (z
))
3816 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3820 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3821 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3825 scm_positive_p (SCM x
)
3827 if (SCM_I_INUMP (x
))
3828 return scm_from_bool (SCM_I_INUM (x
) > 0);
3829 else if (SCM_BIGP (x
))
3831 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3832 scm_remember_upto_here_1 (x
);
3833 return scm_from_bool (sgn
> 0);
3835 else if (SCM_REALP (x
))
3836 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3837 else if (SCM_FRACTIONP (x
))
3838 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3840 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3844 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3845 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3849 scm_negative_p (SCM x
)
3851 if (SCM_I_INUMP (x
))
3852 return scm_from_bool (SCM_I_INUM (x
) < 0);
3853 else if (SCM_BIGP (x
))
3855 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3856 scm_remember_upto_here_1 (x
);
3857 return scm_from_bool (sgn
< 0);
3859 else if (SCM_REALP (x
))
3860 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3861 else if (SCM_FRACTIONP (x
))
3862 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3864 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3868 /* scm_min and scm_max return an inexact when either argument is inexact, as
3869 required by r5rs. On that basis, for exact/inexact combinations the
3870 exact is converted to inexact to compare and possibly return. This is
3871 unlike scm_less_p above which takes some trouble to preserve all bits in
3872 its test, such trouble is not required for min and max. */
3874 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3875 (SCM x
, SCM y
, SCM rest
),
3876 "Return the maximum of all parameter values.")
3877 #define FUNC_NAME s_scm_i_max
3879 while (!scm_is_null (rest
))
3880 { x
= scm_max (x
, y
);
3882 rest
= scm_cdr (rest
);
3884 return scm_max (x
, y
);
3888 #define s_max s_scm_i_max
3889 #define g_max g_scm_i_max
3892 scm_max (SCM x
, SCM y
)
3897 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3898 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3901 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3904 if (SCM_I_INUMP (x
))
3906 scm_t_inum xx
= SCM_I_INUM (x
);
3907 if (SCM_I_INUMP (y
))
3909 scm_t_inum yy
= SCM_I_INUM (y
);
3910 return (xx
< yy
) ? y
: x
;
3912 else if (SCM_BIGP (y
))
3914 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3915 scm_remember_upto_here_1 (y
);
3916 return (sgn
< 0) ? x
: y
;
3918 else if (SCM_REALP (y
))
3921 /* if y==NaN then ">" is false and we return NaN */
3922 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3924 else if (SCM_FRACTIONP (y
))
3927 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3930 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3932 else if (SCM_BIGP (x
))
3934 if (SCM_I_INUMP (y
))
3936 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3937 scm_remember_upto_here_1 (x
);
3938 return (sgn
< 0) ? y
: x
;
3940 else if (SCM_BIGP (y
))
3942 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3943 scm_remember_upto_here_2 (x
, y
);
3944 return (cmp
> 0) ? x
: y
;
3946 else if (SCM_REALP (y
))
3948 /* if y==NaN then xx>yy is false, so we return the NaN y */
3951 xx
= scm_i_big2dbl (x
);
3952 yy
= SCM_REAL_VALUE (y
);
3953 return (xx
> yy
? scm_from_double (xx
) : y
);
3955 else if (SCM_FRACTIONP (y
))
3960 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3962 else if (SCM_REALP (x
))
3964 if (SCM_I_INUMP (y
))
3966 double z
= SCM_I_INUM (y
);
3967 /* if x==NaN then "<" is false and we return NaN */
3968 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3970 else if (SCM_BIGP (y
))
3975 else if (SCM_REALP (y
))
3977 /* if x==NaN then our explicit check means we return NaN
3978 if y==NaN then ">" is false and we return NaN
3979 calling isnan is unavoidable, since it's the only way to know
3980 which of x or y causes any compares to be false */
3981 double xx
= SCM_REAL_VALUE (x
);
3982 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3984 else if (SCM_FRACTIONP (y
))
3986 double yy
= scm_i_fraction2double (y
);
3987 double xx
= SCM_REAL_VALUE (x
);
3988 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3991 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3993 else if (SCM_FRACTIONP (x
))
3995 if (SCM_I_INUMP (y
))
3999 else if (SCM_BIGP (y
))
4003 else if (SCM_REALP (y
))
4005 double xx
= scm_i_fraction2double (x
);
4006 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4008 else if (SCM_FRACTIONP (y
))
4013 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4016 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4020 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4021 (SCM x
, SCM y
, SCM rest
),
4022 "Return the minimum of all parameter values.")
4023 #define FUNC_NAME s_scm_i_min
4025 while (!scm_is_null (rest
))
4026 { x
= scm_min (x
, y
);
4028 rest
= scm_cdr (rest
);
4030 return scm_min (x
, y
);
4034 #define s_min s_scm_i_min
4035 #define g_min g_scm_i_min
4038 scm_min (SCM x
, SCM y
)
4043 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4044 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4047 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4050 if (SCM_I_INUMP (x
))
4052 scm_t_inum xx
= SCM_I_INUM (x
);
4053 if (SCM_I_INUMP (y
))
4055 scm_t_inum yy
= SCM_I_INUM (y
);
4056 return (xx
< yy
) ? x
: y
;
4058 else if (SCM_BIGP (y
))
4060 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4061 scm_remember_upto_here_1 (y
);
4062 return (sgn
< 0) ? y
: x
;
4064 else if (SCM_REALP (y
))
4067 /* if y==NaN then "<" is false and we return NaN */
4068 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4070 else if (SCM_FRACTIONP (y
))
4073 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4076 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4078 else if (SCM_BIGP (x
))
4080 if (SCM_I_INUMP (y
))
4082 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4083 scm_remember_upto_here_1 (x
);
4084 return (sgn
< 0) ? x
: y
;
4086 else if (SCM_BIGP (y
))
4088 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4089 scm_remember_upto_here_2 (x
, y
);
4090 return (cmp
> 0) ? y
: x
;
4092 else if (SCM_REALP (y
))
4094 /* if y==NaN then xx<yy is false, so we return the NaN y */
4097 xx
= scm_i_big2dbl (x
);
4098 yy
= SCM_REAL_VALUE (y
);
4099 return (xx
< yy
? scm_from_double (xx
) : y
);
4101 else if (SCM_FRACTIONP (y
))
4106 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4108 else if (SCM_REALP (x
))
4110 if (SCM_I_INUMP (y
))
4112 double z
= SCM_I_INUM (y
);
4113 /* if x==NaN then "<" is false and we return NaN */
4114 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4116 else if (SCM_BIGP (y
))
4121 else if (SCM_REALP (y
))
4123 /* if x==NaN then our explicit check means we return NaN
4124 if y==NaN then "<" is false and we return NaN
4125 calling isnan is unavoidable, since it's the only way to know
4126 which of x or y causes any compares to be false */
4127 double xx
= SCM_REAL_VALUE (x
);
4128 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4130 else if (SCM_FRACTIONP (y
))
4132 double yy
= scm_i_fraction2double (y
);
4133 double xx
= SCM_REAL_VALUE (x
);
4134 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4137 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4139 else if (SCM_FRACTIONP (x
))
4141 if (SCM_I_INUMP (y
))
4145 else if (SCM_BIGP (y
))
4149 else if (SCM_REALP (y
))
4151 double xx
= scm_i_fraction2double (x
);
4152 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4154 else if (SCM_FRACTIONP (y
))
4159 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4162 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4166 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4167 (SCM x
, SCM y
, SCM rest
),
4168 "Return the sum of all parameter values. Return 0 if called without\n"
4170 #define FUNC_NAME s_scm_i_sum
4172 while (!scm_is_null (rest
))
4173 { x
= scm_sum (x
, y
);
4175 rest
= scm_cdr (rest
);
4177 return scm_sum (x
, y
);
4181 #define s_sum s_scm_i_sum
4182 #define g_sum g_scm_i_sum
4185 scm_sum (SCM x
, SCM y
)
4187 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4189 if (SCM_NUMBERP (x
)) return x
;
4190 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4191 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4194 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4196 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4198 scm_t_inum xx
= SCM_I_INUM (x
);
4199 scm_t_inum yy
= SCM_I_INUM (y
);
4200 scm_t_inum z
= xx
+ yy
;
4201 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
4203 else if (SCM_BIGP (y
))
4208 else if (SCM_REALP (y
))
4210 scm_t_inum xx
= SCM_I_INUM (x
);
4211 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4213 else if (SCM_COMPLEXP (y
))
4215 scm_t_inum xx
= SCM_I_INUM (x
);
4216 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4217 SCM_COMPLEX_IMAG (y
));
4219 else if (SCM_FRACTIONP (y
))
4220 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4221 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4222 SCM_FRACTION_DENOMINATOR (y
));
4224 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4225 } else if (SCM_BIGP (x
))
4227 if (SCM_I_INUMP (y
))
4232 inum
= SCM_I_INUM (y
);
4235 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4238 SCM result
= scm_i_mkbig ();
4239 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4240 scm_remember_upto_here_1 (x
);
4241 /* we know the result will have to be a bignum */
4244 return scm_i_normbig (result
);
4248 SCM result
= scm_i_mkbig ();
4249 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4250 scm_remember_upto_here_1 (x
);
4251 /* we know the result will have to be a bignum */
4254 return scm_i_normbig (result
);
4257 else if (SCM_BIGP (y
))
4259 SCM result
= scm_i_mkbig ();
4260 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4261 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4262 mpz_add (SCM_I_BIG_MPZ (result
),
4265 scm_remember_upto_here_2 (x
, y
);
4266 /* we know the result will have to be a bignum */
4269 return scm_i_normbig (result
);
4271 else if (SCM_REALP (y
))
4273 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4274 scm_remember_upto_here_1 (x
);
4275 return scm_from_double (result
);
4277 else if (SCM_COMPLEXP (y
))
4279 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4280 + SCM_COMPLEX_REAL (y
));
4281 scm_remember_upto_here_1 (x
);
4282 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4284 else if (SCM_FRACTIONP (y
))
4285 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4286 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4287 SCM_FRACTION_DENOMINATOR (y
));
4289 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4291 else if (SCM_REALP (x
))
4293 if (SCM_I_INUMP (y
))
4294 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4295 else if (SCM_BIGP (y
))
4297 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4298 scm_remember_upto_here_1 (y
);
4299 return scm_from_double (result
);
4301 else if (SCM_REALP (y
))
4302 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4303 else if (SCM_COMPLEXP (y
))
4304 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4305 SCM_COMPLEX_IMAG (y
));
4306 else if (SCM_FRACTIONP (y
))
4307 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4309 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4311 else if (SCM_COMPLEXP (x
))
4313 if (SCM_I_INUMP (y
))
4314 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4315 SCM_COMPLEX_IMAG (x
));
4316 else if (SCM_BIGP (y
))
4318 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4319 + SCM_COMPLEX_REAL (x
));
4320 scm_remember_upto_here_1 (y
);
4321 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4323 else if (SCM_REALP (y
))
4324 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4325 SCM_COMPLEX_IMAG (x
));
4326 else if (SCM_COMPLEXP (y
))
4327 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4328 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4329 else if (SCM_FRACTIONP (y
))
4330 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4331 SCM_COMPLEX_IMAG (x
));
4333 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4335 else if (SCM_FRACTIONP (x
))
4337 if (SCM_I_INUMP (y
))
4338 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4339 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4340 SCM_FRACTION_DENOMINATOR (x
));
4341 else if (SCM_BIGP (y
))
4342 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4343 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4344 SCM_FRACTION_DENOMINATOR (x
));
4345 else if (SCM_REALP (y
))
4346 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4347 else if (SCM_COMPLEXP (y
))
4348 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4349 SCM_COMPLEX_IMAG (y
));
4350 else if (SCM_FRACTIONP (y
))
4351 /* a/b + c/d = (ad + bc) / bd */
4352 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4353 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4354 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4356 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4359 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4363 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4365 "Return @math{@var{x}+1}.")
4366 #define FUNC_NAME s_scm_oneplus
4368 return scm_sum (x
, SCM_INUM1
);
4373 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4374 (SCM x
, SCM y
, SCM rest
),
4375 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4376 "the sum of all but the first argument are subtracted from the first\n"
4378 #define FUNC_NAME s_scm_i_difference
4380 while (!scm_is_null (rest
))
4381 { x
= scm_difference (x
, y
);
4383 rest
= scm_cdr (rest
);
4385 return scm_difference (x
, y
);
4389 #define s_difference s_scm_i_difference
4390 #define g_difference g_scm_i_difference
4393 scm_difference (SCM x
, SCM y
)
4394 #define FUNC_NAME s_difference
4396 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4399 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4401 if (SCM_I_INUMP (x
))
4403 scm_t_inum xx
= -SCM_I_INUM (x
);
4404 if (SCM_FIXABLE (xx
))
4405 return SCM_I_MAKINUM (xx
);
4407 return scm_i_inum2big (xx
);
4409 else if (SCM_BIGP (x
))
4410 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4411 bignum, but negating that gives a fixnum. */
4412 return scm_i_normbig (scm_i_clonebig (x
, 0));
4413 else if (SCM_REALP (x
))
4414 return scm_from_double (-SCM_REAL_VALUE (x
));
4415 else if (SCM_COMPLEXP (x
))
4416 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4417 -SCM_COMPLEX_IMAG (x
));
4418 else if (SCM_FRACTIONP (x
))
4419 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4420 SCM_FRACTION_DENOMINATOR (x
));
4422 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4425 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4427 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4429 scm_t_inum xx
= SCM_I_INUM (x
);
4430 scm_t_inum yy
= SCM_I_INUM (y
);
4431 scm_t_inum z
= xx
- yy
;
4432 if (SCM_FIXABLE (z
))
4433 return SCM_I_MAKINUM (z
);
4435 return scm_i_inum2big (z
);
4437 else if (SCM_BIGP (y
))
4439 /* inum-x - big-y */
4440 scm_t_inum xx
= SCM_I_INUM (x
);
4444 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4445 bignum, but negating that gives a fixnum. */
4446 return scm_i_normbig (scm_i_clonebig (y
, 0));
4450 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4451 SCM result
= scm_i_mkbig ();
4454 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4457 /* x - y == -(y + -x) */
4458 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4459 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4461 scm_remember_upto_here_1 (y
);
4463 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4464 /* we know the result will have to be a bignum */
4467 return scm_i_normbig (result
);
4470 else if (SCM_REALP (y
))
4472 scm_t_inum xx
= SCM_I_INUM (x
);
4473 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4475 else if (SCM_COMPLEXP (y
))
4477 scm_t_inum xx
= SCM_I_INUM (x
);
4478 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4479 - SCM_COMPLEX_IMAG (y
));
4481 else if (SCM_FRACTIONP (y
))
4482 /* a - b/c = (ac - b) / c */
4483 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4484 SCM_FRACTION_NUMERATOR (y
)),
4485 SCM_FRACTION_DENOMINATOR (y
));
4487 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4489 else if (SCM_BIGP (x
))
4491 if (SCM_I_INUMP (y
))
4493 /* big-x - inum-y */
4494 scm_t_inum yy
= SCM_I_INUM (y
);
4495 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4497 scm_remember_upto_here_1 (x
);
4499 return (SCM_FIXABLE (-yy
) ?
4500 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
4503 SCM result
= scm_i_mkbig ();
4506 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4508 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4509 scm_remember_upto_here_1 (x
);
4511 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4512 /* we know the result will have to be a bignum */
4515 return scm_i_normbig (result
);
4518 else if (SCM_BIGP (y
))
4520 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4521 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4522 SCM result
= scm_i_mkbig ();
4523 mpz_sub (SCM_I_BIG_MPZ (result
),
4526 scm_remember_upto_here_2 (x
, y
);
4527 /* we know the result will have to be a bignum */
4528 if ((sgn_x
== 1) && (sgn_y
== -1))
4530 if ((sgn_x
== -1) && (sgn_y
== 1))
4532 return scm_i_normbig (result
);
4534 else if (SCM_REALP (y
))
4536 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4537 scm_remember_upto_here_1 (x
);
4538 return scm_from_double (result
);
4540 else if (SCM_COMPLEXP (y
))
4542 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4543 - SCM_COMPLEX_REAL (y
));
4544 scm_remember_upto_here_1 (x
);
4545 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4547 else if (SCM_FRACTIONP (y
))
4548 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4549 SCM_FRACTION_NUMERATOR (y
)),
4550 SCM_FRACTION_DENOMINATOR (y
));
4551 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4553 else if (SCM_REALP (x
))
4555 if (SCM_I_INUMP (y
))
4556 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4557 else if (SCM_BIGP (y
))
4559 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4560 scm_remember_upto_here_1 (x
);
4561 return scm_from_double (result
);
4563 else if (SCM_REALP (y
))
4564 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4565 else if (SCM_COMPLEXP (y
))
4566 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4567 -SCM_COMPLEX_IMAG (y
));
4568 else if (SCM_FRACTIONP (y
))
4569 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4571 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4573 else if (SCM_COMPLEXP (x
))
4575 if (SCM_I_INUMP (y
))
4576 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4577 SCM_COMPLEX_IMAG (x
));
4578 else if (SCM_BIGP (y
))
4580 double real_part
= (SCM_COMPLEX_REAL (x
)
4581 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4582 scm_remember_upto_here_1 (x
);
4583 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4585 else if (SCM_REALP (y
))
4586 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4587 SCM_COMPLEX_IMAG (x
));
4588 else if (SCM_COMPLEXP (y
))
4589 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4590 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4591 else if (SCM_FRACTIONP (y
))
4592 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4593 SCM_COMPLEX_IMAG (x
));
4595 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4597 else if (SCM_FRACTIONP (x
))
4599 if (SCM_I_INUMP (y
))
4600 /* a/b - c = (a - cb) / b */
4601 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4602 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4603 SCM_FRACTION_DENOMINATOR (x
));
4604 else if (SCM_BIGP (y
))
4605 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4606 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4607 SCM_FRACTION_DENOMINATOR (x
));
4608 else if (SCM_REALP (y
))
4609 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4610 else if (SCM_COMPLEXP (y
))
4611 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4612 -SCM_COMPLEX_IMAG (y
));
4613 else if (SCM_FRACTIONP (y
))
4614 /* a/b - c/d = (ad - bc) / bd */
4615 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4616 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4617 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4619 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4622 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4627 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4629 "Return @math{@var{x}-1}.")
4630 #define FUNC_NAME s_scm_oneminus
4632 return scm_difference (x
, SCM_INUM1
);
4637 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4638 (SCM x
, SCM y
, SCM rest
),
4639 "Return the product of all arguments. If called without arguments,\n"
4641 #define FUNC_NAME s_scm_i_product
4643 while (!scm_is_null (rest
))
4644 { x
= scm_product (x
, y
);
4646 rest
= scm_cdr (rest
);
4648 return scm_product (x
, y
);
4652 #define s_product s_scm_i_product
4653 #define g_product g_scm_i_product
4656 scm_product (SCM x
, SCM y
)
4658 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4661 return SCM_I_MAKINUM (1L);
4662 else if (SCM_NUMBERP (x
))
4665 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4668 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4673 xx
= SCM_I_INUM (x
);
4677 case 0: return x
; break;
4678 case 1: return y
; break;
4680 * The following case (x = -1) is important for more than
4681 * just optimization. It handles the case of negating
4682 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
4683 * which is a bignum that must be changed back into a fixnum.
4684 * Failure to do so will cause the following to return #f:
4685 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
4688 return scm_difference(y
, SCM_UNDEFINED
);
4692 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4694 scm_t_inum yy
= SCM_I_INUM (y
);
4695 scm_t_inum kk
= xx
* yy
;
4696 SCM k
= SCM_I_MAKINUM (kk
);
4697 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4701 SCM result
= scm_i_inum2big (xx
);
4702 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4703 return scm_i_normbig (result
);
4706 else if (SCM_BIGP (y
))
4708 SCM result
= scm_i_mkbig ();
4709 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4710 scm_remember_upto_here_1 (y
);
4713 else if (SCM_REALP (y
))
4714 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4715 else if (SCM_COMPLEXP (y
))
4716 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4717 xx
* SCM_COMPLEX_IMAG (y
));
4718 else if (SCM_FRACTIONP (y
))
4719 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4720 SCM_FRACTION_DENOMINATOR (y
));
4722 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4724 else if (SCM_BIGP (x
))
4726 if (SCM_I_INUMP (y
))
4731 else if (SCM_BIGP (y
))
4733 SCM result
= scm_i_mkbig ();
4734 mpz_mul (SCM_I_BIG_MPZ (result
),
4737 scm_remember_upto_here_2 (x
, y
);
4740 else if (SCM_REALP (y
))
4742 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4743 scm_remember_upto_here_1 (x
);
4744 return scm_from_double (result
);
4746 else if (SCM_COMPLEXP (y
))
4748 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4749 scm_remember_upto_here_1 (x
);
4750 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4751 z
* SCM_COMPLEX_IMAG (y
));
4753 else if (SCM_FRACTIONP (y
))
4754 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4755 SCM_FRACTION_DENOMINATOR (y
));
4757 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4759 else if (SCM_REALP (x
))
4761 if (SCM_I_INUMP (y
))
4763 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4764 if (scm_is_eq (y
, SCM_INUM0
))
4766 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4768 else if (SCM_BIGP (y
))
4770 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4771 scm_remember_upto_here_1 (y
);
4772 return scm_from_double (result
);
4774 else if (SCM_REALP (y
))
4775 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4776 else if (SCM_COMPLEXP (y
))
4777 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4778 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4779 else if (SCM_FRACTIONP (y
))
4780 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4782 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4784 else if (SCM_COMPLEXP (x
))
4786 if (SCM_I_INUMP (y
))
4788 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4789 if (scm_is_eq (y
, SCM_INUM0
))
4791 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4792 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4794 else if (SCM_BIGP (y
))
4796 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4797 scm_remember_upto_here_1 (y
);
4798 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4799 z
* SCM_COMPLEX_IMAG (x
));
4801 else if (SCM_REALP (y
))
4802 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4803 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4804 else if (SCM_COMPLEXP (y
))
4806 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4807 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4808 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4809 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4811 else if (SCM_FRACTIONP (y
))
4813 double yy
= scm_i_fraction2double (y
);
4814 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4815 yy
* SCM_COMPLEX_IMAG (x
));
4818 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4820 else if (SCM_FRACTIONP (x
))
4822 if (SCM_I_INUMP (y
))
4823 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4824 SCM_FRACTION_DENOMINATOR (x
));
4825 else if (SCM_BIGP (y
))
4826 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4827 SCM_FRACTION_DENOMINATOR (x
));
4828 else if (SCM_REALP (y
))
4829 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4830 else if (SCM_COMPLEXP (y
))
4832 double xx
= scm_i_fraction2double (x
);
4833 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4834 xx
* SCM_COMPLEX_IMAG (y
));
4836 else if (SCM_FRACTIONP (y
))
4837 /* a/b * c/d = ac / bd */
4838 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4839 SCM_FRACTION_NUMERATOR (y
)),
4840 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4841 SCM_FRACTION_DENOMINATOR (y
)));
4843 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4846 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4849 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4850 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4851 #define ALLOW_DIVIDE_BY_ZERO
4852 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4855 /* The code below for complex division is adapted from the GNU
4856 libstdc++, which adapted it from f2c's libF77, and is subject to
4859 /****************************************************************
4860 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4862 Permission to use, copy, modify, and distribute this software
4863 and its documentation for any purpose and without fee is hereby
4864 granted, provided that the above copyright notice appear in all
4865 copies and that both that the copyright notice and this
4866 permission notice and warranty disclaimer appear in supporting
4867 documentation, and that the names of AT&T Bell Laboratories or
4868 Bellcore or any of their entities not be used in advertising or
4869 publicity pertaining to distribution of the software without
4870 specific, written prior permission.
4872 AT&T and Bellcore disclaim all warranties with regard to this
4873 software, including all implied warranties of merchantability
4874 and fitness. In no event shall AT&T or Bellcore be liable for
4875 any special, indirect or consequential damages or any damages
4876 whatsoever resulting from loss of use, data or profits, whether
4877 in an action of contract, negligence or other tortious action,
4878 arising out of or in connection with the use or performance of
4880 ****************************************************************/
4882 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4883 (SCM x
, SCM y
, SCM rest
),
4884 "Divide the first argument by the product of the remaining\n"
4885 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4887 #define FUNC_NAME s_scm_i_divide
4889 while (!scm_is_null (rest
))
4890 { x
= scm_divide (x
, y
);
4892 rest
= scm_cdr (rest
);
4894 return scm_divide (x
, y
);
4898 #define s_divide s_scm_i_divide
4899 #define g_divide g_scm_i_divide
4902 do_divide (SCM x
, SCM y
, int inexact
)
4903 #define FUNC_NAME s_divide
4907 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4910 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4911 else if (SCM_I_INUMP (x
))
4913 scm_t_inum xx
= SCM_I_INUM (x
);
4914 if (xx
== 1 || xx
== -1)
4916 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4918 scm_num_overflow (s_divide
);
4923 return scm_from_double (1.0 / (double) xx
);
4924 else return scm_i_make_ratio (SCM_INUM1
, x
);
4927 else if (SCM_BIGP (x
))
4930 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4931 else return scm_i_make_ratio (SCM_INUM1
, x
);
4933 else if (SCM_REALP (x
))
4935 double xx
= SCM_REAL_VALUE (x
);
4936 #ifndef ALLOW_DIVIDE_BY_ZERO
4938 scm_num_overflow (s_divide
);
4941 return scm_from_double (1.0 / xx
);
4943 else if (SCM_COMPLEXP (x
))
4945 double r
= SCM_COMPLEX_REAL (x
);
4946 double i
= SCM_COMPLEX_IMAG (x
);
4947 if (fabs(r
) <= fabs(i
))
4950 double d
= i
* (1.0 + t
* t
);
4951 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4956 double d
= r
* (1.0 + t
* t
);
4957 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4960 else if (SCM_FRACTIONP (x
))
4961 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4962 SCM_FRACTION_NUMERATOR (x
));
4964 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4967 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4969 scm_t_inum xx
= SCM_I_INUM (x
);
4970 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4972 scm_t_inum yy
= SCM_I_INUM (y
);
4975 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4976 scm_num_overflow (s_divide
);
4978 return scm_from_double ((double) xx
/ (double) yy
);
4981 else if (xx
% yy
!= 0)
4984 return scm_from_double ((double) xx
/ (double) yy
);
4985 else return scm_i_make_ratio (x
, y
);
4989 scm_t_inum z
= xx
/ yy
;
4990 if (SCM_FIXABLE (z
))
4991 return SCM_I_MAKINUM (z
);
4993 return scm_i_inum2big (z
);
4996 else if (SCM_BIGP (y
))
4999 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
5000 else return scm_i_make_ratio (x
, y
);
5002 else if (SCM_REALP (y
))
5004 double yy
= SCM_REAL_VALUE (y
);
5005 #ifndef ALLOW_DIVIDE_BY_ZERO
5007 scm_num_overflow (s_divide
);
5010 return scm_from_double ((double) xx
/ yy
);
5012 else if (SCM_COMPLEXP (y
))
5015 complex_div
: /* y _must_ be a complex number */
5017 double r
= SCM_COMPLEX_REAL (y
);
5018 double i
= SCM_COMPLEX_IMAG (y
);
5019 if (fabs(r
) <= fabs(i
))
5022 double d
= i
* (1.0 + t
* t
);
5023 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5028 double d
= r
* (1.0 + t
* t
);
5029 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5033 else if (SCM_FRACTIONP (y
))
5034 /* a / b/c = ac / b */
5035 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5036 SCM_FRACTION_NUMERATOR (y
));
5038 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5040 else if (SCM_BIGP (x
))
5042 if (SCM_I_INUMP (y
))
5044 scm_t_inum yy
= SCM_I_INUM (y
);
5047 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5048 scm_num_overflow (s_divide
);
5050 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5051 scm_remember_upto_here_1 (x
);
5052 return (sgn
== 0) ? scm_nan () : scm_inf ();
5059 /* FIXME: HMM, what are the relative performance issues here?
5060 We need to test. Is it faster on average to test
5061 divisible_p, then perform whichever operation, or is it
5062 faster to perform the integer div opportunistically and
5063 switch to real if there's a remainder? For now we take the
5064 middle ground: test, then if divisible, use the faster div
5067 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
5068 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5072 SCM result
= scm_i_mkbig ();
5073 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5074 scm_remember_upto_here_1 (x
);
5076 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5077 return scm_i_normbig (result
);
5082 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5083 else return scm_i_make_ratio (x
, y
);
5087 else if (SCM_BIGP (y
))
5092 /* It's easily possible for the ratio x/y to fit a double
5093 but one or both x and y be too big to fit a double,
5094 hence the use of mpq_get_d rather than converting and
5097 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5098 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5099 return scm_from_double (mpq_get_d (q
));
5103 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5107 SCM result
= scm_i_mkbig ();
5108 mpz_divexact (SCM_I_BIG_MPZ (result
),
5111 scm_remember_upto_here_2 (x
, y
);
5112 return scm_i_normbig (result
);
5115 return scm_i_make_ratio (x
, y
);
5118 else if (SCM_REALP (y
))
5120 double yy
= SCM_REAL_VALUE (y
);
5121 #ifndef ALLOW_DIVIDE_BY_ZERO
5123 scm_num_overflow (s_divide
);
5126 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5128 else if (SCM_COMPLEXP (y
))
5130 a
= scm_i_big2dbl (x
);
5133 else if (SCM_FRACTIONP (y
))
5134 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5135 SCM_FRACTION_NUMERATOR (y
));
5137 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5139 else if (SCM_REALP (x
))
5141 double rx
= SCM_REAL_VALUE (x
);
5142 if (SCM_I_INUMP (y
))
5144 scm_t_inum yy
= SCM_I_INUM (y
);
5145 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5147 scm_num_overflow (s_divide
);
5150 return scm_from_double (rx
/ (double) yy
);
5152 else if (SCM_BIGP (y
))
5154 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5155 scm_remember_upto_here_1 (y
);
5156 return scm_from_double (rx
/ dby
);
5158 else if (SCM_REALP (y
))
5160 double yy
= SCM_REAL_VALUE (y
);
5161 #ifndef ALLOW_DIVIDE_BY_ZERO
5163 scm_num_overflow (s_divide
);
5166 return scm_from_double (rx
/ yy
);
5168 else if (SCM_COMPLEXP (y
))
5173 else if (SCM_FRACTIONP (y
))
5174 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5176 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5178 else if (SCM_COMPLEXP (x
))
5180 double rx
= SCM_COMPLEX_REAL (x
);
5181 double ix
= SCM_COMPLEX_IMAG (x
);
5182 if (SCM_I_INUMP (y
))
5184 scm_t_inum yy
= SCM_I_INUM (y
);
5185 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5187 scm_num_overflow (s_divide
);
5192 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5195 else if (SCM_BIGP (y
))
5197 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5198 scm_remember_upto_here_1 (y
);
5199 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5201 else if (SCM_REALP (y
))
5203 double yy
= SCM_REAL_VALUE (y
);
5204 #ifndef ALLOW_DIVIDE_BY_ZERO
5206 scm_num_overflow (s_divide
);
5209 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5211 else if (SCM_COMPLEXP (y
))
5213 double ry
= SCM_COMPLEX_REAL (y
);
5214 double iy
= SCM_COMPLEX_IMAG (y
);
5215 if (fabs(ry
) <= fabs(iy
))
5218 double d
= iy
* (1.0 + t
* t
);
5219 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5224 double d
= ry
* (1.0 + t
* t
);
5225 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5228 else if (SCM_FRACTIONP (y
))
5230 double yy
= scm_i_fraction2double (y
);
5231 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5234 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5236 else if (SCM_FRACTIONP (x
))
5238 if (SCM_I_INUMP (y
))
5240 scm_t_inum yy
= SCM_I_INUM (y
);
5241 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5243 scm_num_overflow (s_divide
);
5246 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5247 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5249 else if (SCM_BIGP (y
))
5251 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5252 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5254 else if (SCM_REALP (y
))
5256 double yy
= SCM_REAL_VALUE (y
);
5257 #ifndef ALLOW_DIVIDE_BY_ZERO
5259 scm_num_overflow (s_divide
);
5262 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5264 else if (SCM_COMPLEXP (y
))
5266 a
= scm_i_fraction2double (x
);
5269 else if (SCM_FRACTIONP (y
))
5270 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5271 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5273 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5276 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5280 scm_divide (SCM x
, SCM y
)
5282 return do_divide (x
, y
, 0);
5285 static SCM
scm_divide2real (SCM x
, SCM y
)
5287 return do_divide (x
, y
, 1);
5293 scm_c_truncate (double x
)
5304 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5305 half-way case (ie. when x is an integer plus 0.5) going upwards.
5306 Then half-way cases are identified and adjusted down if the
5307 round-upwards didn't give the desired even integer.
5309 "plus_half == result" identifies a half-way case. If plus_half, which is
5310 x + 0.5, is an integer then x must be an integer plus 0.5.
5312 An odd "result" value is identified with result/2 != floor(result/2).
5313 This is done with plus_half, since that value is ready for use sooner in
5314 a pipelined cpu, and we're already requiring plus_half == result.
5316 Note however that we need to be careful when x is big and already an
5317 integer. In that case "x+0.5" may round to an adjacent integer, causing
5318 us to return such a value, incorrectly. For instance if the hardware is
5319 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5320 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5321 returned. Or if the hardware is in round-upwards mode, then other bigger
5322 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5323 representable value, 2^128+2^76 (or whatever), again incorrect.
5325 These bad roundings of x+0.5 are avoided by testing at the start whether
5326 x is already an integer. If it is then clearly that's the desired result
5327 already. And if it's not then the exponent must be small enough to allow
5328 an 0.5 to be represented, and hence added without a bad rounding. */
5331 scm_c_round (double x
)
5333 double plus_half
, result
;
5338 plus_half
= x
+ 0.5;
5339 result
= floor (plus_half
);
5340 /* Adjust so that the rounding is towards even. */
5341 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5346 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5348 "Round the number @var{x} towards zero.")
5349 #define FUNC_NAME s_scm_truncate_number
5351 if (scm_is_false (scm_negative_p (x
)))
5352 return scm_floor (x
);
5354 return scm_ceiling (x
);
5358 static SCM exactly_one_half
;
5360 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5362 "Round the number @var{x} towards the nearest integer. "
5363 "When it is exactly halfway between two integers, "
5364 "round towards the even one.")
5365 #define FUNC_NAME s_scm_round_number
5367 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5369 else if (SCM_REALP (x
))
5370 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5373 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5374 single quotient+remainder division then examining to see which way
5375 the rounding should go. */
5376 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5377 SCM result
= scm_floor (plus_half
);
5378 /* Adjust so that the rounding is towards even. */
5379 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5380 && scm_is_true (scm_odd_p (result
)))
5381 return scm_difference (result
, SCM_INUM1
);
5388 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5390 "Round the number @var{x} towards minus infinity.")
5391 #define FUNC_NAME s_scm_floor
5393 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5395 else if (SCM_REALP (x
))
5396 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5397 else if (SCM_FRACTIONP (x
))
5399 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5400 SCM_FRACTION_DENOMINATOR (x
));
5401 if (scm_is_false (scm_negative_p (x
)))
5403 /* For positive x, rounding towards zero is correct. */
5408 /* For negative x, we need to return q-1 unless x is an
5409 integer. But fractions are never integer, per our
5411 return scm_difference (q
, SCM_INUM1
);
5415 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5419 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5421 "Round the number @var{x} towards infinity.")
5422 #define FUNC_NAME s_scm_ceiling
5424 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5426 else if (SCM_REALP (x
))
5427 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5428 else if (SCM_FRACTIONP (x
))
5430 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5431 SCM_FRACTION_DENOMINATOR (x
));
5432 if (scm_is_false (scm_positive_p (x
)))
5434 /* For negative x, rounding towards zero is correct. */
5439 /* For positive x, we need to return q+1 unless x is an
5440 integer. But fractions are never integer, per our
5442 return scm_sum (q
, SCM_INUM1
);
5446 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5450 /* sin/cos/tan/asin/acos/atan
5451 sinh/cosh/tanh/asinh/acosh/atanh
5452 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5453 Written by Jerry D. Hedden, (C) FSF.
5454 See the file `COPYING' for terms applying to this program. */
5456 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5458 "Return @var{x} raised to the power of @var{y}.")
5459 #define FUNC_NAME s_scm_expt
5461 if (scm_is_integer (y
))
5463 if (scm_is_true (scm_exact_p (y
)))
5464 return scm_integer_expt (x
, y
);
5467 /* Here we handle the case where the exponent is an inexact
5468 integer. We make the exponent exact in order to use
5469 scm_integer_expt, and thus avoid the spurious imaginary
5470 parts that may result from round-off errors in the general
5471 e^(y log x) method below (for example when squaring a large
5472 negative number). In this case, we must return an inexact
5473 result for correctness. We also make the base inexact so
5474 that scm_integer_expt will use fast inexact arithmetic
5475 internally. Note that making the base inexact is not
5476 sufficient to guarantee an inexact result, because
5477 scm_integer_expt will return an exact 1 when the exponent
5478 is 0, even if the base is inexact. */
5479 return scm_exact_to_inexact
5480 (scm_integer_expt (scm_exact_to_inexact (x
),
5481 scm_inexact_to_exact (y
)));
5484 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5486 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5489 return scm_exp (scm_product (scm_log (x
), y
));
5493 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5495 "Compute the sine of @var{z}.")
5496 #define FUNC_NAME s_scm_sin
5498 if (scm_is_real (z
))
5499 return scm_from_double (sin (scm_to_double (z
)));
5500 else if (SCM_COMPLEXP (z
))
5502 x
= SCM_COMPLEX_REAL (z
);
5503 y
= SCM_COMPLEX_IMAG (z
);
5504 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5505 cos (x
) * sinh (y
));
5508 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5512 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5514 "Compute the cosine of @var{z}.")
5515 #define FUNC_NAME s_scm_cos
5517 if (scm_is_real (z
))
5518 return scm_from_double (cos (scm_to_double (z
)));
5519 else if (SCM_COMPLEXP (z
))
5521 x
= SCM_COMPLEX_REAL (z
);
5522 y
= SCM_COMPLEX_IMAG (z
);
5523 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5524 -sin (x
) * sinh (y
));
5527 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5531 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5533 "Compute the tangent of @var{z}.")
5534 #define FUNC_NAME s_scm_tan
5536 if (scm_is_real (z
))
5537 return scm_from_double (tan (scm_to_double (z
)));
5538 else if (SCM_COMPLEXP (z
))
5540 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5541 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5542 w
= cos (x
) + cosh (y
);
5543 #ifndef ALLOW_DIVIDE_BY_ZERO
5545 scm_num_overflow (s_scm_tan
);
5547 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5550 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5554 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5556 "Compute the hyperbolic sine of @var{z}.")
5557 #define FUNC_NAME s_scm_sinh
5559 if (scm_is_real (z
))
5560 return scm_from_double (sinh (scm_to_double (z
)));
5561 else if (SCM_COMPLEXP (z
))
5563 x
= SCM_COMPLEX_REAL (z
);
5564 y
= SCM_COMPLEX_IMAG (z
);
5565 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5566 cosh (x
) * sin (y
));
5569 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5573 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5575 "Compute the hyperbolic cosine of @var{z}.")
5576 #define FUNC_NAME s_scm_cosh
5578 if (scm_is_real (z
))
5579 return scm_from_double (cosh (scm_to_double (z
)));
5580 else if (SCM_COMPLEXP (z
))
5582 x
= SCM_COMPLEX_REAL (z
);
5583 y
= SCM_COMPLEX_IMAG (z
);
5584 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5585 sinh (x
) * sin (y
));
5588 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5592 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5594 "Compute the hyperbolic tangent of @var{z}.")
5595 #define FUNC_NAME s_scm_tanh
5597 if (scm_is_real (z
))
5598 return scm_from_double (tanh (scm_to_double (z
)));
5599 else if (SCM_COMPLEXP (z
))
5601 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5602 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5603 w
= cosh (x
) + cos (y
);
5604 #ifndef ALLOW_DIVIDE_BY_ZERO
5606 scm_num_overflow (s_scm_tanh
);
5608 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5611 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5615 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5617 "Compute the arc sine of @var{z}.")
5618 #define FUNC_NAME s_scm_asin
5620 if (scm_is_real (z
))
5622 double w
= scm_to_double (z
);
5623 if (w
>= -1.0 && w
<= 1.0)
5624 return scm_from_double (asin (w
));
5626 return scm_product (scm_c_make_rectangular (0, -1),
5627 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5629 else if (SCM_COMPLEXP (z
))
5631 x
= SCM_COMPLEX_REAL (z
);
5632 y
= SCM_COMPLEX_IMAG (z
);
5633 return scm_product (scm_c_make_rectangular (0, -1),
5634 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5637 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5641 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5643 "Compute the arc cosine of @var{z}.")
5644 #define FUNC_NAME s_scm_acos
5646 if (scm_is_real (z
))
5648 double w
= scm_to_double (z
);
5649 if (w
>= -1.0 && w
<= 1.0)
5650 return scm_from_double (acos (w
));
5652 return scm_sum (scm_from_double (acos (0.0)),
5653 scm_product (scm_c_make_rectangular (0, 1),
5654 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5656 else if (SCM_COMPLEXP (z
))
5658 x
= SCM_COMPLEX_REAL (z
);
5659 y
= SCM_COMPLEX_IMAG (z
);
5660 return scm_sum (scm_from_double (acos (0.0)),
5661 scm_product (scm_c_make_rectangular (0, 1),
5662 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5665 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5669 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5671 "With one argument, compute the arc tangent of @var{z}.\n"
5672 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5673 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5674 #define FUNC_NAME s_scm_atan
5678 if (scm_is_real (z
))
5679 return scm_from_double (atan (scm_to_double (z
)));
5680 else if (SCM_COMPLEXP (z
))
5683 v
= SCM_COMPLEX_REAL (z
);
5684 w
= SCM_COMPLEX_IMAG (z
);
5685 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5686 scm_c_make_rectangular (v
, w
+ 1.0))),
5687 scm_c_make_rectangular (0, 2));
5690 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5692 else if (scm_is_real (z
))
5694 if (scm_is_real (y
))
5695 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5697 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5700 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5704 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5706 "Compute the inverse hyperbolic sine of @var{z}.")
5707 #define FUNC_NAME s_scm_sys_asinh
5709 if (scm_is_real (z
))
5710 return scm_from_double (asinh (scm_to_double (z
)));
5711 else if (scm_is_number (z
))
5712 return scm_log (scm_sum (z
,
5713 scm_sqrt (scm_sum (scm_product (z
, z
),
5716 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5720 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5722 "Compute the inverse hyperbolic cosine of @var{z}.")
5723 #define FUNC_NAME s_scm_sys_acosh
5725 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5726 return scm_from_double (acosh (scm_to_double (z
)));
5727 else if (scm_is_number (z
))
5728 return scm_log (scm_sum (z
,
5729 scm_sqrt (scm_difference (scm_product (z
, z
),
5732 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5736 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5738 "Compute the inverse hyperbolic tangent of @var{z}.")
5739 #define FUNC_NAME s_scm_sys_atanh
5741 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5742 return scm_from_double (atanh (scm_to_double (z
)));
5743 else if (scm_is_number (z
))
5744 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
5745 scm_difference (SCM_INUM1
, z
))),
5748 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5753 scm_c_make_rectangular (double re
, double im
)
5756 return scm_from_double (re
);
5761 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5763 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
5764 SCM_COMPLEX_REAL (z
) = re
;
5765 SCM_COMPLEX_IMAG (z
) = im
;
5770 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5771 (SCM real_part
, SCM imaginary_part
),
5772 "Return a complex number constructed of the given @var{real-part} "
5773 "and @var{imaginary-part} parts.")
5774 #define FUNC_NAME s_scm_make_rectangular
5776 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5777 SCM_ARG1
, FUNC_NAME
, "real");
5778 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5779 SCM_ARG2
, FUNC_NAME
, "real");
5780 return scm_c_make_rectangular (scm_to_double (real_part
),
5781 scm_to_double (imaginary_part
));
5786 scm_c_make_polar (double mag
, double ang
)
5790 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5791 use it on Glibc-based systems that have it (it's a GNU extension). See
5792 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5794 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5795 sincos (ang
, &s
, &c
);
5800 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5803 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5805 "Return the complex number @var{x} * e^(i * @var{y}).")
5806 #define FUNC_NAME s_scm_make_polar
5808 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5809 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5810 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5815 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5816 /* "Return the real part of the number @var{z}."
5819 scm_real_part (SCM z
)
5821 if (SCM_I_INUMP (z
))
5823 else if (SCM_BIGP (z
))
5825 else if (SCM_REALP (z
))
5827 else if (SCM_COMPLEXP (z
))
5828 return scm_from_double (SCM_COMPLEX_REAL (z
));
5829 else if (SCM_FRACTIONP (z
))
5832 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5836 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5837 /* "Return the imaginary part of the number @var{z}."
5840 scm_imag_part (SCM z
)
5842 if (SCM_I_INUMP (z
))
5844 else if (SCM_BIGP (z
))
5846 else if (SCM_REALP (z
))
5848 else if (SCM_COMPLEXP (z
))
5849 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5850 else if (SCM_FRACTIONP (z
))
5853 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5856 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5857 /* "Return the numerator of the number @var{z}."
5860 scm_numerator (SCM z
)
5862 if (SCM_I_INUMP (z
))
5864 else if (SCM_BIGP (z
))
5866 else if (SCM_FRACTIONP (z
))
5867 return SCM_FRACTION_NUMERATOR (z
);
5868 else if (SCM_REALP (z
))
5869 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5871 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5875 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5876 /* "Return the denominator of the number @var{z}."
5879 scm_denominator (SCM z
)
5881 if (SCM_I_INUMP (z
))
5883 else if (SCM_BIGP (z
))
5885 else if (SCM_FRACTIONP (z
))
5886 return SCM_FRACTION_DENOMINATOR (z
);
5887 else if (SCM_REALP (z
))
5888 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5890 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5893 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5894 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5895 * "@code{abs} for real arguments, but also allows complex numbers."
5898 scm_magnitude (SCM z
)
5900 if (SCM_I_INUMP (z
))
5902 scm_t_inum zz
= SCM_I_INUM (z
);
5905 else if (SCM_POSFIXABLE (-zz
))
5906 return SCM_I_MAKINUM (-zz
);
5908 return scm_i_inum2big (-zz
);
5910 else if (SCM_BIGP (z
))
5912 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5913 scm_remember_upto_here_1 (z
);
5915 return scm_i_clonebig (z
, 0);
5919 else if (SCM_REALP (z
))
5920 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5921 else if (SCM_COMPLEXP (z
))
5922 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5923 else if (SCM_FRACTIONP (z
))
5925 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5927 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5928 SCM_FRACTION_DENOMINATOR (z
));
5931 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5935 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5936 /* "Return the angle of the complex number @var{z}."
5941 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5942 flo0 to save allocating a new flonum with scm_from_double each time.
5943 But if atan2 follows the floating point rounding mode, then the value
5944 is not a constant. Maybe it'd be close enough though. */
5945 if (SCM_I_INUMP (z
))
5947 if (SCM_I_INUM (z
) >= 0)
5950 return scm_from_double (atan2 (0.0, -1.0));
5952 else if (SCM_BIGP (z
))
5954 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5955 scm_remember_upto_here_1 (z
);
5957 return scm_from_double (atan2 (0.0, -1.0));
5961 else if (SCM_REALP (z
))
5963 if (SCM_REAL_VALUE (z
) >= 0)
5966 return scm_from_double (atan2 (0.0, -1.0));
5968 else if (SCM_COMPLEXP (z
))
5969 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5970 else if (SCM_FRACTIONP (z
))
5972 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5974 else return scm_from_double (atan2 (0.0, -1.0));
5977 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5981 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5982 /* Convert the number @var{x} to its inexact representation.\n"
5985 scm_exact_to_inexact (SCM z
)
5987 if (SCM_I_INUMP (z
))
5988 return scm_from_double ((double) SCM_I_INUM (z
));
5989 else if (SCM_BIGP (z
))
5990 return scm_from_double (scm_i_big2dbl (z
));
5991 else if (SCM_FRACTIONP (z
))
5992 return scm_from_double (scm_i_fraction2double (z
));
5993 else if (SCM_INEXACTP (z
))
5996 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
6000 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
6002 "Return an exact number that is numerically closest to @var{z}.")
6003 #define FUNC_NAME s_scm_inexact_to_exact
6005 if (SCM_I_INUMP (z
))
6007 else if (SCM_BIGP (z
))
6009 else if (SCM_REALP (z
))
6011 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
6012 SCM_OUT_OF_RANGE (1, z
);
6019 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6020 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6021 scm_i_mpz2num (mpq_denref (frac
)));
6023 /* When scm_i_make_ratio throws, we leak the memory allocated
6030 else if (SCM_FRACTIONP (z
))
6033 SCM_WRONG_TYPE_ARG (1, z
);
6037 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6039 "Returns the @emph{simplest} rational number differing\n"
6040 "from @var{x} by no more than @var{eps}.\n"
6042 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6043 "exact result when both its arguments are exact. Thus, you might need\n"
6044 "to use @code{inexact->exact} on the arguments.\n"
6047 "(rationalize (inexact->exact 1.2) 1/100)\n"
6050 #define FUNC_NAME s_scm_rationalize
6052 if (SCM_I_INUMP (x
))
6054 else if (SCM_BIGP (x
))
6056 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6058 /* Use continued fractions to find closest ratio. All
6059 arithmetic is done with exact numbers.
6062 SCM ex
= scm_inexact_to_exact (x
);
6063 SCM int_part
= scm_floor (ex
);
6065 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
6066 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
6070 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6073 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6074 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6076 /* We stop after a million iterations just to be absolutely sure
6077 that we don't go into an infinite loop. The process normally
6078 converges after less than a dozen iterations.
6081 eps
= scm_abs (eps
);
6082 while (++i
< 1000000)
6084 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6085 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6086 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6088 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6089 eps
))) /* abs(x-a/b) <= eps */
6091 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6092 if (scm_is_false (scm_exact_p (x
))
6093 || scm_is_false (scm_exact_p (eps
)))
6094 return scm_exact_to_inexact (res
);
6098 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6100 tt
= scm_floor (rx
); /* tt = floor (rx) */
6106 scm_num_overflow (s_scm_rationalize
);
6109 SCM_WRONG_TYPE_ARG (1, x
);
6113 /* conversion functions */
6116 scm_is_integer (SCM val
)
6118 return scm_is_true (scm_integer_p (val
));
6122 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6124 if (SCM_I_INUMP (val
))
6126 scm_t_signed_bits n
= SCM_I_INUM (val
);
6127 return n
>= min
&& n
<= max
;
6129 else if (SCM_BIGP (val
))
6131 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6133 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6135 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6137 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6138 return n
>= min
&& n
<= max
;
6148 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6149 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6152 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6153 SCM_I_BIG_MPZ (val
));
6155 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6167 return n
>= min
&& n
<= max
;
6175 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6177 if (SCM_I_INUMP (val
))
6179 scm_t_signed_bits n
= SCM_I_INUM (val
);
6180 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6182 else if (SCM_BIGP (val
))
6184 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6186 else if (max
<= ULONG_MAX
)
6188 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6190 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6191 return n
>= min
&& n
<= max
;
6201 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6204 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6205 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6208 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6209 SCM_I_BIG_MPZ (val
));
6211 return n
>= min
&& n
<= max
;
6219 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6221 scm_error (scm_out_of_range_key
,
6223 "Value out of range ~S to ~S: ~S",
6224 scm_list_3 (min
, max
, bad_val
),
6225 scm_list_1 (bad_val
));
6228 #define TYPE scm_t_intmax
6229 #define TYPE_MIN min
6230 #define TYPE_MAX max
6231 #define SIZEOF_TYPE 0
6232 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6233 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6234 #include "libguile/conv-integer.i.c"
6236 #define TYPE scm_t_uintmax
6237 #define TYPE_MIN min
6238 #define TYPE_MAX max
6239 #define SIZEOF_TYPE 0
6240 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6241 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6242 #include "libguile/conv-uinteger.i.c"
6244 #define TYPE scm_t_int8
6245 #define TYPE_MIN SCM_T_INT8_MIN
6246 #define TYPE_MAX SCM_T_INT8_MAX
6247 #define SIZEOF_TYPE 1
6248 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6249 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6250 #include "libguile/conv-integer.i.c"
6252 #define TYPE scm_t_uint8
6254 #define TYPE_MAX SCM_T_UINT8_MAX
6255 #define SIZEOF_TYPE 1
6256 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6257 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6258 #include "libguile/conv-uinteger.i.c"
6260 #define TYPE scm_t_int16
6261 #define TYPE_MIN SCM_T_INT16_MIN
6262 #define TYPE_MAX SCM_T_INT16_MAX
6263 #define SIZEOF_TYPE 2
6264 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6265 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6266 #include "libguile/conv-integer.i.c"
6268 #define TYPE scm_t_uint16
6270 #define TYPE_MAX SCM_T_UINT16_MAX
6271 #define SIZEOF_TYPE 2
6272 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6273 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6274 #include "libguile/conv-uinteger.i.c"
6276 #define TYPE scm_t_int32
6277 #define TYPE_MIN SCM_T_INT32_MIN
6278 #define TYPE_MAX SCM_T_INT32_MAX
6279 #define SIZEOF_TYPE 4
6280 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6281 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6282 #include "libguile/conv-integer.i.c"
6284 #define TYPE scm_t_uint32
6286 #define TYPE_MAX SCM_T_UINT32_MAX
6287 #define SIZEOF_TYPE 4
6288 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6289 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6290 #include "libguile/conv-uinteger.i.c"
6292 #define TYPE scm_t_wchar
6293 #define TYPE_MIN (scm_t_int32)-1
6294 #define TYPE_MAX (scm_t_int32)0x10ffff
6295 #define SIZEOF_TYPE 4
6296 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6297 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6298 #include "libguile/conv-integer.i.c"
6300 #define TYPE scm_t_int64
6301 #define TYPE_MIN SCM_T_INT64_MIN
6302 #define TYPE_MAX SCM_T_INT64_MAX
6303 #define SIZEOF_TYPE 8
6304 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6305 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6306 #include "libguile/conv-integer.i.c"
6308 #define TYPE scm_t_uint64
6310 #define TYPE_MAX SCM_T_UINT64_MAX
6311 #define SIZEOF_TYPE 8
6312 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6313 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6314 #include "libguile/conv-uinteger.i.c"
6317 scm_to_mpz (SCM val
, mpz_t rop
)
6319 if (SCM_I_INUMP (val
))
6320 mpz_set_si (rop
, SCM_I_INUM (val
));
6321 else if (SCM_BIGP (val
))
6322 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6324 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6328 scm_from_mpz (mpz_t val
)
6330 return scm_i_mpz2num (val
);
6334 scm_is_real (SCM val
)
6336 return scm_is_true (scm_real_p (val
));
6340 scm_is_rational (SCM val
)
6342 return scm_is_true (scm_rational_p (val
));
6346 scm_to_double (SCM val
)
6348 if (SCM_I_INUMP (val
))
6349 return SCM_I_INUM (val
);
6350 else if (SCM_BIGP (val
))
6351 return scm_i_big2dbl (val
);
6352 else if (SCM_FRACTIONP (val
))
6353 return scm_i_fraction2double (val
);
6354 else if (SCM_REALP (val
))
6355 return SCM_REAL_VALUE (val
);
6357 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6361 scm_from_double (double val
)
6365 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
6367 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
6368 SCM_REAL_VALUE (z
) = val
;
6373 #if SCM_ENABLE_DEPRECATED == 1
6376 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
6378 scm_c_issue_deprecation_warning
6379 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6383 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6387 scm_out_of_range (NULL
, num
);
6390 return scm_to_double (num
);
6394 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
6396 scm_c_issue_deprecation_warning
6397 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6401 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6405 scm_out_of_range (NULL
, num
);
6408 return scm_to_double (num
);
6414 scm_is_complex (SCM val
)
6416 return scm_is_true (scm_complex_p (val
));
6420 scm_c_real_part (SCM z
)
6422 if (SCM_COMPLEXP (z
))
6423 return SCM_COMPLEX_REAL (z
);
6426 /* Use the scm_real_part to get proper error checking and
6429 return scm_to_double (scm_real_part (z
));
6434 scm_c_imag_part (SCM z
)
6436 if (SCM_COMPLEXP (z
))
6437 return SCM_COMPLEX_IMAG (z
);
6440 /* Use the scm_imag_part to get proper error checking and
6441 dispatching. The result will almost always be 0.0, but not
6444 return scm_to_double (scm_imag_part (z
));
6449 scm_c_magnitude (SCM z
)
6451 return scm_to_double (scm_magnitude (z
));
6457 return scm_to_double (scm_angle (z
));
6461 scm_is_number (SCM z
)
6463 return scm_is_true (scm_number_p (z
));
6467 /* In the following functions we dispatch to the real-arg funcs like log()
6468 when we know the arg is real, instead of just handing everything to
6469 clog() for instance. This is in case clog() doesn't optimize for a
6470 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6471 well use it to go straight to the applicable C func. */
6473 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6475 "Return the natural logarithm of @var{z}.")
6476 #define FUNC_NAME s_scm_log
6478 if (SCM_COMPLEXP (z
))
6480 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6481 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6483 double re
= SCM_COMPLEX_REAL (z
);
6484 double im
= SCM_COMPLEX_IMAG (z
);
6485 return scm_c_make_rectangular (log (hypot (re
, im
)),
6491 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6492 although the value itself overflows. */
6493 double re
= scm_to_double (z
);
6494 double l
= log (fabs (re
));
6496 return scm_from_double (l
);
6498 return scm_c_make_rectangular (l
, M_PI
);
6504 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6506 "Return the base 10 logarithm of @var{z}.")
6507 #define FUNC_NAME s_scm_log10
6509 if (SCM_COMPLEXP (z
))
6511 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6512 clog() and a multiply by M_LOG10E, rather than the fallback
6513 log10+hypot+atan2.) */
6514 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6515 && defined SCM_COMPLEX_VALUE
6516 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6518 double re
= SCM_COMPLEX_REAL (z
);
6519 double im
= SCM_COMPLEX_IMAG (z
);
6520 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6521 M_LOG10E
* atan2 (im
, re
));
6526 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6527 although the value itself overflows. */
6528 double re
= scm_to_double (z
);
6529 double l
= log10 (fabs (re
));
6531 return scm_from_double (l
);
6533 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6539 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6541 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6542 "base of natural logarithms (2.71828@dots{}).")
6543 #define FUNC_NAME s_scm_exp
6545 if (SCM_COMPLEXP (z
))
6547 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6548 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6550 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6551 SCM_COMPLEX_IMAG (z
));
6556 /* When z is a negative bignum the conversion to double overflows,
6557 giving -infinity, but that's ok, the exp is still 0.0. */
6558 return scm_from_double (exp (scm_to_double (z
)));
6564 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6566 "Return the square root of @var{z}. Of the two possible roots\n"
6567 "(positive and negative), the one with the a positive real part\n"
6568 "is returned, or if that's zero then a positive imaginary part.\n"
6572 "(sqrt 9.0) @result{} 3.0\n"
6573 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6574 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6575 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6577 #define FUNC_NAME s_scm_sqrt
6579 if (SCM_COMPLEXP (x
))
6581 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6582 && defined SCM_COMPLEX_VALUE
6583 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6585 double re
= SCM_COMPLEX_REAL (x
);
6586 double im
= SCM_COMPLEX_IMAG (x
);
6587 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6588 0.5 * atan2 (im
, re
));
6593 double xx
= scm_to_double (x
);
6595 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6597 return scm_from_double (sqrt (xx
));
6609 mpz_init_set_si (z_negative_one
, -1);
6611 /* It may be possible to tune the performance of some algorithms by using
6612 * the following constants to avoid the creation of bignums. Please, before
6613 * using these values, remember the two rules of program optimization:
6614 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6615 scm_c_define ("most-positive-fixnum",
6616 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6617 scm_c_define ("most-negative-fixnum",
6618 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6620 scm_add_feature ("complex");
6621 scm_add_feature ("inexact");
6622 flo0
= scm_from_double (0.0);
6624 /* determine floating point precision */
6625 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6627 init_dblprec(&scm_dblprec
[i
-2],i
);
6628 init_fx_radix(fx_per_radix
[i
-2],i
);
6631 /* hard code precision for base 10 if the preprocessor tells us to... */
6632 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6635 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
6636 #include "libguile/numbers.x"