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 ***/
3253 scm_bigequal (SCM x
, SCM y
)
3255 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3256 scm_remember_upto_here_2 (x
, y
);
3257 return scm_from_bool (0 == result
);
3261 scm_real_equalp (SCM x
, SCM y
)
3263 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3267 scm_complex_equalp (SCM x
, SCM y
)
3269 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3270 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3274 scm_i_fraction_equalp (SCM x
, SCM y
)
3276 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3277 SCM_FRACTION_NUMERATOR (y
)))
3278 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3279 SCM_FRACTION_DENOMINATOR (y
))))
3286 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3288 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3290 #define FUNC_NAME s_scm_number_p
3292 return scm_from_bool (SCM_NUMBERP (x
));
3296 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3298 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3299 "otherwise. Note that the sets of real, rational and integer\n"
3300 "values form subsets of the set of complex numbers, i. e. the\n"
3301 "predicate will also be fulfilled if @var{x} is a real,\n"
3302 "rational or integer number.")
3303 #define FUNC_NAME s_scm_complex_p
3305 /* all numbers are complex. */
3306 return scm_number_p (x
);
3310 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3312 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3313 "otherwise. Note that the set of integer values forms a subset of\n"
3314 "the set of real numbers, i. e. the predicate will also be\n"
3315 "fulfilled if @var{x} is an integer number.")
3316 #define FUNC_NAME s_scm_real_p
3318 /* we can't represent irrational numbers. */
3319 return scm_rational_p (x
);
3323 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3325 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3326 "otherwise. Note that the set of integer values forms a subset of\n"
3327 "the set of rational numbers, i. e. the predicate will also be\n"
3328 "fulfilled if @var{x} is an integer number.")
3329 #define FUNC_NAME s_scm_rational_p
3331 if (SCM_I_INUMP (x
))
3333 else if (SCM_IMP (x
))
3335 else if (SCM_BIGP (x
))
3337 else if (SCM_FRACTIONP (x
))
3339 else if (SCM_REALP (x
))
3340 /* due to their limited precision, all floating point numbers are
3341 rational as well. */
3348 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3350 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3352 #define FUNC_NAME s_scm_integer_p
3355 if (SCM_I_INUMP (x
))
3361 if (!SCM_INEXACTP (x
))
3363 if (SCM_COMPLEXP (x
))
3365 r
= SCM_REAL_VALUE (x
);
3375 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3376 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3377 (SCM x
, SCM y
, SCM rest
),
3378 "Return @code{#t} if all parameters are numerically equal.")
3379 #define FUNC_NAME s_scm_i_num_eq_p
3381 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3383 while (!scm_is_null (rest
))
3385 if (scm_is_false (scm_num_eq_p (x
, y
)))
3389 rest
= scm_cdr (rest
);
3391 return scm_num_eq_p (x
, y
);
3395 scm_num_eq_p (SCM x
, SCM y
)
3398 if (SCM_I_INUMP (x
))
3400 scm_t_signed_bits xx
= SCM_I_INUM (x
);
3401 if (SCM_I_INUMP (y
))
3403 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3404 return scm_from_bool (xx
== yy
);
3406 else if (SCM_BIGP (y
))
3408 else if (SCM_REALP (y
))
3410 /* On a 32-bit system an inum fits a double, we can cast the inum
3411 to a double and compare.
3413 But on a 64-bit system an inum is bigger than a double and
3414 casting it to a double (call that dxx) will round. dxx is at
3415 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3416 an integer and fits a long. So we cast yy to a long and
3417 compare with plain xx.
3419 An alternative (for any size system actually) would be to check
3420 yy is an integer (with floor) and is in range of an inum
3421 (compare against appropriate powers of 2) then test
3422 xx==(scm_t_signed_bits)yy. It's just a matter of which
3423 casts/comparisons might be fastest or easiest for the cpu. */
3425 double yy
= SCM_REAL_VALUE (y
);
3426 return scm_from_bool ((double) xx
== yy
3427 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3428 || xx
== (scm_t_signed_bits
) yy
));
3430 else if (SCM_COMPLEXP (y
))
3431 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3432 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3433 else if (SCM_FRACTIONP (y
))
3436 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3438 else if (SCM_BIGP (x
))
3440 if (SCM_I_INUMP (y
))
3442 else if (SCM_BIGP (y
))
3444 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3445 scm_remember_upto_here_2 (x
, y
);
3446 return scm_from_bool (0 == cmp
);
3448 else if (SCM_REALP (y
))
3451 if (isnan (SCM_REAL_VALUE (y
)))
3453 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3454 scm_remember_upto_here_1 (x
);
3455 return scm_from_bool (0 == cmp
);
3457 else if (SCM_COMPLEXP (y
))
3460 if (0.0 != SCM_COMPLEX_IMAG (y
))
3462 if (isnan (SCM_COMPLEX_REAL (y
)))
3464 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3465 scm_remember_upto_here_1 (x
);
3466 return scm_from_bool (0 == cmp
);
3468 else if (SCM_FRACTIONP (y
))
3471 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3473 else if (SCM_REALP (x
))
3475 double xx
= SCM_REAL_VALUE (x
);
3476 if (SCM_I_INUMP (y
))
3478 /* see comments with inum/real above */
3479 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3480 return scm_from_bool (xx
== (double) yy
3481 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3482 || (scm_t_signed_bits
) xx
== yy
));
3484 else if (SCM_BIGP (y
))
3487 if (isnan (SCM_REAL_VALUE (x
)))
3489 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3490 scm_remember_upto_here_1 (y
);
3491 return scm_from_bool (0 == cmp
);
3493 else if (SCM_REALP (y
))
3494 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3495 else if (SCM_COMPLEXP (y
))
3496 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3497 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3498 else if (SCM_FRACTIONP (y
))
3500 double xx
= SCM_REAL_VALUE (x
);
3504 return scm_from_bool (xx
< 0.0);
3505 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3509 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3511 else if (SCM_COMPLEXP (x
))
3513 if (SCM_I_INUMP (y
))
3514 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3515 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3516 else if (SCM_BIGP (y
))
3519 if (0.0 != SCM_COMPLEX_IMAG (x
))
3521 if (isnan (SCM_COMPLEX_REAL (x
)))
3523 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3524 scm_remember_upto_here_1 (y
);
3525 return scm_from_bool (0 == cmp
);
3527 else if (SCM_REALP (y
))
3528 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3529 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3530 else if (SCM_COMPLEXP (y
))
3531 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3532 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3533 else if (SCM_FRACTIONP (y
))
3536 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3538 xx
= SCM_COMPLEX_REAL (x
);
3542 return scm_from_bool (xx
< 0.0);
3543 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3547 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3549 else if (SCM_FRACTIONP (x
))
3551 if (SCM_I_INUMP (y
))
3553 else if (SCM_BIGP (y
))
3555 else if (SCM_REALP (y
))
3557 double yy
= SCM_REAL_VALUE (y
);
3561 return scm_from_bool (0.0 < yy
);
3562 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3565 else if (SCM_COMPLEXP (y
))
3568 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3570 yy
= SCM_COMPLEX_REAL (y
);
3574 return scm_from_bool (0.0 < yy
);
3575 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3578 else if (SCM_FRACTIONP (y
))
3579 return scm_i_fraction_equalp (x
, y
);
3581 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3584 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3588 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3589 done are good for inums, but for bignums an answer can almost always be
3590 had by just examining a few high bits of the operands, as done by GMP in
3591 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3592 of the float exponent to take into account. */
3594 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3595 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3596 (SCM x
, SCM y
, SCM rest
),
3597 "Return @code{#t} if the list of parameters is monotonically\n"
3599 #define FUNC_NAME s_scm_i_num_less_p
3601 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3603 while (!scm_is_null (rest
))
3605 if (scm_is_false (scm_less_p (x
, y
)))
3609 rest
= scm_cdr (rest
);
3611 return scm_less_p (x
, y
);
3615 scm_less_p (SCM x
, SCM y
)
3618 if (SCM_I_INUMP (x
))
3620 scm_t_inum xx
= SCM_I_INUM (x
);
3621 if (SCM_I_INUMP (y
))
3623 scm_t_inum yy
= SCM_I_INUM (y
);
3624 return scm_from_bool (xx
< yy
);
3626 else if (SCM_BIGP (y
))
3628 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3629 scm_remember_upto_here_1 (y
);
3630 return scm_from_bool (sgn
> 0);
3632 else if (SCM_REALP (y
))
3633 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3634 else if (SCM_FRACTIONP (y
))
3636 /* "x < a/b" becomes "x*b < a" */
3638 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3639 y
= SCM_FRACTION_NUMERATOR (y
);
3643 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3645 else if (SCM_BIGP (x
))
3647 if (SCM_I_INUMP (y
))
3649 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3650 scm_remember_upto_here_1 (x
);
3651 return scm_from_bool (sgn
< 0);
3653 else if (SCM_BIGP (y
))
3655 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3656 scm_remember_upto_here_2 (x
, y
);
3657 return scm_from_bool (cmp
< 0);
3659 else if (SCM_REALP (y
))
3662 if (isnan (SCM_REAL_VALUE (y
)))
3664 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3665 scm_remember_upto_here_1 (x
);
3666 return scm_from_bool (cmp
< 0);
3668 else if (SCM_FRACTIONP (y
))
3671 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3673 else if (SCM_REALP (x
))
3675 if (SCM_I_INUMP (y
))
3676 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3677 else if (SCM_BIGP (y
))
3680 if (isnan (SCM_REAL_VALUE (x
)))
3682 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3683 scm_remember_upto_here_1 (y
);
3684 return scm_from_bool (cmp
> 0);
3686 else if (SCM_REALP (y
))
3687 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3688 else if (SCM_FRACTIONP (y
))
3690 double xx
= SCM_REAL_VALUE (x
);
3694 return scm_from_bool (xx
< 0.0);
3695 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3699 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3701 else if (SCM_FRACTIONP (x
))
3703 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3705 /* "a/b < y" becomes "a < y*b" */
3706 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3707 x
= SCM_FRACTION_NUMERATOR (x
);
3710 else if (SCM_REALP (y
))
3712 double yy
= SCM_REAL_VALUE (y
);
3716 return scm_from_bool (0.0 < yy
);
3717 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3720 else if (SCM_FRACTIONP (y
))
3722 /* "a/b < c/d" becomes "a*d < c*b" */
3723 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3724 SCM_FRACTION_DENOMINATOR (y
));
3725 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3726 SCM_FRACTION_DENOMINATOR (x
));
3732 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3735 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3739 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3740 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3741 (SCM x
, SCM y
, SCM rest
),
3742 "Return @code{#t} if the list of parameters is monotonically\n"
3744 #define FUNC_NAME s_scm_i_num_gr_p
3746 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3748 while (!scm_is_null (rest
))
3750 if (scm_is_false (scm_gr_p (x
, y
)))
3754 rest
= scm_cdr (rest
);
3756 return scm_gr_p (x
, y
);
3759 #define FUNC_NAME s_scm_i_num_gr_p
3761 scm_gr_p (SCM x
, SCM y
)
3763 if (!SCM_NUMBERP (x
))
3764 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3765 else if (!SCM_NUMBERP (y
))
3766 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3768 return scm_less_p (y
, x
);
3773 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3774 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3775 (SCM x
, SCM y
, SCM rest
),
3776 "Return @code{#t} if the list of parameters is monotonically\n"
3778 #define FUNC_NAME s_scm_i_num_leq_p
3780 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3782 while (!scm_is_null (rest
))
3784 if (scm_is_false (scm_leq_p (x
, y
)))
3788 rest
= scm_cdr (rest
);
3790 return scm_leq_p (x
, y
);
3793 #define FUNC_NAME s_scm_i_num_leq_p
3795 scm_leq_p (SCM x
, SCM y
)
3797 if (!SCM_NUMBERP (x
))
3798 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3799 else if (!SCM_NUMBERP (y
))
3800 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3801 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3804 return scm_not (scm_less_p (y
, x
));
3809 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3810 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3811 (SCM x
, SCM y
, SCM rest
),
3812 "Return @code{#t} if the list of parameters is monotonically\n"
3814 #define FUNC_NAME s_scm_i_num_geq_p
3816 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3818 while (!scm_is_null (rest
))
3820 if (scm_is_false (scm_geq_p (x
, y
)))
3824 rest
= scm_cdr (rest
);
3826 return scm_geq_p (x
, y
);
3829 #define FUNC_NAME s_scm_i_num_geq_p
3831 scm_geq_p (SCM x
, SCM y
)
3833 if (!SCM_NUMBERP (x
))
3834 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3835 else if (!SCM_NUMBERP (y
))
3836 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3837 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3840 return scm_not (scm_less_p (x
, y
));
3845 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3846 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3852 if (SCM_I_INUMP (z
))
3853 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3854 else if (SCM_BIGP (z
))
3856 else if (SCM_REALP (z
))
3857 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3858 else if (SCM_COMPLEXP (z
))
3859 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3860 && SCM_COMPLEX_IMAG (z
) == 0.0);
3861 else if (SCM_FRACTIONP (z
))
3864 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3868 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3869 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3873 scm_positive_p (SCM x
)
3875 if (SCM_I_INUMP (x
))
3876 return scm_from_bool (SCM_I_INUM (x
) > 0);
3877 else if (SCM_BIGP (x
))
3879 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3880 scm_remember_upto_here_1 (x
);
3881 return scm_from_bool (sgn
> 0);
3883 else if (SCM_REALP (x
))
3884 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3885 else if (SCM_FRACTIONP (x
))
3886 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3888 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3892 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3893 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3897 scm_negative_p (SCM x
)
3899 if (SCM_I_INUMP (x
))
3900 return scm_from_bool (SCM_I_INUM (x
) < 0);
3901 else if (SCM_BIGP (x
))
3903 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3904 scm_remember_upto_here_1 (x
);
3905 return scm_from_bool (sgn
< 0);
3907 else if (SCM_REALP (x
))
3908 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3909 else if (SCM_FRACTIONP (x
))
3910 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3912 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3916 /* scm_min and scm_max return an inexact when either argument is inexact, as
3917 required by r5rs. On that basis, for exact/inexact combinations the
3918 exact is converted to inexact to compare and possibly return. This is
3919 unlike scm_less_p above which takes some trouble to preserve all bits in
3920 its test, such trouble is not required for min and max. */
3922 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3923 (SCM x
, SCM y
, SCM rest
),
3924 "Return the maximum of all parameter values.")
3925 #define FUNC_NAME s_scm_i_max
3927 while (!scm_is_null (rest
))
3928 { x
= scm_max (x
, y
);
3930 rest
= scm_cdr (rest
);
3932 return scm_max (x
, y
);
3936 #define s_max s_scm_i_max
3937 #define g_max g_scm_i_max
3940 scm_max (SCM x
, SCM y
)
3945 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3946 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3949 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3952 if (SCM_I_INUMP (x
))
3954 scm_t_inum xx
= SCM_I_INUM (x
);
3955 if (SCM_I_INUMP (y
))
3957 scm_t_inum yy
= SCM_I_INUM (y
);
3958 return (xx
< yy
) ? y
: x
;
3960 else if (SCM_BIGP (y
))
3962 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3963 scm_remember_upto_here_1 (y
);
3964 return (sgn
< 0) ? x
: y
;
3966 else if (SCM_REALP (y
))
3969 /* if y==NaN then ">" is false and we return NaN */
3970 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3972 else if (SCM_FRACTIONP (y
))
3975 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3978 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3980 else if (SCM_BIGP (x
))
3982 if (SCM_I_INUMP (y
))
3984 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3985 scm_remember_upto_here_1 (x
);
3986 return (sgn
< 0) ? y
: x
;
3988 else if (SCM_BIGP (y
))
3990 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3991 scm_remember_upto_here_2 (x
, y
);
3992 return (cmp
> 0) ? x
: y
;
3994 else if (SCM_REALP (y
))
3996 /* if y==NaN then xx>yy is false, so we return the NaN y */
3999 xx
= scm_i_big2dbl (x
);
4000 yy
= SCM_REAL_VALUE (y
);
4001 return (xx
> yy
? scm_from_double (xx
) : y
);
4003 else if (SCM_FRACTIONP (y
))
4008 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4010 else if (SCM_REALP (x
))
4012 if (SCM_I_INUMP (y
))
4014 double z
= SCM_I_INUM (y
);
4015 /* if x==NaN then "<" is false and we return NaN */
4016 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
4018 else if (SCM_BIGP (y
))
4023 else if (SCM_REALP (y
))
4025 /* if x==NaN then our explicit check means we return NaN
4026 if y==NaN then ">" is false and we return NaN
4027 calling isnan is unavoidable, since it's the only way to know
4028 which of x or y causes any compares to be false */
4029 double xx
= SCM_REAL_VALUE (x
);
4030 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
4032 else if (SCM_FRACTIONP (y
))
4034 double yy
= scm_i_fraction2double (y
);
4035 double xx
= SCM_REAL_VALUE (x
);
4036 return (xx
< yy
) ? scm_from_double (yy
) : x
;
4039 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4041 else if (SCM_FRACTIONP (x
))
4043 if (SCM_I_INUMP (y
))
4047 else if (SCM_BIGP (y
))
4051 else if (SCM_REALP (y
))
4053 double xx
= scm_i_fraction2double (x
);
4054 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4056 else if (SCM_FRACTIONP (y
))
4061 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4064 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4068 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4069 (SCM x
, SCM y
, SCM rest
),
4070 "Return the minimum of all parameter values.")
4071 #define FUNC_NAME s_scm_i_min
4073 while (!scm_is_null (rest
))
4074 { x
= scm_min (x
, y
);
4076 rest
= scm_cdr (rest
);
4078 return scm_min (x
, y
);
4082 #define s_min s_scm_i_min
4083 #define g_min g_scm_i_min
4086 scm_min (SCM x
, SCM y
)
4091 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4092 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4095 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4098 if (SCM_I_INUMP (x
))
4100 scm_t_inum xx
= SCM_I_INUM (x
);
4101 if (SCM_I_INUMP (y
))
4103 scm_t_inum yy
= SCM_I_INUM (y
);
4104 return (xx
< yy
) ? x
: y
;
4106 else if (SCM_BIGP (y
))
4108 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4109 scm_remember_upto_here_1 (y
);
4110 return (sgn
< 0) ? y
: x
;
4112 else if (SCM_REALP (y
))
4115 /* if y==NaN then "<" is false and we return NaN */
4116 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4118 else if (SCM_FRACTIONP (y
))
4121 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4124 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4126 else if (SCM_BIGP (x
))
4128 if (SCM_I_INUMP (y
))
4130 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4131 scm_remember_upto_here_1 (x
);
4132 return (sgn
< 0) ? x
: y
;
4134 else if (SCM_BIGP (y
))
4136 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4137 scm_remember_upto_here_2 (x
, y
);
4138 return (cmp
> 0) ? y
: x
;
4140 else if (SCM_REALP (y
))
4142 /* if y==NaN then xx<yy is false, so we return the NaN y */
4145 xx
= scm_i_big2dbl (x
);
4146 yy
= SCM_REAL_VALUE (y
);
4147 return (xx
< yy
? scm_from_double (xx
) : y
);
4149 else if (SCM_FRACTIONP (y
))
4154 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4156 else if (SCM_REALP (x
))
4158 if (SCM_I_INUMP (y
))
4160 double z
= SCM_I_INUM (y
);
4161 /* if x==NaN then "<" is false and we return NaN */
4162 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4164 else if (SCM_BIGP (y
))
4169 else if (SCM_REALP (y
))
4171 /* if x==NaN then our explicit check means we return NaN
4172 if y==NaN then "<" is false and we return NaN
4173 calling isnan is unavoidable, since it's the only way to know
4174 which of x or y causes any compares to be false */
4175 double xx
= SCM_REAL_VALUE (x
);
4176 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4178 else if (SCM_FRACTIONP (y
))
4180 double yy
= scm_i_fraction2double (y
);
4181 double xx
= SCM_REAL_VALUE (x
);
4182 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4185 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4187 else if (SCM_FRACTIONP (x
))
4189 if (SCM_I_INUMP (y
))
4193 else if (SCM_BIGP (y
))
4197 else if (SCM_REALP (y
))
4199 double xx
= scm_i_fraction2double (x
);
4200 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4202 else if (SCM_FRACTIONP (y
))
4207 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4210 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4214 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4215 (SCM x
, SCM y
, SCM rest
),
4216 "Return the sum of all parameter values. Return 0 if called without\n"
4218 #define FUNC_NAME s_scm_i_sum
4220 while (!scm_is_null (rest
))
4221 { x
= scm_sum (x
, y
);
4223 rest
= scm_cdr (rest
);
4225 return scm_sum (x
, y
);
4229 #define s_sum s_scm_i_sum
4230 #define g_sum g_scm_i_sum
4233 scm_sum (SCM x
, SCM y
)
4235 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4237 if (SCM_NUMBERP (x
)) return x
;
4238 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4239 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4242 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4244 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4246 scm_t_inum xx
= SCM_I_INUM (x
);
4247 scm_t_inum yy
= SCM_I_INUM (y
);
4248 scm_t_inum z
= xx
+ yy
;
4249 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
4251 else if (SCM_BIGP (y
))
4256 else if (SCM_REALP (y
))
4258 scm_t_inum xx
= SCM_I_INUM (x
);
4259 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4261 else if (SCM_COMPLEXP (y
))
4263 scm_t_inum xx
= SCM_I_INUM (x
);
4264 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4265 SCM_COMPLEX_IMAG (y
));
4267 else if (SCM_FRACTIONP (y
))
4268 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4269 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4270 SCM_FRACTION_DENOMINATOR (y
));
4272 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4273 } else if (SCM_BIGP (x
))
4275 if (SCM_I_INUMP (y
))
4280 inum
= SCM_I_INUM (y
);
4283 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4286 SCM result
= scm_i_mkbig ();
4287 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4288 scm_remember_upto_here_1 (x
);
4289 /* we know the result will have to be a bignum */
4292 return scm_i_normbig (result
);
4296 SCM result
= scm_i_mkbig ();
4297 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4298 scm_remember_upto_here_1 (x
);
4299 /* we know the result will have to be a bignum */
4302 return scm_i_normbig (result
);
4305 else if (SCM_BIGP (y
))
4307 SCM result
= scm_i_mkbig ();
4308 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4309 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4310 mpz_add (SCM_I_BIG_MPZ (result
),
4313 scm_remember_upto_here_2 (x
, y
);
4314 /* we know the result will have to be a bignum */
4317 return scm_i_normbig (result
);
4319 else if (SCM_REALP (y
))
4321 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4322 scm_remember_upto_here_1 (x
);
4323 return scm_from_double (result
);
4325 else if (SCM_COMPLEXP (y
))
4327 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4328 + SCM_COMPLEX_REAL (y
));
4329 scm_remember_upto_here_1 (x
);
4330 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4332 else if (SCM_FRACTIONP (y
))
4333 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4334 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4335 SCM_FRACTION_DENOMINATOR (y
));
4337 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4339 else if (SCM_REALP (x
))
4341 if (SCM_I_INUMP (y
))
4342 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4343 else if (SCM_BIGP (y
))
4345 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4346 scm_remember_upto_here_1 (y
);
4347 return scm_from_double (result
);
4349 else if (SCM_REALP (y
))
4350 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4351 else if (SCM_COMPLEXP (y
))
4352 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4353 SCM_COMPLEX_IMAG (y
));
4354 else if (SCM_FRACTIONP (y
))
4355 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4357 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4359 else if (SCM_COMPLEXP (x
))
4361 if (SCM_I_INUMP (y
))
4362 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4363 SCM_COMPLEX_IMAG (x
));
4364 else if (SCM_BIGP (y
))
4366 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4367 + SCM_COMPLEX_REAL (x
));
4368 scm_remember_upto_here_1 (y
);
4369 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4371 else if (SCM_REALP (y
))
4372 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4373 SCM_COMPLEX_IMAG (x
));
4374 else if (SCM_COMPLEXP (y
))
4375 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4376 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4377 else if (SCM_FRACTIONP (y
))
4378 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4379 SCM_COMPLEX_IMAG (x
));
4381 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4383 else if (SCM_FRACTIONP (x
))
4385 if (SCM_I_INUMP (y
))
4386 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4387 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4388 SCM_FRACTION_DENOMINATOR (x
));
4389 else if (SCM_BIGP (y
))
4390 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4391 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4392 SCM_FRACTION_DENOMINATOR (x
));
4393 else if (SCM_REALP (y
))
4394 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4395 else if (SCM_COMPLEXP (y
))
4396 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4397 SCM_COMPLEX_IMAG (y
));
4398 else if (SCM_FRACTIONP (y
))
4399 /* a/b + c/d = (ad + bc) / bd */
4400 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4401 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4402 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4404 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4407 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4411 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4413 "Return @math{@var{x}+1}.")
4414 #define FUNC_NAME s_scm_oneplus
4416 return scm_sum (x
, SCM_INUM1
);
4421 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4422 (SCM x
, SCM y
, SCM rest
),
4423 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4424 "the sum of all but the first argument are subtracted from the first\n"
4426 #define FUNC_NAME s_scm_i_difference
4428 while (!scm_is_null (rest
))
4429 { x
= scm_difference (x
, y
);
4431 rest
= scm_cdr (rest
);
4433 return scm_difference (x
, y
);
4437 #define s_difference s_scm_i_difference
4438 #define g_difference g_scm_i_difference
4441 scm_difference (SCM x
, SCM y
)
4442 #define FUNC_NAME s_difference
4444 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4447 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4449 if (SCM_I_INUMP (x
))
4451 scm_t_inum xx
= -SCM_I_INUM (x
);
4452 if (SCM_FIXABLE (xx
))
4453 return SCM_I_MAKINUM (xx
);
4455 return scm_i_inum2big (xx
);
4457 else if (SCM_BIGP (x
))
4458 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4459 bignum, but negating that gives a fixnum. */
4460 return scm_i_normbig (scm_i_clonebig (x
, 0));
4461 else if (SCM_REALP (x
))
4462 return scm_from_double (-SCM_REAL_VALUE (x
));
4463 else if (SCM_COMPLEXP (x
))
4464 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4465 -SCM_COMPLEX_IMAG (x
));
4466 else if (SCM_FRACTIONP (x
))
4467 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4468 SCM_FRACTION_DENOMINATOR (x
));
4470 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4473 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4475 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4477 scm_t_inum xx
= SCM_I_INUM (x
);
4478 scm_t_inum yy
= SCM_I_INUM (y
);
4479 scm_t_inum z
= xx
- yy
;
4480 if (SCM_FIXABLE (z
))
4481 return SCM_I_MAKINUM (z
);
4483 return scm_i_inum2big (z
);
4485 else if (SCM_BIGP (y
))
4487 /* inum-x - big-y */
4488 scm_t_inum xx
= SCM_I_INUM (x
);
4491 return scm_i_clonebig (y
, 0);
4494 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4495 SCM result
= scm_i_mkbig ();
4498 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4501 /* x - y == -(y + -x) */
4502 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4503 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4505 scm_remember_upto_here_1 (y
);
4507 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4508 /* we know the result will have to be a bignum */
4511 return scm_i_normbig (result
);
4514 else if (SCM_REALP (y
))
4516 scm_t_inum xx
= SCM_I_INUM (x
);
4517 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4519 else if (SCM_COMPLEXP (y
))
4521 scm_t_inum xx
= SCM_I_INUM (x
);
4522 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4523 - SCM_COMPLEX_IMAG (y
));
4525 else if (SCM_FRACTIONP (y
))
4526 /* a - b/c = (ac - b) / c */
4527 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4528 SCM_FRACTION_NUMERATOR (y
)),
4529 SCM_FRACTION_DENOMINATOR (y
));
4531 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4533 else if (SCM_BIGP (x
))
4535 if (SCM_I_INUMP (y
))
4537 /* big-x - inum-y */
4538 scm_t_inum yy
= SCM_I_INUM (y
);
4539 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4541 scm_remember_upto_here_1 (x
);
4543 return (SCM_FIXABLE (-yy
) ?
4544 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
4547 SCM result
= scm_i_mkbig ();
4550 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4552 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4553 scm_remember_upto_here_1 (x
);
4555 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4556 /* we know the result will have to be a bignum */
4559 return scm_i_normbig (result
);
4562 else if (SCM_BIGP (y
))
4564 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4565 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4566 SCM result
= scm_i_mkbig ();
4567 mpz_sub (SCM_I_BIG_MPZ (result
),
4570 scm_remember_upto_here_2 (x
, y
);
4571 /* we know the result will have to be a bignum */
4572 if ((sgn_x
== 1) && (sgn_y
== -1))
4574 if ((sgn_x
== -1) && (sgn_y
== 1))
4576 return scm_i_normbig (result
);
4578 else if (SCM_REALP (y
))
4580 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4581 scm_remember_upto_here_1 (x
);
4582 return scm_from_double (result
);
4584 else if (SCM_COMPLEXP (y
))
4586 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4587 - SCM_COMPLEX_REAL (y
));
4588 scm_remember_upto_here_1 (x
);
4589 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4591 else if (SCM_FRACTIONP (y
))
4592 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4593 SCM_FRACTION_NUMERATOR (y
)),
4594 SCM_FRACTION_DENOMINATOR (y
));
4595 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4597 else if (SCM_REALP (x
))
4599 if (SCM_I_INUMP (y
))
4600 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4601 else if (SCM_BIGP (y
))
4603 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4604 scm_remember_upto_here_1 (x
);
4605 return scm_from_double (result
);
4607 else if (SCM_REALP (y
))
4608 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4609 else if (SCM_COMPLEXP (y
))
4610 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4611 -SCM_COMPLEX_IMAG (y
));
4612 else if (SCM_FRACTIONP (y
))
4613 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4615 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4617 else if (SCM_COMPLEXP (x
))
4619 if (SCM_I_INUMP (y
))
4620 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4621 SCM_COMPLEX_IMAG (x
));
4622 else if (SCM_BIGP (y
))
4624 double real_part
= (SCM_COMPLEX_REAL (x
)
4625 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4626 scm_remember_upto_here_1 (x
);
4627 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4629 else if (SCM_REALP (y
))
4630 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4631 SCM_COMPLEX_IMAG (x
));
4632 else if (SCM_COMPLEXP (y
))
4633 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4634 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4635 else if (SCM_FRACTIONP (y
))
4636 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4637 SCM_COMPLEX_IMAG (x
));
4639 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4641 else if (SCM_FRACTIONP (x
))
4643 if (SCM_I_INUMP (y
))
4644 /* a/b - c = (a - cb) / b */
4645 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4646 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4647 SCM_FRACTION_DENOMINATOR (x
));
4648 else if (SCM_BIGP (y
))
4649 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4650 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4651 SCM_FRACTION_DENOMINATOR (x
));
4652 else if (SCM_REALP (y
))
4653 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4654 else if (SCM_COMPLEXP (y
))
4655 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4656 -SCM_COMPLEX_IMAG (y
));
4657 else if (SCM_FRACTIONP (y
))
4658 /* a/b - c/d = (ad - bc) / bd */
4659 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4660 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4661 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4663 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4666 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4671 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4673 "Return @math{@var{x}-1}.")
4674 #define FUNC_NAME s_scm_oneminus
4676 return scm_difference (x
, SCM_INUM1
);
4681 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4682 (SCM x
, SCM y
, SCM rest
),
4683 "Return the product of all arguments. If called without arguments,\n"
4685 #define FUNC_NAME s_scm_i_product
4687 while (!scm_is_null (rest
))
4688 { x
= scm_product (x
, y
);
4690 rest
= scm_cdr (rest
);
4692 return scm_product (x
, y
);
4696 #define s_product s_scm_i_product
4697 #define g_product g_scm_i_product
4700 scm_product (SCM x
, SCM y
)
4702 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4705 return SCM_I_MAKINUM (1L);
4706 else if (SCM_NUMBERP (x
))
4709 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4712 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4717 xx
= SCM_I_INUM (x
);
4721 case 0: return x
; break;
4722 case 1: return y
; break;
4725 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4727 scm_t_inum yy
= SCM_I_INUM (y
);
4728 scm_t_inum kk
= xx
* yy
;
4729 SCM k
= SCM_I_MAKINUM (kk
);
4730 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4734 SCM result
= scm_i_inum2big (xx
);
4735 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4736 return scm_i_normbig (result
);
4739 else if (SCM_BIGP (y
))
4741 SCM result
= scm_i_mkbig ();
4742 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4743 scm_remember_upto_here_1 (y
);
4746 else if (SCM_REALP (y
))
4747 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4748 else if (SCM_COMPLEXP (y
))
4749 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4750 xx
* SCM_COMPLEX_IMAG (y
));
4751 else if (SCM_FRACTIONP (y
))
4752 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4753 SCM_FRACTION_DENOMINATOR (y
));
4755 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4757 else if (SCM_BIGP (x
))
4759 if (SCM_I_INUMP (y
))
4764 else if (SCM_BIGP (y
))
4766 SCM result
= scm_i_mkbig ();
4767 mpz_mul (SCM_I_BIG_MPZ (result
),
4770 scm_remember_upto_here_2 (x
, y
);
4773 else if (SCM_REALP (y
))
4775 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4776 scm_remember_upto_here_1 (x
);
4777 return scm_from_double (result
);
4779 else if (SCM_COMPLEXP (y
))
4781 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4782 scm_remember_upto_here_1 (x
);
4783 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4784 z
* SCM_COMPLEX_IMAG (y
));
4786 else if (SCM_FRACTIONP (y
))
4787 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4788 SCM_FRACTION_DENOMINATOR (y
));
4790 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4792 else if (SCM_REALP (x
))
4794 if (SCM_I_INUMP (y
))
4796 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4797 if (scm_is_eq (y
, SCM_INUM0
))
4799 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4801 else if (SCM_BIGP (y
))
4803 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4804 scm_remember_upto_here_1 (y
);
4805 return scm_from_double (result
);
4807 else if (SCM_REALP (y
))
4808 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4809 else if (SCM_COMPLEXP (y
))
4810 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4811 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4812 else if (SCM_FRACTIONP (y
))
4813 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4815 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4817 else if (SCM_COMPLEXP (x
))
4819 if (SCM_I_INUMP (y
))
4821 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4822 if (scm_is_eq (y
, SCM_INUM0
))
4824 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4825 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4827 else if (SCM_BIGP (y
))
4829 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4830 scm_remember_upto_here_1 (y
);
4831 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4832 z
* SCM_COMPLEX_IMAG (x
));
4834 else if (SCM_REALP (y
))
4835 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4836 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4837 else if (SCM_COMPLEXP (y
))
4839 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4840 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4841 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4842 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4844 else if (SCM_FRACTIONP (y
))
4846 double yy
= scm_i_fraction2double (y
);
4847 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4848 yy
* SCM_COMPLEX_IMAG (x
));
4851 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4853 else if (SCM_FRACTIONP (x
))
4855 if (SCM_I_INUMP (y
))
4856 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4857 SCM_FRACTION_DENOMINATOR (x
));
4858 else if (SCM_BIGP (y
))
4859 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4860 SCM_FRACTION_DENOMINATOR (x
));
4861 else if (SCM_REALP (y
))
4862 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4863 else if (SCM_COMPLEXP (y
))
4865 double xx
= scm_i_fraction2double (x
);
4866 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4867 xx
* SCM_COMPLEX_IMAG (y
));
4869 else if (SCM_FRACTIONP (y
))
4870 /* a/b * c/d = ac / bd */
4871 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4872 SCM_FRACTION_NUMERATOR (y
)),
4873 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4874 SCM_FRACTION_DENOMINATOR (y
)));
4876 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4879 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4882 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4883 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4884 #define ALLOW_DIVIDE_BY_ZERO
4885 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4888 /* The code below for complex division is adapted from the GNU
4889 libstdc++, which adapted it from f2c's libF77, and is subject to
4892 /****************************************************************
4893 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4895 Permission to use, copy, modify, and distribute this software
4896 and its documentation for any purpose and without fee is hereby
4897 granted, provided that the above copyright notice appear in all
4898 copies and that both that the copyright notice and this
4899 permission notice and warranty disclaimer appear in supporting
4900 documentation, and that the names of AT&T Bell Laboratories or
4901 Bellcore or any of their entities not be used in advertising or
4902 publicity pertaining to distribution of the software without
4903 specific, written prior permission.
4905 AT&T and Bellcore disclaim all warranties with regard to this
4906 software, including all implied warranties of merchantability
4907 and fitness. In no event shall AT&T or Bellcore be liable for
4908 any special, indirect or consequential damages or any damages
4909 whatsoever resulting from loss of use, data or profits, whether
4910 in an action of contract, negligence or other tortious action,
4911 arising out of or in connection with the use or performance of
4913 ****************************************************************/
4915 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4916 (SCM x
, SCM y
, SCM rest
),
4917 "Divide the first argument by the product of the remaining\n"
4918 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4920 #define FUNC_NAME s_scm_i_divide
4922 while (!scm_is_null (rest
))
4923 { x
= scm_divide (x
, y
);
4925 rest
= scm_cdr (rest
);
4927 return scm_divide (x
, y
);
4931 #define s_divide s_scm_i_divide
4932 #define g_divide g_scm_i_divide
4935 do_divide (SCM x
, SCM y
, int inexact
)
4936 #define FUNC_NAME s_divide
4940 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4943 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4944 else if (SCM_I_INUMP (x
))
4946 scm_t_inum xx
= SCM_I_INUM (x
);
4947 if (xx
== 1 || xx
== -1)
4949 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4951 scm_num_overflow (s_divide
);
4956 return scm_from_double (1.0 / (double) xx
);
4957 else return scm_i_make_ratio (SCM_INUM1
, x
);
4960 else if (SCM_BIGP (x
))
4963 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4964 else return scm_i_make_ratio (SCM_INUM1
, x
);
4966 else if (SCM_REALP (x
))
4968 double xx
= SCM_REAL_VALUE (x
);
4969 #ifndef ALLOW_DIVIDE_BY_ZERO
4971 scm_num_overflow (s_divide
);
4974 return scm_from_double (1.0 / xx
);
4976 else if (SCM_COMPLEXP (x
))
4978 double r
= SCM_COMPLEX_REAL (x
);
4979 double i
= SCM_COMPLEX_IMAG (x
);
4980 if (fabs(r
) <= fabs(i
))
4983 double d
= i
* (1.0 + t
* t
);
4984 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4989 double d
= r
* (1.0 + t
* t
);
4990 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4993 else if (SCM_FRACTIONP (x
))
4994 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4995 SCM_FRACTION_NUMERATOR (x
));
4997 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
5000 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5002 scm_t_inum xx
= SCM_I_INUM (x
);
5003 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5005 scm_t_inum yy
= SCM_I_INUM (y
);
5008 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5009 scm_num_overflow (s_divide
);
5011 return scm_from_double ((double) xx
/ (double) yy
);
5014 else if (xx
% yy
!= 0)
5017 return scm_from_double ((double) xx
/ (double) yy
);
5018 else return scm_i_make_ratio (x
, y
);
5022 scm_t_inum z
= xx
/ yy
;
5023 if (SCM_FIXABLE (z
))
5024 return SCM_I_MAKINUM (z
);
5026 return scm_i_inum2big (z
);
5029 else if (SCM_BIGP (y
))
5032 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
5033 else return scm_i_make_ratio (x
, y
);
5035 else if (SCM_REALP (y
))
5037 double yy
= SCM_REAL_VALUE (y
);
5038 #ifndef ALLOW_DIVIDE_BY_ZERO
5040 scm_num_overflow (s_divide
);
5043 return scm_from_double ((double) xx
/ yy
);
5045 else if (SCM_COMPLEXP (y
))
5048 complex_div
: /* y _must_ be a complex number */
5050 double r
= SCM_COMPLEX_REAL (y
);
5051 double i
= SCM_COMPLEX_IMAG (y
);
5052 if (fabs(r
) <= fabs(i
))
5055 double d
= i
* (1.0 + t
* t
);
5056 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5061 double d
= r
* (1.0 + t
* t
);
5062 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5066 else if (SCM_FRACTIONP (y
))
5067 /* a / b/c = ac / b */
5068 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5069 SCM_FRACTION_NUMERATOR (y
));
5071 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5073 else if (SCM_BIGP (x
))
5075 if (SCM_I_INUMP (y
))
5077 scm_t_inum yy
= SCM_I_INUM (y
);
5080 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5081 scm_num_overflow (s_divide
);
5083 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5084 scm_remember_upto_here_1 (x
);
5085 return (sgn
== 0) ? scm_nan () : scm_inf ();
5092 /* FIXME: HMM, what are the relative performance issues here?
5093 We need to test. Is it faster on average to test
5094 divisible_p, then perform whichever operation, or is it
5095 faster to perform the integer div opportunistically and
5096 switch to real if there's a remainder? For now we take the
5097 middle ground: test, then if divisible, use the faster div
5100 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
5101 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5105 SCM result
= scm_i_mkbig ();
5106 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5107 scm_remember_upto_here_1 (x
);
5109 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5110 return scm_i_normbig (result
);
5115 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5116 else return scm_i_make_ratio (x
, y
);
5120 else if (SCM_BIGP (y
))
5125 /* It's easily possible for the ratio x/y to fit a double
5126 but one or both x and y be too big to fit a double,
5127 hence the use of mpq_get_d rather than converting and
5130 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5131 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5132 return scm_from_double (mpq_get_d (q
));
5136 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5140 SCM result
= scm_i_mkbig ();
5141 mpz_divexact (SCM_I_BIG_MPZ (result
),
5144 scm_remember_upto_here_2 (x
, y
);
5145 return scm_i_normbig (result
);
5148 return scm_i_make_ratio (x
, y
);
5151 else if (SCM_REALP (y
))
5153 double yy
= SCM_REAL_VALUE (y
);
5154 #ifndef ALLOW_DIVIDE_BY_ZERO
5156 scm_num_overflow (s_divide
);
5159 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5161 else if (SCM_COMPLEXP (y
))
5163 a
= scm_i_big2dbl (x
);
5166 else if (SCM_FRACTIONP (y
))
5167 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5168 SCM_FRACTION_NUMERATOR (y
));
5170 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5172 else if (SCM_REALP (x
))
5174 double rx
= SCM_REAL_VALUE (x
);
5175 if (SCM_I_INUMP (y
))
5177 scm_t_inum yy
= SCM_I_INUM (y
);
5178 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5180 scm_num_overflow (s_divide
);
5183 return scm_from_double (rx
/ (double) yy
);
5185 else if (SCM_BIGP (y
))
5187 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5188 scm_remember_upto_here_1 (y
);
5189 return scm_from_double (rx
/ dby
);
5191 else if (SCM_REALP (y
))
5193 double yy
= SCM_REAL_VALUE (y
);
5194 #ifndef ALLOW_DIVIDE_BY_ZERO
5196 scm_num_overflow (s_divide
);
5199 return scm_from_double (rx
/ yy
);
5201 else if (SCM_COMPLEXP (y
))
5206 else if (SCM_FRACTIONP (y
))
5207 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5209 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5211 else if (SCM_COMPLEXP (x
))
5213 double rx
= SCM_COMPLEX_REAL (x
);
5214 double ix
= SCM_COMPLEX_IMAG (x
);
5215 if (SCM_I_INUMP (y
))
5217 scm_t_inum yy
= SCM_I_INUM (y
);
5218 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5220 scm_num_overflow (s_divide
);
5225 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5228 else if (SCM_BIGP (y
))
5230 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5231 scm_remember_upto_here_1 (y
);
5232 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5234 else if (SCM_REALP (y
))
5236 double yy
= SCM_REAL_VALUE (y
);
5237 #ifndef ALLOW_DIVIDE_BY_ZERO
5239 scm_num_overflow (s_divide
);
5242 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5244 else if (SCM_COMPLEXP (y
))
5246 double ry
= SCM_COMPLEX_REAL (y
);
5247 double iy
= SCM_COMPLEX_IMAG (y
);
5248 if (fabs(ry
) <= fabs(iy
))
5251 double d
= iy
* (1.0 + t
* t
);
5252 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5257 double d
= ry
* (1.0 + t
* t
);
5258 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5261 else if (SCM_FRACTIONP (y
))
5263 double yy
= scm_i_fraction2double (y
);
5264 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5267 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5269 else if (SCM_FRACTIONP (x
))
5271 if (SCM_I_INUMP (y
))
5273 scm_t_inum yy
= SCM_I_INUM (y
);
5274 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5276 scm_num_overflow (s_divide
);
5279 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5280 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5282 else if (SCM_BIGP (y
))
5284 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5285 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5287 else if (SCM_REALP (y
))
5289 double yy
= SCM_REAL_VALUE (y
);
5290 #ifndef ALLOW_DIVIDE_BY_ZERO
5292 scm_num_overflow (s_divide
);
5295 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5297 else if (SCM_COMPLEXP (y
))
5299 a
= scm_i_fraction2double (x
);
5302 else if (SCM_FRACTIONP (y
))
5303 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5304 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5306 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5309 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5313 scm_divide (SCM x
, SCM y
)
5315 return do_divide (x
, y
, 0);
5318 static SCM
scm_divide2real (SCM x
, SCM y
)
5320 return do_divide (x
, y
, 1);
5326 scm_c_truncate (double x
)
5337 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5338 half-way case (ie. when x is an integer plus 0.5) going upwards.
5339 Then half-way cases are identified and adjusted down if the
5340 round-upwards didn't give the desired even integer.
5342 "plus_half == result" identifies a half-way case. If plus_half, which is
5343 x + 0.5, is an integer then x must be an integer plus 0.5.
5345 An odd "result" value is identified with result/2 != floor(result/2).
5346 This is done with plus_half, since that value is ready for use sooner in
5347 a pipelined cpu, and we're already requiring plus_half == result.
5349 Note however that we need to be careful when x is big and already an
5350 integer. In that case "x+0.5" may round to an adjacent integer, causing
5351 us to return such a value, incorrectly. For instance if the hardware is
5352 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5353 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5354 returned. Or if the hardware is in round-upwards mode, then other bigger
5355 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5356 representable value, 2^128+2^76 (or whatever), again incorrect.
5358 These bad roundings of x+0.5 are avoided by testing at the start whether
5359 x is already an integer. If it is then clearly that's the desired result
5360 already. And if it's not then the exponent must be small enough to allow
5361 an 0.5 to be represented, and hence added without a bad rounding. */
5364 scm_c_round (double x
)
5366 double plus_half
, result
;
5371 plus_half
= x
+ 0.5;
5372 result
= floor (plus_half
);
5373 /* Adjust so that the rounding is towards even. */
5374 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5379 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5381 "Round the number @var{x} towards zero.")
5382 #define FUNC_NAME s_scm_truncate_number
5384 if (scm_is_false (scm_negative_p (x
)))
5385 return scm_floor (x
);
5387 return scm_ceiling (x
);
5391 static SCM exactly_one_half
;
5393 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5395 "Round the number @var{x} towards the nearest integer. "
5396 "When it is exactly halfway between two integers, "
5397 "round towards the even one.")
5398 #define FUNC_NAME s_scm_round_number
5400 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5402 else if (SCM_REALP (x
))
5403 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5406 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5407 single quotient+remainder division then examining to see which way
5408 the rounding should go. */
5409 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5410 SCM result
= scm_floor (plus_half
);
5411 /* Adjust so that the rounding is towards even. */
5412 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5413 && scm_is_true (scm_odd_p (result
)))
5414 return scm_difference (result
, SCM_INUM1
);
5421 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5423 "Round the number @var{x} towards minus infinity.")
5424 #define FUNC_NAME s_scm_floor
5426 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5428 else if (SCM_REALP (x
))
5429 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5430 else if (SCM_FRACTIONP (x
))
5432 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5433 SCM_FRACTION_DENOMINATOR (x
));
5434 if (scm_is_false (scm_negative_p (x
)))
5436 /* For positive x, rounding towards zero is correct. */
5441 /* For negative x, we need to return q-1 unless x is an
5442 integer. But fractions are never integer, per our
5444 return scm_difference (q
, SCM_INUM1
);
5448 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5452 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5454 "Round the number @var{x} towards infinity.")
5455 #define FUNC_NAME s_scm_ceiling
5457 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5459 else if (SCM_REALP (x
))
5460 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5461 else if (SCM_FRACTIONP (x
))
5463 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5464 SCM_FRACTION_DENOMINATOR (x
));
5465 if (scm_is_false (scm_positive_p (x
)))
5467 /* For negative x, rounding towards zero is correct. */
5472 /* For positive x, we need to return q+1 unless x is an
5473 integer. But fractions are never integer, per our
5475 return scm_sum (q
, SCM_INUM1
);
5479 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5483 /* sin/cos/tan/asin/acos/atan
5484 sinh/cosh/tanh/asinh/acosh/atanh
5485 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5486 Written by Jerry D. Hedden, (C) FSF.
5487 See the file `COPYING' for terms applying to this program. */
5489 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5491 "Return @var{x} raised to the power of @var{y}.")
5492 #define FUNC_NAME s_scm_expt
5494 if (scm_is_integer (y
))
5496 if (scm_is_true (scm_exact_p (y
)))
5497 return scm_integer_expt (x
, y
);
5500 /* Here we handle the case where the exponent is an inexact
5501 integer. We make the exponent exact in order to use
5502 scm_integer_expt, and thus avoid the spurious imaginary
5503 parts that may result from round-off errors in the general
5504 e^(y log x) method below (for example when squaring a large
5505 negative number). In this case, we must return an inexact
5506 result for correctness. We also make the base inexact so
5507 that scm_integer_expt will use fast inexact arithmetic
5508 internally. Note that making the base inexact is not
5509 sufficient to guarantee an inexact result, because
5510 scm_integer_expt will return an exact 1 when the exponent
5511 is 0, even if the base is inexact. */
5512 return scm_exact_to_inexact
5513 (scm_integer_expt (scm_exact_to_inexact (x
),
5514 scm_inexact_to_exact (y
)));
5517 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5519 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5522 return scm_exp (scm_product (scm_log (x
), y
));
5526 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5528 "Compute the sine of @var{z}.")
5529 #define FUNC_NAME s_scm_sin
5531 if (scm_is_real (z
))
5532 return scm_from_double (sin (scm_to_double (z
)));
5533 else if (SCM_COMPLEXP (z
))
5535 x
= SCM_COMPLEX_REAL (z
);
5536 y
= SCM_COMPLEX_IMAG (z
);
5537 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5538 cos (x
) * sinh (y
));
5541 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5545 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5547 "Compute the cosine of @var{z}.")
5548 #define FUNC_NAME s_scm_cos
5550 if (scm_is_real (z
))
5551 return scm_from_double (cos (scm_to_double (z
)));
5552 else if (SCM_COMPLEXP (z
))
5554 x
= SCM_COMPLEX_REAL (z
);
5555 y
= SCM_COMPLEX_IMAG (z
);
5556 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5557 -sin (x
) * sinh (y
));
5560 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5564 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5566 "Compute the tangent of @var{z}.")
5567 #define FUNC_NAME s_scm_tan
5569 if (scm_is_real (z
))
5570 return scm_from_double (tan (scm_to_double (z
)));
5571 else if (SCM_COMPLEXP (z
))
5573 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5574 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5575 w
= cos (x
) + cosh (y
);
5576 #ifndef ALLOW_DIVIDE_BY_ZERO
5578 scm_num_overflow (s_scm_tan
);
5580 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5583 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5587 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5589 "Compute the hyperbolic sine of @var{z}.")
5590 #define FUNC_NAME s_scm_sinh
5592 if (scm_is_real (z
))
5593 return scm_from_double (sinh (scm_to_double (z
)));
5594 else if (SCM_COMPLEXP (z
))
5596 x
= SCM_COMPLEX_REAL (z
);
5597 y
= SCM_COMPLEX_IMAG (z
);
5598 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5599 cosh (x
) * sin (y
));
5602 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5606 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5608 "Compute the hyperbolic cosine of @var{z}.")
5609 #define FUNC_NAME s_scm_cosh
5611 if (scm_is_real (z
))
5612 return scm_from_double (cosh (scm_to_double (z
)));
5613 else if (SCM_COMPLEXP (z
))
5615 x
= SCM_COMPLEX_REAL (z
);
5616 y
= SCM_COMPLEX_IMAG (z
);
5617 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5618 sinh (x
) * sin (y
));
5621 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5625 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5627 "Compute the hyperbolic tangent of @var{z}.")
5628 #define FUNC_NAME s_scm_tanh
5630 if (scm_is_real (z
))
5631 return scm_from_double (tanh (scm_to_double (z
)));
5632 else if (SCM_COMPLEXP (z
))
5634 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5635 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5636 w
= cosh (x
) + cos (y
);
5637 #ifndef ALLOW_DIVIDE_BY_ZERO
5639 scm_num_overflow (s_scm_tanh
);
5641 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5644 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5648 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5650 "Compute the arc sine of @var{z}.")
5651 #define FUNC_NAME s_scm_asin
5653 if (scm_is_real (z
))
5655 double w
= scm_to_double (z
);
5656 if (w
>= -1.0 && w
<= 1.0)
5657 return scm_from_double (asin (w
));
5659 return scm_product (scm_c_make_rectangular (0, -1),
5660 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5662 else if (SCM_COMPLEXP (z
))
5664 x
= SCM_COMPLEX_REAL (z
);
5665 y
= SCM_COMPLEX_IMAG (z
);
5666 return scm_product (scm_c_make_rectangular (0, -1),
5667 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5670 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5674 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5676 "Compute the arc cosine of @var{z}.")
5677 #define FUNC_NAME s_scm_acos
5679 if (scm_is_real (z
))
5681 double w
= scm_to_double (z
);
5682 if (w
>= -1.0 && w
<= 1.0)
5683 return scm_from_double (acos (w
));
5685 return scm_sum (scm_from_double (acos (0.0)),
5686 scm_product (scm_c_make_rectangular (0, 1),
5687 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5689 else if (SCM_COMPLEXP (z
))
5691 x
= SCM_COMPLEX_REAL (z
);
5692 y
= SCM_COMPLEX_IMAG (z
);
5693 return scm_sum (scm_from_double (acos (0.0)),
5694 scm_product (scm_c_make_rectangular (0, 1),
5695 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5698 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5702 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5704 "With one argument, compute the arc tangent of @var{z}.\n"
5705 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5706 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5707 #define FUNC_NAME s_scm_atan
5711 if (scm_is_real (z
))
5712 return scm_from_double (atan (scm_to_double (z
)));
5713 else if (SCM_COMPLEXP (z
))
5716 v
= SCM_COMPLEX_REAL (z
);
5717 w
= SCM_COMPLEX_IMAG (z
);
5718 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5719 scm_c_make_rectangular (v
, w
+ 1.0))),
5720 scm_c_make_rectangular (0, 2));
5723 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5725 else if (scm_is_real (z
))
5727 if (scm_is_real (y
))
5728 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5730 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5733 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5737 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5739 "Compute the inverse hyperbolic sine of @var{z}.")
5740 #define FUNC_NAME s_scm_sys_asinh
5742 if (scm_is_real (z
))
5743 return scm_from_double (asinh (scm_to_double (z
)));
5744 else if (scm_is_number (z
))
5745 return scm_log (scm_sum (z
,
5746 scm_sqrt (scm_sum (scm_product (z
, z
),
5749 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5753 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5755 "Compute the inverse hyperbolic cosine of @var{z}.")
5756 #define FUNC_NAME s_scm_sys_acosh
5758 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5759 return scm_from_double (acosh (scm_to_double (z
)));
5760 else if (scm_is_number (z
))
5761 return scm_log (scm_sum (z
,
5762 scm_sqrt (scm_difference (scm_product (z
, z
),
5765 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5769 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5771 "Compute the inverse hyperbolic tangent of @var{z}.")
5772 #define FUNC_NAME s_scm_sys_atanh
5774 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5775 return scm_from_double (atanh (scm_to_double (z
)));
5776 else if (scm_is_number (z
))
5777 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
5778 scm_difference (SCM_INUM1
, z
))),
5781 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5786 scm_c_make_rectangular (double re
, double im
)
5789 return scm_from_double (re
);
5794 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5796 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
5797 SCM_COMPLEX_REAL (z
) = re
;
5798 SCM_COMPLEX_IMAG (z
) = im
;
5803 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5804 (SCM real_part
, SCM imaginary_part
),
5805 "Return a complex number constructed of the given @var{real-part} "
5806 "and @var{imaginary-part} parts.")
5807 #define FUNC_NAME s_scm_make_rectangular
5809 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5810 SCM_ARG1
, FUNC_NAME
, "real");
5811 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5812 SCM_ARG2
, FUNC_NAME
, "real");
5813 return scm_c_make_rectangular (scm_to_double (real_part
),
5814 scm_to_double (imaginary_part
));
5819 scm_c_make_polar (double mag
, double ang
)
5823 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5824 use it on Glibc-based systems that have it (it's a GNU extension). See
5825 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5827 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5828 sincos (ang
, &s
, &c
);
5833 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5836 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5838 "Return the complex number @var{x} * e^(i * @var{y}).")
5839 #define FUNC_NAME s_scm_make_polar
5841 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5842 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5843 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5848 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5849 /* "Return the real part of the number @var{z}."
5852 scm_real_part (SCM z
)
5854 if (SCM_I_INUMP (z
))
5856 else if (SCM_BIGP (z
))
5858 else if (SCM_REALP (z
))
5860 else if (SCM_COMPLEXP (z
))
5861 return scm_from_double (SCM_COMPLEX_REAL (z
));
5862 else if (SCM_FRACTIONP (z
))
5865 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5869 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5870 /* "Return the imaginary part of the number @var{z}."
5873 scm_imag_part (SCM z
)
5875 if (SCM_I_INUMP (z
))
5877 else if (SCM_BIGP (z
))
5879 else if (SCM_REALP (z
))
5881 else if (SCM_COMPLEXP (z
))
5882 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5883 else if (SCM_FRACTIONP (z
))
5886 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5889 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5890 /* "Return the numerator of the number @var{z}."
5893 scm_numerator (SCM z
)
5895 if (SCM_I_INUMP (z
))
5897 else if (SCM_BIGP (z
))
5899 else if (SCM_FRACTIONP (z
))
5900 return SCM_FRACTION_NUMERATOR (z
);
5901 else if (SCM_REALP (z
))
5902 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5904 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5908 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5909 /* "Return the denominator of the number @var{z}."
5912 scm_denominator (SCM z
)
5914 if (SCM_I_INUMP (z
))
5916 else if (SCM_BIGP (z
))
5918 else if (SCM_FRACTIONP (z
))
5919 return SCM_FRACTION_DENOMINATOR (z
);
5920 else if (SCM_REALP (z
))
5921 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5923 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5926 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5927 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5928 * "@code{abs} for real arguments, but also allows complex numbers."
5931 scm_magnitude (SCM z
)
5933 if (SCM_I_INUMP (z
))
5935 scm_t_inum zz
= SCM_I_INUM (z
);
5938 else if (SCM_POSFIXABLE (-zz
))
5939 return SCM_I_MAKINUM (-zz
);
5941 return scm_i_inum2big (-zz
);
5943 else if (SCM_BIGP (z
))
5945 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5946 scm_remember_upto_here_1 (z
);
5948 return scm_i_clonebig (z
, 0);
5952 else if (SCM_REALP (z
))
5953 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5954 else if (SCM_COMPLEXP (z
))
5955 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5956 else if (SCM_FRACTIONP (z
))
5958 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5960 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5961 SCM_FRACTION_DENOMINATOR (z
));
5964 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5968 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5969 /* "Return the angle of the complex number @var{z}."
5974 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5975 flo0 to save allocating a new flonum with scm_from_double each time.
5976 But if atan2 follows the floating point rounding mode, then the value
5977 is not a constant. Maybe it'd be close enough though. */
5978 if (SCM_I_INUMP (z
))
5980 if (SCM_I_INUM (z
) >= 0)
5983 return scm_from_double (atan2 (0.0, -1.0));
5985 else if (SCM_BIGP (z
))
5987 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5988 scm_remember_upto_here_1 (z
);
5990 return scm_from_double (atan2 (0.0, -1.0));
5994 else if (SCM_REALP (z
))
5996 if (SCM_REAL_VALUE (z
) >= 0)
5999 return scm_from_double (atan2 (0.0, -1.0));
6001 else if (SCM_COMPLEXP (z
))
6002 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
6003 else if (SCM_FRACTIONP (z
))
6005 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
6007 else return scm_from_double (atan2 (0.0, -1.0));
6010 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
6014 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
6015 /* Convert the number @var{x} to its inexact representation.\n"
6018 scm_exact_to_inexact (SCM z
)
6020 if (SCM_I_INUMP (z
))
6021 return scm_from_double ((double) SCM_I_INUM (z
));
6022 else if (SCM_BIGP (z
))
6023 return scm_from_double (scm_i_big2dbl (z
));
6024 else if (SCM_FRACTIONP (z
))
6025 return scm_from_double (scm_i_fraction2double (z
));
6026 else if (SCM_INEXACTP (z
))
6029 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
6033 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
6035 "Return an exact number that is numerically closest to @var{z}.")
6036 #define FUNC_NAME s_scm_inexact_to_exact
6038 if (SCM_I_INUMP (z
))
6040 else if (SCM_BIGP (z
))
6042 else if (SCM_REALP (z
))
6044 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
6045 SCM_OUT_OF_RANGE (1, z
);
6052 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6053 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6054 scm_i_mpz2num (mpq_denref (frac
)));
6056 /* When scm_i_make_ratio throws, we leak the memory allocated
6063 else if (SCM_FRACTIONP (z
))
6066 SCM_WRONG_TYPE_ARG (1, z
);
6070 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6072 "Returns the @emph{simplest} rational number differing\n"
6073 "from @var{x} by no more than @var{eps}.\n"
6075 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6076 "exact result when both its arguments are exact. Thus, you might need\n"
6077 "to use @code{inexact->exact} on the arguments.\n"
6080 "(rationalize (inexact->exact 1.2) 1/100)\n"
6083 #define FUNC_NAME s_scm_rationalize
6085 if (SCM_I_INUMP (x
))
6087 else if (SCM_BIGP (x
))
6089 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6091 /* Use continued fractions to find closest ratio. All
6092 arithmetic is done with exact numbers.
6095 SCM ex
= scm_inexact_to_exact (x
);
6096 SCM int_part
= scm_floor (ex
);
6098 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
6099 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
6103 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6106 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6107 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6109 /* We stop after a million iterations just to be absolutely sure
6110 that we don't go into an infinite loop. The process normally
6111 converges after less than a dozen iterations.
6114 eps
= scm_abs (eps
);
6115 while (++i
< 1000000)
6117 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6118 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6119 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6121 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6122 eps
))) /* abs(x-a/b) <= eps */
6124 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6125 if (scm_is_false (scm_exact_p (x
))
6126 || scm_is_false (scm_exact_p (eps
)))
6127 return scm_exact_to_inexact (res
);
6131 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6133 tt
= scm_floor (rx
); /* tt = floor (rx) */
6139 scm_num_overflow (s_scm_rationalize
);
6142 SCM_WRONG_TYPE_ARG (1, x
);
6146 /* conversion functions */
6149 scm_is_integer (SCM val
)
6151 return scm_is_true (scm_integer_p (val
));
6155 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6157 if (SCM_I_INUMP (val
))
6159 scm_t_signed_bits n
= SCM_I_INUM (val
);
6160 return n
>= min
&& n
<= max
;
6162 else if (SCM_BIGP (val
))
6164 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6166 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6168 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6170 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6171 return n
>= min
&& n
<= max
;
6181 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6182 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6185 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6186 SCM_I_BIG_MPZ (val
));
6188 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6200 return n
>= min
&& n
<= max
;
6208 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6210 if (SCM_I_INUMP (val
))
6212 scm_t_signed_bits n
= SCM_I_INUM (val
);
6213 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6215 else if (SCM_BIGP (val
))
6217 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6219 else if (max
<= ULONG_MAX
)
6221 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6223 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6224 return n
>= min
&& n
<= max
;
6234 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6237 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6238 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6241 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6242 SCM_I_BIG_MPZ (val
));
6244 return n
>= min
&& n
<= max
;
6252 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6254 scm_error (scm_out_of_range_key
,
6256 "Value out of range ~S to ~S: ~S",
6257 scm_list_3 (min
, max
, bad_val
),
6258 scm_list_1 (bad_val
));
6261 #define TYPE scm_t_intmax
6262 #define TYPE_MIN min
6263 #define TYPE_MAX max
6264 #define SIZEOF_TYPE 0
6265 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6266 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6267 #include "libguile/conv-integer.i.c"
6269 #define TYPE scm_t_uintmax
6270 #define TYPE_MIN min
6271 #define TYPE_MAX max
6272 #define SIZEOF_TYPE 0
6273 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6274 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6275 #include "libguile/conv-uinteger.i.c"
6277 #define TYPE scm_t_int8
6278 #define TYPE_MIN SCM_T_INT8_MIN
6279 #define TYPE_MAX SCM_T_INT8_MAX
6280 #define SIZEOF_TYPE 1
6281 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6282 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6283 #include "libguile/conv-integer.i.c"
6285 #define TYPE scm_t_uint8
6287 #define TYPE_MAX SCM_T_UINT8_MAX
6288 #define SIZEOF_TYPE 1
6289 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6290 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6291 #include "libguile/conv-uinteger.i.c"
6293 #define TYPE scm_t_int16
6294 #define TYPE_MIN SCM_T_INT16_MIN
6295 #define TYPE_MAX SCM_T_INT16_MAX
6296 #define SIZEOF_TYPE 2
6297 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6298 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6299 #include "libguile/conv-integer.i.c"
6301 #define TYPE scm_t_uint16
6303 #define TYPE_MAX SCM_T_UINT16_MAX
6304 #define SIZEOF_TYPE 2
6305 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6306 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6307 #include "libguile/conv-uinteger.i.c"
6309 #define TYPE scm_t_int32
6310 #define TYPE_MIN SCM_T_INT32_MIN
6311 #define TYPE_MAX SCM_T_INT32_MAX
6312 #define SIZEOF_TYPE 4
6313 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6314 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6315 #include "libguile/conv-integer.i.c"
6317 #define TYPE scm_t_uint32
6319 #define TYPE_MAX SCM_T_UINT32_MAX
6320 #define SIZEOF_TYPE 4
6321 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6322 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6323 #include "libguile/conv-uinteger.i.c"
6325 #define TYPE scm_t_wchar
6326 #define TYPE_MIN (scm_t_int32)-1
6327 #define TYPE_MAX (scm_t_int32)0x10ffff
6328 #define SIZEOF_TYPE 4
6329 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6330 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6331 #include "libguile/conv-integer.i.c"
6333 #define TYPE scm_t_int64
6334 #define TYPE_MIN SCM_T_INT64_MIN
6335 #define TYPE_MAX SCM_T_INT64_MAX
6336 #define SIZEOF_TYPE 8
6337 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6338 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6339 #include "libguile/conv-integer.i.c"
6341 #define TYPE scm_t_uint64
6343 #define TYPE_MAX SCM_T_UINT64_MAX
6344 #define SIZEOF_TYPE 8
6345 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6346 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6347 #include "libguile/conv-uinteger.i.c"
6350 scm_to_mpz (SCM val
, mpz_t rop
)
6352 if (SCM_I_INUMP (val
))
6353 mpz_set_si (rop
, SCM_I_INUM (val
));
6354 else if (SCM_BIGP (val
))
6355 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6357 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6361 scm_from_mpz (mpz_t val
)
6363 return scm_i_mpz2num (val
);
6367 scm_is_real (SCM val
)
6369 return scm_is_true (scm_real_p (val
));
6373 scm_is_rational (SCM val
)
6375 return scm_is_true (scm_rational_p (val
));
6379 scm_to_double (SCM val
)
6381 if (SCM_I_INUMP (val
))
6382 return SCM_I_INUM (val
);
6383 else if (SCM_BIGP (val
))
6384 return scm_i_big2dbl (val
);
6385 else if (SCM_FRACTIONP (val
))
6386 return scm_i_fraction2double (val
);
6387 else if (SCM_REALP (val
))
6388 return SCM_REAL_VALUE (val
);
6390 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6394 scm_from_double (double val
)
6398 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
6400 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
6401 SCM_REAL_VALUE (z
) = val
;
6406 #if SCM_ENABLE_DEPRECATED == 1
6409 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
6411 scm_c_issue_deprecation_warning
6412 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6416 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6420 scm_out_of_range (NULL
, num
);
6423 return scm_to_double (num
);
6427 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
6429 scm_c_issue_deprecation_warning
6430 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6434 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6438 scm_out_of_range (NULL
, num
);
6441 return scm_to_double (num
);
6447 scm_is_complex (SCM val
)
6449 return scm_is_true (scm_complex_p (val
));
6453 scm_c_real_part (SCM z
)
6455 if (SCM_COMPLEXP (z
))
6456 return SCM_COMPLEX_REAL (z
);
6459 /* Use the scm_real_part to get proper error checking and
6462 return scm_to_double (scm_real_part (z
));
6467 scm_c_imag_part (SCM z
)
6469 if (SCM_COMPLEXP (z
))
6470 return SCM_COMPLEX_IMAG (z
);
6473 /* Use the scm_imag_part to get proper error checking and
6474 dispatching. The result will almost always be 0.0, but not
6477 return scm_to_double (scm_imag_part (z
));
6482 scm_c_magnitude (SCM z
)
6484 return scm_to_double (scm_magnitude (z
));
6490 return scm_to_double (scm_angle (z
));
6494 scm_is_number (SCM z
)
6496 return scm_is_true (scm_number_p (z
));
6500 /* In the following functions we dispatch to the real-arg funcs like log()
6501 when we know the arg is real, instead of just handing everything to
6502 clog() for instance. This is in case clog() doesn't optimize for a
6503 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6504 well use it to go straight to the applicable C func. */
6506 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6508 "Return the natural logarithm of @var{z}.")
6509 #define FUNC_NAME s_scm_log
6511 if (SCM_COMPLEXP (z
))
6513 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6514 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6516 double re
= SCM_COMPLEX_REAL (z
);
6517 double im
= SCM_COMPLEX_IMAG (z
);
6518 return scm_c_make_rectangular (log (hypot (re
, im
)),
6524 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6525 although the value itself overflows. */
6526 double re
= scm_to_double (z
);
6527 double l
= log (fabs (re
));
6529 return scm_from_double (l
);
6531 return scm_c_make_rectangular (l
, M_PI
);
6537 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6539 "Return the base 10 logarithm of @var{z}.")
6540 #define FUNC_NAME s_scm_log10
6542 if (SCM_COMPLEXP (z
))
6544 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6545 clog() and a multiply by M_LOG10E, rather than the fallback
6546 log10+hypot+atan2.) */
6547 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6548 && defined SCM_COMPLEX_VALUE
6549 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6551 double re
= SCM_COMPLEX_REAL (z
);
6552 double im
= SCM_COMPLEX_IMAG (z
);
6553 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6554 M_LOG10E
* atan2 (im
, re
));
6559 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6560 although the value itself overflows. */
6561 double re
= scm_to_double (z
);
6562 double l
= log10 (fabs (re
));
6564 return scm_from_double (l
);
6566 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6572 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6574 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6575 "base of natural logarithms (2.71828@dots{}).")
6576 #define FUNC_NAME s_scm_exp
6578 if (SCM_COMPLEXP (z
))
6580 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6581 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6583 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6584 SCM_COMPLEX_IMAG (z
));
6589 /* When z is a negative bignum the conversion to double overflows,
6590 giving -infinity, but that's ok, the exp is still 0.0. */
6591 return scm_from_double (exp (scm_to_double (z
)));
6597 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6599 "Return the square root of @var{z}. Of the two possible roots\n"
6600 "(positive and negative), the one with the a positive real part\n"
6601 "is returned, or if that's zero then a positive imaginary part.\n"
6605 "(sqrt 9.0) @result{} 3.0\n"
6606 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6607 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6608 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6610 #define FUNC_NAME s_scm_sqrt
6612 if (SCM_COMPLEXP (x
))
6614 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6615 && defined SCM_COMPLEX_VALUE
6616 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6618 double re
= SCM_COMPLEX_REAL (x
);
6619 double im
= SCM_COMPLEX_IMAG (x
);
6620 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6621 0.5 * atan2 (im
, re
));
6626 double xx
= scm_to_double (x
);
6628 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6630 return scm_from_double (sqrt (xx
));
6642 mpz_init_set_si (z_negative_one
, -1);
6644 /* It may be possible to tune the performance of some algorithms by using
6645 * the following constants to avoid the creation of bignums. Please, before
6646 * using these values, remember the two rules of program optimization:
6647 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6648 scm_c_define ("most-positive-fixnum",
6649 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6650 scm_c_define ("most-negative-fixnum",
6651 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6653 scm_add_feature ("complex");
6654 scm_add_feature ("inexact");
6655 flo0
= scm_from_double (0.0);
6657 /* determine floating point precision */
6658 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6660 init_dblprec(&scm_dblprec
[i
-2],i
);
6661 init_fx_radix(fx_per_radix
[i
-2],i
);
6664 /* hard code precision for base 10 if the preprocessor tells us to... */
6665 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6668 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
6669 #include "libguile/numbers.x"