1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010 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
82 Wonder if this might be faster for some of our code? A switch on
83 the numtag would jump directly to the right case, and the
84 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
86 #define SCM_I_NUMTAG_NOTNUM 0
87 #define SCM_I_NUMTAG_INUM 1
88 #define SCM_I_NUMTAG_BIG scm_tc16_big
89 #define SCM_I_NUMTAG_REAL scm_tc16_real
90 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
91 #define SCM_I_NUMTAG(x) \
92 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
93 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
94 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
95 : SCM_I_NUMTAG_NOTNUM)))
97 /* the macro above will not work as is with fractions */
102 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
104 /* FLOBUFLEN is the maximum number of characters neccessary for the
105 * printed or scm_string representation of an inexact number.
107 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
110 #if !defined (HAVE_ASINH)
111 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
113 #if !defined (HAVE_ACOSH)
114 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
116 #if !defined (HAVE_ATANH)
117 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
120 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
121 an explicit check. In some future gmp (don't know what version number),
122 mpz_cmp_d is supposed to do this itself. */
124 #define xmpz_cmp_d(z, d) \
125 (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
127 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
131 #if defined (GUILE_I)
132 #if HAVE_COMPLEX_DOUBLE
134 /* For an SCM object Z which is a complex number (ie. satisfies
135 SCM_COMPLEXP), return its value as a C level "complex double". */
136 #define SCM_COMPLEX_VALUE(z) \
137 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
139 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
141 /* Convert a C "complex double" to an SCM value. */
143 scm_from_complex_double (complex double z
)
145 return scm_c_make_rectangular (creal (z
), cimag (z
));
148 #endif /* HAVE_COMPLEX_DOUBLE */
153 static mpz_t z_negative_one
;
156 /* Clear the `mpz_t' embedded in bignum PTR. */
158 finalize_bignum (GC_PTR ptr
, GC_PTR data
)
162 bignum
= PTR2SCM (ptr
);
163 mpz_clear (SCM_I_BIG_MPZ (bignum
));
166 /* Return a new uninitialized bignum. */
171 GC_finalization_proc prev_finalizer
;
172 GC_PTR prev_finalizer_data
;
174 /* Allocate one word for the type tag and enough room for an `mpz_t'. */
175 p
= scm_gc_malloc_pointerless (sizeof (scm_t_bits
) + sizeof (mpz_t
),
179 GC_REGISTER_FINALIZER_NO_ORDER (p
, finalize_bignum
, NULL
,
181 &prev_finalizer_data
);
190 /* Return a newly created bignum. */
191 SCM z
= make_bignum ();
192 mpz_init (SCM_I_BIG_MPZ (z
));
197 scm_i_long2big (long x
)
199 /* Return a newly created bignum initialized to X. */
200 SCM z
= make_bignum ();
201 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
206 scm_i_ulong2big (unsigned long x
)
208 /* Return a newly created bignum initialized to X. */
209 SCM z
= make_bignum ();
210 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
215 scm_i_clonebig (SCM src_big
, int same_sign_p
)
217 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
218 SCM z
= make_bignum ();
219 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
221 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
226 scm_i_bigcmp (SCM x
, SCM y
)
228 /* Return neg if x < y, pos if x > y, and 0 if x == y */
229 /* presume we already know x and y are bignums */
230 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
231 scm_remember_upto_here_2 (x
, y
);
236 scm_i_dbl2big (double d
)
238 /* results are only defined if d is an integer */
239 SCM z
= make_bignum ();
240 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
244 /* Convert a integer in double representation to a SCM number. */
247 scm_i_dbl2num (double u
)
249 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
250 powers of 2, so there's no rounding when making "double" values
251 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
252 get rounded on a 64-bit machine, hence the "+1".
254 The use of floor() to force to an integer value ensures we get a
255 "numerically closest" value without depending on how a
256 double->long cast or how mpz_set_d will round. For reference,
257 double->long probably follows the hardware rounding mode,
258 mpz_set_d truncates towards zero. */
260 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
261 representable as a double? */
263 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
264 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
265 return SCM_I_MAKINUM ((long) u
);
267 return scm_i_dbl2big (u
);
270 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
271 with R5RS exact->inexact.
273 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
274 (ie. truncate towards zero), then adjust to get the closest double by
275 examining the next lower bit and adding 1 (to the absolute value) if
278 Bignums exactly half way between representable doubles are rounded to the
279 next higher absolute value (ie. away from zero). This seems like an
280 adequate interpretation of R5RS "numerically closest", and it's easier
281 and faster than a full "nearest-even" style.
283 The bit test must be done on the absolute value of the mpz_t, which means
284 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
285 negatives as twos complement.
287 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
288 following the hardware rounding mode, but applied to the absolute value
289 of the mpz_t operand. This is not what we want so we put the high
290 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
291 mpz_get_d is supposed to always truncate towards zero.
293 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
294 is a slowdown. It'd be faster to pick out the relevant high bits with
295 mpz_getlimbn if we could be bothered coding that, and if the new
296 truncating gmp doesn't come out. */
299 scm_i_big2dbl (SCM b
)
304 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
308 /* Current GMP, eg. 4.1.3, force truncation towards zero */
310 if (bits
> DBL_MANT_DIG
)
312 size_t shift
= bits
- DBL_MANT_DIG
;
313 mpz_init2 (tmp
, DBL_MANT_DIG
);
314 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
315 result
= ldexp (mpz_get_d (tmp
), shift
);
320 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
325 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
328 if (bits
> DBL_MANT_DIG
)
330 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
331 /* test bit number "pos" in absolute value */
332 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
333 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
335 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
339 scm_remember_upto_here_1 (b
);
344 scm_i_normbig (SCM b
)
346 /* convert a big back to a fixnum if it'll fit */
347 /* presume b is a bignum */
348 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
350 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
351 if (SCM_FIXABLE (val
))
352 b
= SCM_I_MAKINUM (val
);
357 static SCM_C_INLINE_KEYWORD SCM
358 scm_i_mpz2num (mpz_t b
)
360 /* convert a mpz number to a SCM number. */
361 if (mpz_fits_slong_p (b
))
363 long val
= mpz_get_si (b
);
364 if (SCM_FIXABLE (val
))
365 return SCM_I_MAKINUM (val
);
369 SCM z
= make_bignum ();
370 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
375 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
376 static SCM
scm_divide2real (SCM x
, SCM y
);
379 scm_i_make_ratio (SCM numerator
, SCM denominator
)
380 #define FUNC_NAME "make-ratio"
382 /* First make sure the arguments are proper.
384 if (SCM_I_INUMP (denominator
))
386 if (scm_is_eq (denominator
, SCM_INUM0
))
387 scm_num_overflow ("make-ratio");
388 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
393 if (!(SCM_BIGP(denominator
)))
394 SCM_WRONG_TYPE_ARG (2, denominator
);
396 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
397 SCM_WRONG_TYPE_ARG (1, numerator
);
399 /* Then flip signs so that the denominator is positive.
401 if (scm_is_true (scm_negative_p (denominator
)))
403 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
404 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
407 /* Now consider for each of the four fixnum/bignum combinations
408 whether the rational number is really an integer.
410 if (SCM_I_INUMP (numerator
))
412 long x
= SCM_I_INUM (numerator
);
413 if (scm_is_eq (numerator
, SCM_INUM0
))
415 if (SCM_I_INUMP (denominator
))
418 y
= SCM_I_INUM (denominator
);
420 return SCM_I_MAKINUM(1);
422 return SCM_I_MAKINUM (x
/ y
);
426 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
427 of that value for the denominator, as a bignum. Apart from
428 that case, abs(bignum) > abs(inum) so inum/bignum is not an
430 if (x
== SCM_MOST_NEGATIVE_FIXNUM
431 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
432 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
433 return SCM_I_MAKINUM(-1);
436 else if (SCM_BIGP (numerator
))
438 if (SCM_I_INUMP (denominator
))
440 long yy
= SCM_I_INUM (denominator
);
441 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
442 return scm_divide (numerator
, denominator
);
446 if (scm_is_eq (numerator
, denominator
))
447 return SCM_I_MAKINUM(1);
448 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
449 SCM_I_BIG_MPZ (denominator
)))
450 return scm_divide(numerator
, denominator
);
454 /* No, it's a proper fraction.
457 SCM divisor
= scm_gcd (numerator
, denominator
);
458 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
460 numerator
= scm_divide (numerator
, divisor
);
461 denominator
= scm_divide (denominator
, divisor
);
464 return scm_double_cell (scm_tc16_fraction
,
465 SCM_UNPACK (numerator
),
466 SCM_UNPACK (denominator
), 0);
472 scm_i_fraction2double (SCM z
)
474 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
475 SCM_FRACTION_DENOMINATOR (z
)));
478 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
480 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
482 #define FUNC_NAME s_scm_exact_p
488 if (SCM_FRACTIONP (x
))
492 SCM_WRONG_TYPE_ARG (1, x
);
497 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
499 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
501 #define FUNC_NAME s_scm_odd_p
505 long val
= SCM_I_INUM (n
);
506 return scm_from_bool ((val
& 1L) != 0);
508 else if (SCM_BIGP (n
))
510 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
511 scm_remember_upto_here_1 (n
);
512 return scm_from_bool (odd_p
);
514 else if (scm_is_true (scm_inf_p (n
)))
516 else if (SCM_REALP (n
))
518 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
524 SCM_WRONG_TYPE_ARG (1, n
);
527 SCM_WRONG_TYPE_ARG (1, n
);
532 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
534 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
536 #define FUNC_NAME s_scm_even_p
540 long val
= SCM_I_INUM (n
);
541 return scm_from_bool ((val
& 1L) == 0);
543 else if (SCM_BIGP (n
))
545 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
546 scm_remember_upto_here_1 (n
);
547 return scm_from_bool (even_p
);
549 else if (scm_is_true (scm_inf_p (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
);
566 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
568 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
569 "or @samp{-inf.0}, @code{#f} otherwise.")
570 #define FUNC_NAME s_scm_inf_p
573 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
574 else if (SCM_COMPLEXP (x
))
575 return scm_from_bool (isinf (SCM_COMPLEX_REAL (x
))
576 || isinf (SCM_COMPLEX_IMAG (x
)));
582 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
584 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
586 #define FUNC_NAME s_scm_nan_p
589 return scm_from_bool (isnan (SCM_REAL_VALUE (n
)));
590 else if (SCM_COMPLEXP (n
))
591 return scm_from_bool (isnan (SCM_COMPLEX_REAL (n
))
592 || isnan (SCM_COMPLEX_IMAG (n
)));
598 /* Guile's idea of infinity. */
599 static double guile_Inf
;
601 /* Guile's idea of not a number. */
602 static double guile_NaN
;
605 guile_ieee_init (void)
607 /* Some version of gcc on some old version of Linux used to crash when
608 trying to make Inf and NaN. */
611 /* C99 INFINITY, when available.
612 FIXME: The standard allows for INFINITY to be something that overflows
613 at compile time. We ought to have a configure test to check for that
614 before trying to use it. (But in practice we believe this is not a
615 problem on any system guile is likely to target.) */
616 guile_Inf
= INFINITY
;
617 #elif defined HAVE_DINFINITY
619 extern unsigned int DINFINITY
[2];
620 guile_Inf
= (*((double *) (DINFINITY
)));
627 if (guile_Inf
== tmp
)
634 /* C99 NAN, when available */
636 #elif defined HAVE_DQNAN
639 extern unsigned int DQNAN
[2];
640 guile_NaN
= (*((double *)(DQNAN
)));
643 guile_NaN
= guile_Inf
/ guile_Inf
;
647 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
650 #define FUNC_NAME s_scm_inf
652 static int initialized
= 0;
658 return scm_from_double (guile_Inf
);
662 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
665 #define FUNC_NAME s_scm_nan
667 static int initialized
= 0;
673 return scm_from_double (guile_NaN
);
678 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
680 "Return the absolute value of @var{x}.")
685 long int xx
= SCM_I_INUM (x
);
688 else if (SCM_POSFIXABLE (-xx
))
689 return SCM_I_MAKINUM (-xx
);
691 return scm_i_long2big (-xx
);
693 else if (SCM_BIGP (x
))
695 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
697 return scm_i_clonebig (x
, 0);
701 else if (SCM_REALP (x
))
703 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
704 double xx
= SCM_REAL_VALUE (x
);
706 return scm_from_double (-xx
);
710 else if (SCM_FRACTIONP (x
))
712 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
714 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
715 SCM_FRACTION_DENOMINATOR (x
));
718 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
723 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
724 /* "Return the quotient of the numbers @var{x} and @var{y}."
727 scm_quotient (SCM x
, SCM y
)
731 long xx
= SCM_I_INUM (x
);
734 long yy
= SCM_I_INUM (y
);
736 scm_num_overflow (s_quotient
);
741 return SCM_I_MAKINUM (z
);
743 return scm_i_long2big (z
);
746 else if (SCM_BIGP (y
))
748 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
749 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
750 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
752 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
753 scm_remember_upto_here_1 (y
);
754 return SCM_I_MAKINUM (-1);
757 return SCM_I_MAKINUM (0);
760 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
762 else if (SCM_BIGP (x
))
766 long yy
= SCM_I_INUM (y
);
768 scm_num_overflow (s_quotient
);
773 SCM result
= scm_i_mkbig ();
776 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
779 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
782 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
783 scm_remember_upto_here_1 (x
);
784 return scm_i_normbig (result
);
787 else if (SCM_BIGP (y
))
789 SCM result
= scm_i_mkbig ();
790 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
793 scm_remember_upto_here_2 (x
, y
);
794 return scm_i_normbig (result
);
797 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
800 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
803 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
804 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
806 * "(remainder 13 4) @result{} 1\n"
807 * "(remainder -13 4) @result{} -1\n"
811 scm_remainder (SCM x
, SCM y
)
817 long yy
= SCM_I_INUM (y
);
819 scm_num_overflow (s_remainder
);
822 long z
= SCM_I_INUM (x
) % yy
;
823 return SCM_I_MAKINUM (z
);
826 else if (SCM_BIGP (y
))
828 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
829 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
830 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
832 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
833 scm_remember_upto_here_1 (y
);
834 return SCM_I_MAKINUM (0);
840 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
842 else if (SCM_BIGP (x
))
846 long yy
= SCM_I_INUM (y
);
848 scm_num_overflow (s_remainder
);
851 SCM result
= scm_i_mkbig ();
854 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
855 scm_remember_upto_here_1 (x
);
856 return scm_i_normbig (result
);
859 else if (SCM_BIGP (y
))
861 SCM result
= scm_i_mkbig ();
862 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
865 scm_remember_upto_here_2 (x
, y
);
866 return scm_i_normbig (result
);
869 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
872 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
876 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
877 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
879 * "(modulo 13 4) @result{} 1\n"
880 * "(modulo -13 4) @result{} 3\n"
884 scm_modulo (SCM x
, SCM y
)
888 long xx
= SCM_I_INUM (x
);
891 long yy
= SCM_I_INUM (y
);
893 scm_num_overflow (s_modulo
);
896 /* C99 specifies that "%" is the remainder corresponding to a
897 quotient rounded towards zero, and that's also traditional
898 for machine division, so z here should be well defined. */
916 return SCM_I_MAKINUM (result
);
919 else if (SCM_BIGP (y
))
921 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
928 SCM pos_y
= scm_i_clonebig (y
, 0);
929 /* do this after the last scm_op */
930 mpz_init_set_si (z_x
, xx
);
931 result
= pos_y
; /* re-use this bignum */
932 mpz_mod (SCM_I_BIG_MPZ (result
),
934 SCM_I_BIG_MPZ (pos_y
));
935 scm_remember_upto_here_1 (pos_y
);
939 result
= scm_i_mkbig ();
940 /* do this after the last scm_op */
941 mpz_init_set_si (z_x
, xx
);
942 mpz_mod (SCM_I_BIG_MPZ (result
),
945 scm_remember_upto_here_1 (y
);
948 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
949 mpz_add (SCM_I_BIG_MPZ (result
),
951 SCM_I_BIG_MPZ (result
));
952 scm_remember_upto_here_1 (y
);
953 /* and do this before the next one */
955 return scm_i_normbig (result
);
959 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
961 else if (SCM_BIGP (x
))
965 long yy
= SCM_I_INUM (y
);
967 scm_num_overflow (s_modulo
);
970 SCM result
= scm_i_mkbig ();
971 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
973 (yy
< 0) ? - yy
: yy
);
974 scm_remember_upto_here_1 (x
);
975 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
976 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
977 SCM_I_BIG_MPZ (result
),
979 return scm_i_normbig (result
);
982 else if (SCM_BIGP (y
))
985 SCM result
= scm_i_mkbig ();
986 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
987 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
988 mpz_mod (SCM_I_BIG_MPZ (result
),
990 SCM_I_BIG_MPZ (pos_y
));
992 scm_remember_upto_here_1 (x
);
993 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
994 mpz_add (SCM_I_BIG_MPZ (result
),
996 SCM_I_BIG_MPZ (result
));
997 scm_remember_upto_here_2 (y
, pos_y
);
998 return scm_i_normbig (result
);
1002 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1005 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1008 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1009 (SCM x
, SCM y
, SCM rest
),
1010 "Return the greatest common divisor of all parameter values.\n"
1011 "If called without arguments, 0 is returned.")
1012 #define FUNC_NAME s_scm_i_gcd
1014 while (!scm_is_null (rest
))
1015 { x
= scm_gcd (x
, y
);
1017 rest
= scm_cdr (rest
);
1019 return scm_gcd (x
, y
);
1023 #define s_gcd s_scm_i_gcd
1024 #define g_gcd g_scm_i_gcd
1027 scm_gcd (SCM x
, SCM y
)
1030 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1032 if (SCM_I_INUMP (x
))
1034 if (SCM_I_INUMP (y
))
1036 long xx
= SCM_I_INUM (x
);
1037 long yy
= SCM_I_INUM (y
);
1038 long u
= xx
< 0 ? -xx
: xx
;
1039 long v
= yy
< 0 ? -yy
: yy
;
1049 /* Determine a common factor 2^k */
1050 while (!(1 & (u
| v
)))
1056 /* Now, any factor 2^n can be eliminated */
1076 return (SCM_POSFIXABLE (result
)
1077 ? SCM_I_MAKINUM (result
)
1078 : scm_i_long2big (result
));
1080 else if (SCM_BIGP (y
))
1086 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1088 else if (SCM_BIGP (x
))
1090 if (SCM_I_INUMP (y
))
1092 unsigned long result
;
1095 yy
= SCM_I_INUM (y
);
1100 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1101 scm_remember_upto_here_1 (x
);
1102 return (SCM_POSFIXABLE (result
)
1103 ? SCM_I_MAKINUM (result
)
1104 : scm_from_ulong (result
));
1106 else if (SCM_BIGP (y
))
1108 SCM result
= scm_i_mkbig ();
1109 mpz_gcd (SCM_I_BIG_MPZ (result
),
1112 scm_remember_upto_here_2 (x
, y
);
1113 return scm_i_normbig (result
);
1116 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1119 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1122 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1123 (SCM x
, SCM y
, SCM rest
),
1124 "Return the least common multiple of the arguments.\n"
1125 "If called without arguments, 1 is returned.")
1126 #define FUNC_NAME s_scm_i_lcm
1128 while (!scm_is_null (rest
))
1129 { x
= scm_lcm (x
, y
);
1131 rest
= scm_cdr (rest
);
1133 return scm_lcm (x
, y
);
1137 #define s_lcm s_scm_i_lcm
1138 #define g_lcm g_scm_i_lcm
1141 scm_lcm (SCM n1
, SCM n2
)
1143 if (SCM_UNBNDP (n2
))
1145 if (SCM_UNBNDP (n1
))
1146 return SCM_I_MAKINUM (1L);
1147 n2
= SCM_I_MAKINUM (1L);
1150 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1151 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1152 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1153 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1155 if (SCM_I_INUMP (n1
))
1157 if (SCM_I_INUMP (n2
))
1159 SCM d
= scm_gcd (n1
, n2
);
1160 if (scm_is_eq (d
, SCM_INUM0
))
1163 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1167 /* inum n1, big n2 */
1170 SCM result
= scm_i_mkbig ();
1171 long nn1
= SCM_I_INUM (n1
);
1172 if (nn1
== 0) return SCM_INUM0
;
1173 if (nn1
< 0) nn1
= - nn1
;
1174 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1175 scm_remember_upto_here_1 (n2
);
1183 if (SCM_I_INUMP (n2
))
1190 SCM result
= scm_i_mkbig ();
1191 mpz_lcm(SCM_I_BIG_MPZ (result
),
1193 SCM_I_BIG_MPZ (n2
));
1194 scm_remember_upto_here_2(n1
, n2
);
1195 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1201 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1206 + + + x (map digit:logand X Y)
1207 + - + x (map digit:logand X (lognot (+ -1 Y)))
1208 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1209 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1214 + + + (map digit:logior X Y)
1215 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1216 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1217 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1222 + + + (map digit:logxor X Y)
1223 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1224 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1225 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1230 + + (any digit:logand X Y)
1231 + - (any digit:logand X (lognot (+ -1 Y)))
1232 - + (any digit:logand (lognot (+ -1 X)) Y)
1237 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1238 (SCM x
, SCM y
, SCM rest
),
1239 "Return the bitwise AND of the integer arguments.\n\n"
1241 "(logand) @result{} -1\n"
1242 "(logand 7) @result{} 7\n"
1243 "(logand #b111 #b011 #b001) @result{} 1\n"
1245 #define FUNC_NAME s_scm_i_logand
1247 while (!scm_is_null (rest
))
1248 { x
= scm_logand (x
, y
);
1250 rest
= scm_cdr (rest
);
1252 return scm_logand (x
, y
);
1256 #define s_scm_logand s_scm_i_logand
1258 SCM
scm_logand (SCM n1
, SCM n2
)
1259 #define FUNC_NAME s_scm_logand
1263 if (SCM_UNBNDP (n2
))
1265 if (SCM_UNBNDP (n1
))
1266 return SCM_I_MAKINUM (-1);
1267 else if (!SCM_NUMBERP (n1
))
1268 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1269 else if (SCM_NUMBERP (n1
))
1272 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1275 if (SCM_I_INUMP (n1
))
1277 nn1
= SCM_I_INUM (n1
);
1278 if (SCM_I_INUMP (n2
))
1280 long nn2
= SCM_I_INUM (n2
);
1281 return SCM_I_MAKINUM (nn1
& nn2
);
1283 else if SCM_BIGP (n2
)
1289 SCM result_z
= scm_i_mkbig ();
1291 mpz_init_set_si (nn1_z
, nn1
);
1292 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1293 scm_remember_upto_here_1 (n2
);
1295 return scm_i_normbig (result_z
);
1299 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1301 else if (SCM_BIGP (n1
))
1303 if (SCM_I_INUMP (n2
))
1306 nn1
= SCM_I_INUM (n1
);
1309 else if (SCM_BIGP (n2
))
1311 SCM result_z
= scm_i_mkbig ();
1312 mpz_and (SCM_I_BIG_MPZ (result_z
),
1314 SCM_I_BIG_MPZ (n2
));
1315 scm_remember_upto_here_2 (n1
, n2
);
1316 return scm_i_normbig (result_z
);
1319 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1322 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1327 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1328 (SCM x
, SCM y
, SCM rest
),
1329 "Return the bitwise OR of the integer arguments.\n\n"
1331 "(logior) @result{} 0\n"
1332 "(logior 7) @result{} 7\n"
1333 "(logior #b000 #b001 #b011) @result{} 3\n"
1335 #define FUNC_NAME s_scm_i_logior
1337 while (!scm_is_null (rest
))
1338 { x
= scm_logior (x
, y
);
1340 rest
= scm_cdr (rest
);
1342 return scm_logior (x
, y
);
1346 #define s_scm_logior s_scm_i_logior
1348 SCM
scm_logior (SCM n1
, SCM n2
)
1349 #define FUNC_NAME s_scm_logior
1353 if (SCM_UNBNDP (n2
))
1355 if (SCM_UNBNDP (n1
))
1357 else if (SCM_NUMBERP (n1
))
1360 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1363 if (SCM_I_INUMP (n1
))
1365 nn1
= SCM_I_INUM (n1
);
1366 if (SCM_I_INUMP (n2
))
1368 long nn2
= SCM_I_INUM (n2
);
1369 return SCM_I_MAKINUM (nn1
| nn2
);
1371 else if (SCM_BIGP (n2
))
1377 SCM result_z
= scm_i_mkbig ();
1379 mpz_init_set_si (nn1_z
, nn1
);
1380 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1381 scm_remember_upto_here_1 (n2
);
1383 return scm_i_normbig (result_z
);
1387 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1389 else if (SCM_BIGP (n1
))
1391 if (SCM_I_INUMP (n2
))
1394 nn1
= SCM_I_INUM (n1
);
1397 else if (SCM_BIGP (n2
))
1399 SCM result_z
= scm_i_mkbig ();
1400 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1402 SCM_I_BIG_MPZ (n2
));
1403 scm_remember_upto_here_2 (n1
, n2
);
1404 return scm_i_normbig (result_z
);
1407 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1410 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1415 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1416 (SCM x
, SCM y
, SCM rest
),
1417 "Return the bitwise XOR of the integer arguments. A bit is\n"
1418 "set in the result if it is set in an odd number of arguments.\n"
1420 "(logxor) @result{} 0\n"
1421 "(logxor 7) @result{} 7\n"
1422 "(logxor #b000 #b001 #b011) @result{} 2\n"
1423 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1425 #define FUNC_NAME s_scm_i_logxor
1427 while (!scm_is_null (rest
))
1428 { x
= scm_logxor (x
, y
);
1430 rest
= scm_cdr (rest
);
1432 return scm_logxor (x
, y
);
1436 #define s_scm_logxor s_scm_i_logxor
1438 SCM
scm_logxor (SCM n1
, SCM n2
)
1439 #define FUNC_NAME s_scm_logxor
1443 if (SCM_UNBNDP (n2
))
1445 if (SCM_UNBNDP (n1
))
1447 else if (SCM_NUMBERP (n1
))
1450 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1453 if (SCM_I_INUMP (n1
))
1455 nn1
= SCM_I_INUM (n1
);
1456 if (SCM_I_INUMP (n2
))
1458 long nn2
= SCM_I_INUM (n2
);
1459 return SCM_I_MAKINUM (nn1
^ nn2
);
1461 else if (SCM_BIGP (n2
))
1465 SCM result_z
= scm_i_mkbig ();
1467 mpz_init_set_si (nn1_z
, nn1
);
1468 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1469 scm_remember_upto_here_1 (n2
);
1471 return scm_i_normbig (result_z
);
1475 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1477 else if (SCM_BIGP (n1
))
1479 if (SCM_I_INUMP (n2
))
1482 nn1
= SCM_I_INUM (n1
);
1485 else if (SCM_BIGP (n2
))
1487 SCM result_z
= scm_i_mkbig ();
1488 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1490 SCM_I_BIG_MPZ (n2
));
1491 scm_remember_upto_here_2 (n1
, n2
);
1492 return scm_i_normbig (result_z
);
1495 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1498 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1503 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1505 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1506 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1507 "without actually calculating the @code{logand}, just testing\n"
1511 "(logtest #b0100 #b1011) @result{} #f\n"
1512 "(logtest #b0100 #b0111) @result{} #t\n"
1514 #define FUNC_NAME s_scm_logtest
1518 if (SCM_I_INUMP (j
))
1520 nj
= SCM_I_INUM (j
);
1521 if (SCM_I_INUMP (k
))
1523 long nk
= SCM_I_INUM (k
);
1524 return scm_from_bool (nj
& nk
);
1526 else if (SCM_BIGP (k
))
1534 mpz_init_set_si (nj_z
, nj
);
1535 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1536 scm_remember_upto_here_1 (k
);
1537 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1543 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1545 else if (SCM_BIGP (j
))
1547 if (SCM_I_INUMP (k
))
1550 nj
= SCM_I_INUM (j
);
1553 else if (SCM_BIGP (k
))
1557 mpz_init (result_z
);
1561 scm_remember_upto_here_2 (j
, k
);
1562 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1563 mpz_clear (result_z
);
1567 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1570 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1575 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1577 "Test whether bit number @var{index} in @var{j} is set.\n"
1578 "@var{index} starts from 0 for the least significant bit.\n"
1581 "(logbit? 0 #b1101) @result{} #t\n"
1582 "(logbit? 1 #b1101) @result{} #f\n"
1583 "(logbit? 2 #b1101) @result{} #t\n"
1584 "(logbit? 3 #b1101) @result{} #t\n"
1585 "(logbit? 4 #b1101) @result{} #f\n"
1587 #define FUNC_NAME s_scm_logbit_p
1589 unsigned long int iindex
;
1590 iindex
= scm_to_ulong (index
);
1592 if (SCM_I_INUMP (j
))
1594 /* bits above what's in an inum follow the sign bit */
1595 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1596 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1598 else if (SCM_BIGP (j
))
1600 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1601 scm_remember_upto_here_1 (j
);
1602 return scm_from_bool (val
);
1605 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1610 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1612 "Return the integer which is the ones-complement of the integer\n"
1616 "(number->string (lognot #b10000000) 2)\n"
1617 " @result{} \"-10000001\"\n"
1618 "(number->string (lognot #b0) 2)\n"
1619 " @result{} \"-1\"\n"
1621 #define FUNC_NAME s_scm_lognot
1623 if (SCM_I_INUMP (n
)) {
1624 /* No overflow here, just need to toggle all the bits making up the inum.
1625 Enhancement: No need to strip the tag and add it back, could just xor
1626 a block of 1 bits, if that worked with the various debug versions of
1628 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1630 } else if (SCM_BIGP (n
)) {
1631 SCM result
= scm_i_mkbig ();
1632 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1633 scm_remember_upto_here_1 (n
);
1637 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1642 /* returns 0 if IN is not an integer. OUT must already be
1645 coerce_to_big (SCM in
, mpz_t out
)
1648 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1649 else if (SCM_I_INUMP (in
))
1650 mpz_set_si (out
, SCM_I_INUM (in
));
1657 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1658 (SCM n
, SCM k
, SCM m
),
1659 "Return @var{n} raised to the integer exponent\n"
1660 "@var{k}, modulo @var{m}.\n"
1663 "(modulo-expt 2 3 5)\n"
1666 #define FUNC_NAME s_scm_modulo_expt
1672 /* There are two classes of error we might encounter --
1673 1) Math errors, which we'll report by calling scm_num_overflow,
1675 2) wrong-type errors, which of course we'll report by calling
1677 We don't report those errors immediately, however; instead we do
1678 some cleanup first. These variables tell us which error (if
1679 any) we should report after cleaning up.
1681 int report_overflow
= 0;
1683 int position_of_wrong_type
= 0;
1684 SCM value_of_wrong_type
= SCM_INUM0
;
1686 SCM result
= SCM_UNDEFINED
;
1692 if (scm_is_eq (m
, SCM_INUM0
))
1694 report_overflow
= 1;
1698 if (!coerce_to_big (n
, n_tmp
))
1700 value_of_wrong_type
= n
;
1701 position_of_wrong_type
= 1;
1705 if (!coerce_to_big (k
, k_tmp
))
1707 value_of_wrong_type
= k
;
1708 position_of_wrong_type
= 2;
1712 if (!coerce_to_big (m
, m_tmp
))
1714 value_of_wrong_type
= m
;
1715 position_of_wrong_type
= 3;
1719 /* if the exponent K is negative, and we simply call mpz_powm, we
1720 will get a divide-by-zero exception when an inverse 1/n mod m
1721 doesn't exist (or is not unique). Since exceptions are hard to
1722 handle, we'll attempt the inversion "by hand" -- that way, we get
1723 a simple failure code, which is easy to handle. */
1725 if (-1 == mpz_sgn (k_tmp
))
1727 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1729 report_overflow
= 1;
1732 mpz_neg (k_tmp
, k_tmp
);
1735 result
= scm_i_mkbig ();
1736 mpz_powm (SCM_I_BIG_MPZ (result
),
1741 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1742 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1749 if (report_overflow
)
1750 scm_num_overflow (FUNC_NAME
);
1752 if (position_of_wrong_type
)
1753 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1754 value_of_wrong_type
);
1756 return scm_i_normbig (result
);
1760 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1762 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1763 "exact integer, @var{n} can be any number.\n"
1765 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1766 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1767 "includes @math{0^0} is 1.\n"
1770 "(integer-expt 2 5) @result{} 32\n"
1771 "(integer-expt -3 3) @result{} -27\n"
1772 "(integer-expt 5 -3) @result{} 1/125\n"
1773 "(integer-expt 0 0) @result{} 1\n"
1775 #define FUNC_NAME s_scm_integer_expt
1778 SCM z_i2
= SCM_BOOL_F
;
1780 SCM acc
= SCM_I_MAKINUM (1L);
1782 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1784 /* 0^0 == 1 according to R5RS */
1785 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1786 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1787 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1788 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1790 if (SCM_I_INUMP (k
))
1791 i2
= SCM_I_INUM (k
);
1792 else if (SCM_BIGP (k
))
1794 z_i2
= scm_i_clonebig (k
, 1);
1795 scm_remember_upto_here_1 (k
);
1799 SCM_WRONG_TYPE_ARG (2, k
);
1803 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1805 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1806 n
= scm_divide (n
, SCM_UNDEFINED
);
1810 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1814 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1816 return scm_product (acc
, n
);
1818 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1819 acc
= scm_product (acc
, n
);
1820 n
= scm_product (n
, n
);
1821 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1829 n
= scm_divide (n
, SCM_UNDEFINED
);
1836 return scm_product (acc
, n
);
1838 acc
= scm_product (acc
, n
);
1839 n
= scm_product (n
, n
);
1846 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1848 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1849 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1851 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1852 "@var{cnt} is negative it's a division, rounded towards negative\n"
1853 "infinity. (Note that this is not the same rounding as\n"
1854 "@code{quotient} does.)\n"
1856 "With @var{n} viewed as an infinite precision twos complement,\n"
1857 "@code{ash} means a left shift introducing zero bits, or a right\n"
1858 "shift dropping bits.\n"
1861 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1862 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1864 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1865 "(ash -23 -2) @result{} -6\n"
1867 #define FUNC_NAME s_scm_ash
1870 bits_to_shift
= scm_to_long (cnt
);
1872 if (SCM_I_INUMP (n
))
1874 long nn
= SCM_I_INUM (n
);
1876 if (bits_to_shift
> 0)
1878 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1879 overflow a non-zero fixnum. For smaller shifts we check the
1880 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1881 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1882 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1888 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1890 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1893 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1897 SCM result
= scm_i_long2big (nn
);
1898 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1905 bits_to_shift
= -bits_to_shift
;
1906 if (bits_to_shift
>= SCM_LONG_BIT
)
1907 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1909 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1913 else if (SCM_BIGP (n
))
1917 if (bits_to_shift
== 0)
1920 result
= scm_i_mkbig ();
1921 if (bits_to_shift
>= 0)
1923 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1929 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1930 we have to allocate a bignum even if the result is going to be a
1932 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1934 return scm_i_normbig (result
);
1940 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1946 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1947 (SCM n
, SCM start
, SCM end
),
1948 "Return the integer composed of the @var{start} (inclusive)\n"
1949 "through @var{end} (exclusive) bits of @var{n}. The\n"
1950 "@var{start}th bit becomes the 0-th bit in the result.\n"
1953 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1954 " @result{} \"1010\"\n"
1955 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1956 " @result{} \"10110\"\n"
1958 #define FUNC_NAME s_scm_bit_extract
1960 unsigned long int istart
, iend
, bits
;
1961 istart
= scm_to_ulong (start
);
1962 iend
= scm_to_ulong (end
);
1963 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1965 /* how many bits to keep */
1966 bits
= iend
- istart
;
1968 if (SCM_I_INUMP (n
))
1970 long int in
= SCM_I_INUM (n
);
1972 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1973 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1974 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1976 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1978 /* Since we emulate two's complement encoded numbers, this
1979 * special case requires us to produce a result that has
1980 * more bits than can be stored in a fixnum.
1982 SCM result
= scm_i_long2big (in
);
1983 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1988 /* mask down to requisite bits */
1989 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1990 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1992 else if (SCM_BIGP (n
))
1997 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2001 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2002 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2003 such bits into a ulong. */
2004 result
= scm_i_mkbig ();
2005 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2006 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2007 result
= scm_i_normbig (result
);
2009 scm_remember_upto_here_1 (n
);
2013 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2018 static const char scm_logtab
[] = {
2019 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2022 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2024 "Return the number of bits in integer @var{n}. If integer is\n"
2025 "positive, the 1-bits in its binary representation are counted.\n"
2026 "If negative, the 0-bits in its two's-complement binary\n"
2027 "representation are counted. If 0, 0 is returned.\n"
2030 "(logcount #b10101010)\n"
2037 #define FUNC_NAME s_scm_logcount
2039 if (SCM_I_INUMP (n
))
2041 unsigned long int c
= 0;
2042 long int nn
= SCM_I_INUM (n
);
2047 c
+= scm_logtab
[15 & nn
];
2050 return SCM_I_MAKINUM (c
);
2052 else if (SCM_BIGP (n
))
2054 unsigned long count
;
2055 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2056 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2058 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2059 scm_remember_upto_here_1 (n
);
2060 return SCM_I_MAKINUM (count
);
2063 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2068 static const char scm_ilentab
[] = {
2069 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2073 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2075 "Return the number of bits necessary to represent @var{n}.\n"
2078 "(integer-length #b10101010)\n"
2080 "(integer-length 0)\n"
2082 "(integer-length #b1111)\n"
2085 #define FUNC_NAME s_scm_integer_length
2087 if (SCM_I_INUMP (n
))
2089 unsigned long int c
= 0;
2091 long int nn
= SCM_I_INUM (n
);
2097 l
= scm_ilentab
[15 & nn
];
2100 return SCM_I_MAKINUM (c
- 4 + l
);
2102 else if (SCM_BIGP (n
))
2104 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2105 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2106 1 too big, so check for that and adjust. */
2107 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2108 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2109 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2110 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2112 scm_remember_upto_here_1 (n
);
2113 return SCM_I_MAKINUM (size
);
2116 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2120 /*** NUMBERS -> STRINGS ***/
2121 #define SCM_MAX_DBL_PREC 60
2122 #define SCM_MAX_DBL_RADIX 36
2124 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2125 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2126 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2129 void init_dblprec(int *prec
, int radix
) {
2130 /* determine floating point precision by adding successively
2131 smaller increments to 1.0 until it is considered == 1.0 */
2132 double f
= ((double)1.0)/radix
;
2133 double fsum
= 1.0 + f
;
2138 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2150 void init_fx_radix(double *fx_list
, int radix
)
2152 /* initialize a per-radix list of tolerances. When added
2153 to a number < 1.0, we can determine if we should raund
2154 up and quit converting a number to a string. */
2158 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2159 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2162 /* use this array as a way to generate a single digit */
2163 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2166 idbl2str (double f
, char *a
, int radix
)
2168 int efmt
, dpt
, d
, i
, wp
;
2170 #ifdef DBL_MIN_10_EXP
2173 #endif /* DBL_MIN_10_EXP */
2178 radix
> SCM_MAX_DBL_RADIX
)
2180 /* revert to existing behavior */
2184 wp
= scm_dblprec
[radix
-2];
2185 fx
= fx_per_radix
[radix
-2];
2189 #ifdef HAVE_COPYSIGN
2190 double sgn
= copysign (1.0, f
);
2195 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2201 strcpy (a
, "-inf.0");
2203 strcpy (a
, "+inf.0");
2208 strcpy (a
, "+nan.0");
2218 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2219 make-uniform-vector, from causing infinite loops. */
2220 /* just do the checking...if it passes, we do the conversion for our
2221 radix again below */
2228 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2236 while (f_cpy
> 10.0)
2239 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2260 if (f
+ fx
[wp
] >= radix
)
2267 /* adding 9999 makes this equivalent to abs(x) % 3 */
2268 dpt
= (exp
+ 9999) % 3;
2272 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2294 a
[ch
++] = number_chars
[d
];
2297 if (f
+ fx
[wp
] >= 1.0)
2299 a
[ch
- 1] = number_chars
[d
+1];
2311 if ((dpt
> 4) && (exp
> 6))
2313 d
= (a
[0] == '-' ? 2 : 1);
2314 for (i
= ch
++; i
> d
; i
--)
2327 if (a
[ch
- 1] == '.')
2328 a
[ch
++] = '0'; /* trailing zero */
2337 for (i
= radix
; i
<= exp
; i
*= radix
);
2338 for (i
/= radix
; i
; i
/= radix
)
2340 a
[ch
++] = number_chars
[exp
/ i
];
2349 icmplx2str (double real
, double imag
, char *str
, int radix
)
2353 i
= idbl2str (real
, str
, radix
);
2356 /* Don't output a '+' for negative numbers or for Inf and
2357 NaN. They will provide their own sign. */
2358 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
2360 i
+= idbl2str (imag
, &str
[i
], radix
);
2367 iflo2str (SCM flt
, char *str
, int radix
)
2370 if (SCM_REALP (flt
))
2371 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2373 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2378 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2379 characters in the result.
2381 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2383 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2388 return scm_iuint2str (-num
, rad
, p
) + 1;
2391 return scm_iuint2str (num
, rad
, p
);
2394 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2395 characters in the result.
2397 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2399 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2403 scm_t_uintmax n
= num
;
2405 for (n
/= rad
; n
> 0; n
/= rad
)
2415 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2420 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2422 "Return a string holding the external representation of the\n"
2423 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2424 "inexact, a radix of 10 will be used.")
2425 #define FUNC_NAME s_scm_number_to_string
2429 if (SCM_UNBNDP (radix
))
2432 base
= scm_to_signed_integer (radix
, 2, 36);
2434 if (SCM_I_INUMP (n
))
2436 char num_buf
[SCM_INTBUFLEN
];
2437 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2438 return scm_from_locale_stringn (num_buf
, length
);
2440 else if (SCM_BIGP (n
))
2442 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2443 scm_remember_upto_here_1 (n
);
2444 return scm_take_locale_string (str
);
2446 else if (SCM_FRACTIONP (n
))
2448 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2449 scm_from_locale_string ("/"),
2450 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2452 else if (SCM_INEXACTP (n
))
2454 char num_buf
[FLOBUFLEN
];
2455 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2458 SCM_WRONG_TYPE_ARG (1, n
);
2463 /* These print routines used to be stubbed here so that scm_repl.c
2464 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2467 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2469 char num_buf
[FLOBUFLEN
];
2470 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2475 scm_i_print_double (double val
, SCM port
)
2477 char num_buf
[FLOBUFLEN
];
2478 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2482 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2485 char num_buf
[FLOBUFLEN
];
2486 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2491 scm_i_print_complex (double real
, double imag
, SCM port
)
2493 char num_buf
[FLOBUFLEN
];
2494 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2498 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2501 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2502 scm_lfwrite_str (str
, port
);
2503 scm_remember_upto_here_1 (str
);
2508 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2510 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2511 scm_remember_upto_here_1 (exp
);
2512 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2516 /*** END nums->strs ***/
2519 /*** STRINGS -> NUMBERS ***/
2521 /* The following functions implement the conversion from strings to numbers.
2522 * The implementation somehow follows the grammar for numbers as it is given
2523 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2524 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2525 * points should be noted about the implementation:
2526 * * Each function keeps a local index variable 'idx' that points at the
2527 * current position within the parsed string. The global index is only
2528 * updated if the function could parse the corresponding syntactic unit
2530 * * Similarly, the functions keep track of indicators of inexactness ('#',
2531 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2532 * global exactness information is only updated after each part has been
2533 * successfully parsed.
2534 * * Sequences of digits are parsed into temporary variables holding fixnums.
2535 * Only if these fixnums would overflow, the result variables are updated
2536 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2537 * the temporary variables holding the fixnums are cleared, and the process
2538 * starts over again. If for example fixnums were able to store five decimal
2539 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2540 * and the result was computed as 12345 * 100000 + 67890. In other words,
2541 * only every five digits two bignum operations were performed.
2544 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2546 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2548 /* In non ASCII-style encodings the following macro might not work. */
2549 #define XDIGIT2UINT(d) \
2550 (uc_is_property_decimal_digit ((int) (unsigned char) d) \
2552 : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
2555 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2556 unsigned int radix
, enum t_exactness
*p_exactness
)
2558 unsigned int idx
= *p_idx
;
2559 unsigned int hash_seen
= 0;
2560 scm_t_bits shift
= 1;
2562 unsigned int digit_value
;
2565 size_t len
= scm_i_string_length (mem
);
2570 c
= scm_i_string_ref (mem
, idx
);
2571 if (!uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2573 digit_value
= XDIGIT2UINT (c
);
2574 if (digit_value
>= radix
)
2578 result
= SCM_I_MAKINUM (digit_value
);
2581 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2582 if (uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2586 digit_value
= XDIGIT2UINT (c
);
2587 if (digit_value
>= radix
)
2599 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2601 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2603 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2610 shift
= shift
* radix
;
2611 add
= add
* radix
+ digit_value
;
2616 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2618 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2622 *p_exactness
= INEXACT
;
2628 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2629 * covers the parts of the rules that start at a potential point. The value
2630 * of the digits up to the point have been parsed by the caller and are given
2631 * in variable result. The content of *p_exactness indicates, whether a hash
2632 * has already been seen in the digits before the point.
2635 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2638 mem2decimal_from_point (SCM result
, SCM mem
,
2639 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2641 unsigned int idx
= *p_idx
;
2642 enum t_exactness x
= *p_exactness
;
2643 size_t len
= scm_i_string_length (mem
);
2648 if (scm_i_string_ref (mem
, idx
) == '.')
2650 scm_t_bits shift
= 1;
2652 unsigned int digit_value
;
2653 SCM big_shift
= SCM_I_MAKINUM (1);
2658 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2659 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2664 digit_value
= DIGIT2UINT (c
);
2675 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2677 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2678 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2680 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2688 add
= add
* 10 + digit_value
;
2694 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2695 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2696 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2699 result
= scm_divide (result
, big_shift
);
2701 /* We've seen a decimal point, thus the value is implicitly inexact. */
2713 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2715 switch (scm_i_string_ref (mem
, idx
))
2727 c
= scm_i_string_ref (mem
, idx
);
2735 c
= scm_i_string_ref (mem
, idx
);
2744 c
= scm_i_string_ref (mem
, idx
);
2749 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2753 exponent
= DIGIT2UINT (c
);
2756 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2757 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2760 if (exponent
<= SCM_MAXEXP
)
2761 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2767 if (exponent
> SCM_MAXEXP
)
2769 size_t exp_len
= idx
- start
;
2770 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2771 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2772 scm_out_of_range ("string->number", exp_num
);
2775 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2777 result
= scm_product (result
, e
);
2779 result
= scm_divide2real (result
, e
);
2781 /* We've seen an exponent, thus the value is implicitly inexact. */
2799 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2802 mem2ureal (SCM mem
, unsigned int *p_idx
,
2803 unsigned int radix
, enum t_exactness
*p_exactness
)
2805 unsigned int idx
= *p_idx
;
2807 size_t len
= scm_i_string_length (mem
);
2809 /* Start off believing that the number will be exact. This changes
2810 to INEXACT if we see a decimal point or a hash. */
2811 enum t_exactness x
= EXACT
;
2816 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2822 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2824 /* Cobble up the fractional part. We might want to set the
2825 NaN's mantissa from it. */
2827 mem2uinteger (mem
, &idx
, 10, &x
);
2832 if (scm_i_string_ref (mem
, idx
) == '.')
2836 else if (idx
+ 1 == len
)
2838 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2841 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2848 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2849 if (scm_is_false (uinteger
))
2854 else if (scm_i_string_ref (mem
, idx
) == '/')
2862 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2863 if (scm_is_false (divisor
))
2866 /* both are int/big here, I assume */
2867 result
= scm_i_make_ratio (uinteger
, divisor
);
2869 else if (radix
== 10)
2871 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2872 if (scm_is_false (result
))
2881 /* Update *p_exactness if the number just read was inexact. This is
2882 important for complex numbers, so that a complex number is
2883 treated as inexact overall if either its real or imaginary part
2889 /* When returning an inexact zero, make sure it is represented as a
2890 floating point value so that we can change its sign.
2892 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2893 result
= scm_from_double (0.0);
2899 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2902 mem2complex (SCM mem
, unsigned int idx
,
2903 unsigned int radix
, enum t_exactness
*p_exactness
)
2908 size_t len
= scm_i_string_length (mem
);
2913 c
= scm_i_string_ref (mem
, idx
);
2928 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2929 if (scm_is_false (ureal
))
2931 /* input must be either +i or -i */
2936 if (scm_i_string_ref (mem
, idx
) == 'i'
2937 || scm_i_string_ref (mem
, idx
) == 'I')
2943 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2950 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2951 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2956 c
= scm_i_string_ref (mem
, idx
);
2960 /* either +<ureal>i or -<ureal>i */
2967 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2970 /* polar input: <real>@<real>. */
2981 c
= scm_i_string_ref (mem
, idx
);
2999 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3000 if (scm_is_false (angle
))
3005 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3006 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3008 result
= scm_make_polar (ureal
, angle
);
3013 /* expecting input matching <real>[+-]<ureal>?i */
3020 int sign
= (c
== '+') ? 1 : -1;
3021 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3023 if (scm_is_false (imag
))
3024 imag
= SCM_I_MAKINUM (sign
);
3025 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
3026 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3030 if (scm_i_string_ref (mem
, idx
) != 'i'
3031 && scm_i_string_ref (mem
, idx
) != 'I')
3038 return scm_make_rectangular (ureal
, imag
);
3047 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3049 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3052 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3054 unsigned int idx
= 0;
3055 unsigned int radix
= NO_RADIX
;
3056 enum t_exactness forced_x
= NO_EXACTNESS
;
3057 enum t_exactness implicit_x
= EXACT
;
3059 size_t len
= scm_i_string_length (mem
);
3061 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3062 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3064 switch (scm_i_string_ref (mem
, idx
+ 1))
3067 if (radix
!= NO_RADIX
)
3072 if (radix
!= NO_RADIX
)
3077 if (forced_x
!= NO_EXACTNESS
)
3082 if (forced_x
!= NO_EXACTNESS
)
3087 if (radix
!= NO_RADIX
)
3092 if (radix
!= NO_RADIX
)
3102 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3103 if (radix
== NO_RADIX
)
3104 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3106 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3108 if (scm_is_false (result
))
3114 if (SCM_INEXACTP (result
))
3115 return scm_inexact_to_exact (result
);
3119 if (SCM_INEXACTP (result
))
3122 return scm_exact_to_inexact (result
);
3125 if (implicit_x
== INEXACT
)
3127 if (SCM_INEXACTP (result
))
3130 return scm_exact_to_inexact (result
);
3138 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3139 unsigned int default_radix
)
3141 SCM str
= scm_from_locale_stringn (mem
, len
);
3143 return scm_i_string_to_number (str
, default_radix
);
3147 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3148 (SCM string
, SCM radix
),
3149 "Return a number of the maximally precise representation\n"
3150 "expressed by the given @var{string}. @var{radix} must be an\n"
3151 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3152 "is a default radix that may be overridden by an explicit radix\n"
3153 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3154 "supplied, then the default radix is 10. If string is not a\n"
3155 "syntactically valid notation for a number, then\n"
3156 "@code{string->number} returns @code{#f}.")
3157 #define FUNC_NAME s_scm_string_to_number
3161 SCM_VALIDATE_STRING (1, string
);
3163 if (SCM_UNBNDP (radix
))
3166 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3168 answer
= scm_i_string_to_number (string
, base
);
3169 scm_remember_upto_here_1 (string
);
3175 /*** END strs->nums ***/
3179 scm_bigequal (SCM x
, SCM y
)
3181 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3182 scm_remember_upto_here_2 (x
, y
);
3183 return scm_from_bool (0 == result
);
3187 scm_real_equalp (SCM x
, SCM y
)
3189 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3193 scm_complex_equalp (SCM x
, SCM y
)
3195 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3196 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3200 scm_i_fraction_equalp (SCM x
, SCM y
)
3202 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3203 SCM_FRACTION_NUMERATOR (y
)))
3204 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3205 SCM_FRACTION_DENOMINATOR (y
))))
3212 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3214 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3216 #define FUNC_NAME s_scm_number_p
3218 return scm_from_bool (SCM_NUMBERP (x
));
3222 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3224 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3225 "otherwise. Note that the sets of real, rational and integer\n"
3226 "values form subsets of the set of complex numbers, i. e. the\n"
3227 "predicate will also be fulfilled if @var{x} is a real,\n"
3228 "rational or integer number.")
3229 #define FUNC_NAME s_scm_complex_p
3231 /* all numbers are complex. */
3232 return scm_number_p (x
);
3236 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3238 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3239 "otherwise. Note that the set of integer values forms a subset of\n"
3240 "the set of real numbers, i. e. the predicate will also be\n"
3241 "fulfilled if @var{x} is an integer number.")
3242 #define FUNC_NAME s_scm_real_p
3244 /* we can't represent irrational numbers. */
3245 return scm_rational_p (x
);
3249 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3251 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3252 "otherwise. Note that the set of integer values forms a subset of\n"
3253 "the set of rational numbers, i. e. the predicate will also be\n"
3254 "fulfilled if @var{x} is an integer number.")
3255 #define FUNC_NAME s_scm_rational_p
3257 if (SCM_I_INUMP (x
))
3259 else if (SCM_IMP (x
))
3261 else if (SCM_BIGP (x
))
3263 else if (SCM_FRACTIONP (x
))
3265 else if (SCM_REALP (x
))
3266 /* due to their limited precision, all floating point numbers are
3267 rational as well. */
3274 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3276 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3278 #define FUNC_NAME s_scm_integer_p
3281 if (SCM_I_INUMP (x
))
3287 if (!SCM_INEXACTP (x
))
3289 if (SCM_COMPLEXP (x
))
3291 r
= SCM_REAL_VALUE (x
);
3292 /* +/-inf passes r==floor(r), making those #t */
3300 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3302 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3304 #define FUNC_NAME s_scm_inexact_p
3306 if (SCM_INEXACTP (x
))
3308 if (SCM_NUMBERP (x
))
3310 SCM_WRONG_TYPE_ARG (1, x
);
3315 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3316 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3317 (SCM x
, SCM y
, SCM rest
),
3318 "Return @code{#t} if all parameters are numerically equal.")
3319 #define FUNC_NAME s_scm_i_num_eq_p
3321 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3323 while (!scm_is_null (rest
))
3325 if (scm_is_false (scm_num_eq_p (x
, y
)))
3329 rest
= scm_cdr (rest
);
3331 return scm_num_eq_p (x
, y
);
3335 scm_num_eq_p (SCM x
, SCM y
)
3338 if (SCM_I_INUMP (x
))
3340 long xx
= SCM_I_INUM (x
);
3341 if (SCM_I_INUMP (y
))
3343 long yy
= SCM_I_INUM (y
);
3344 return scm_from_bool (xx
== yy
);
3346 else if (SCM_BIGP (y
))
3348 else if (SCM_REALP (y
))
3350 /* On a 32-bit system an inum fits a double, we can cast the inum
3351 to a double and compare.
3353 But on a 64-bit system an inum is bigger than a double and
3354 casting it to a double (call that dxx) will round. dxx is at
3355 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3356 an integer and fits a long. So we cast yy to a long and
3357 compare with plain xx.
3359 An alternative (for any size system actually) would be to check
3360 yy is an integer (with floor) and is in range of an inum
3361 (compare against appropriate powers of 2) then test
3362 xx==(long)yy. It's just a matter of which casts/comparisons
3363 might be fastest or easiest for the cpu. */
3365 double yy
= SCM_REAL_VALUE (y
);
3366 return scm_from_bool ((double) xx
== yy
3367 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3368 || xx
== (long) yy
));
3370 else if (SCM_COMPLEXP (y
))
3371 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3372 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3373 else if (SCM_FRACTIONP (y
))
3376 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3378 else if (SCM_BIGP (x
))
3380 if (SCM_I_INUMP (y
))
3382 else if (SCM_BIGP (y
))
3384 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3385 scm_remember_upto_here_2 (x
, y
);
3386 return scm_from_bool (0 == cmp
);
3388 else if (SCM_REALP (y
))
3391 if (isnan (SCM_REAL_VALUE (y
)))
3393 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3394 scm_remember_upto_here_1 (x
);
3395 return scm_from_bool (0 == cmp
);
3397 else if (SCM_COMPLEXP (y
))
3400 if (0.0 != SCM_COMPLEX_IMAG (y
))
3402 if (isnan (SCM_COMPLEX_REAL (y
)))
3404 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3405 scm_remember_upto_here_1 (x
);
3406 return scm_from_bool (0 == cmp
);
3408 else if (SCM_FRACTIONP (y
))
3411 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3413 else if (SCM_REALP (x
))
3415 double xx
= SCM_REAL_VALUE (x
);
3416 if (SCM_I_INUMP (y
))
3418 /* see comments with inum/real above */
3419 long yy
= SCM_I_INUM (y
);
3420 return scm_from_bool (xx
== (double) yy
3421 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3422 || (long) xx
== yy
));
3424 else if (SCM_BIGP (y
))
3427 if (isnan (SCM_REAL_VALUE (x
)))
3429 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3430 scm_remember_upto_here_1 (y
);
3431 return scm_from_bool (0 == cmp
);
3433 else if (SCM_REALP (y
))
3434 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3435 else if (SCM_COMPLEXP (y
))
3436 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3437 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3438 else if (SCM_FRACTIONP (y
))
3440 double xx
= SCM_REAL_VALUE (x
);
3444 return scm_from_bool (xx
< 0.0);
3445 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3449 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3451 else if (SCM_COMPLEXP (x
))
3453 if (SCM_I_INUMP (y
))
3454 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3455 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3456 else if (SCM_BIGP (y
))
3459 if (0.0 != SCM_COMPLEX_IMAG (x
))
3461 if (isnan (SCM_COMPLEX_REAL (x
)))
3463 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3464 scm_remember_upto_here_1 (y
);
3465 return scm_from_bool (0 == cmp
);
3467 else if (SCM_REALP (y
))
3468 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3469 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3470 else if (SCM_COMPLEXP (y
))
3471 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3472 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3473 else if (SCM_FRACTIONP (y
))
3476 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3478 xx
= SCM_COMPLEX_REAL (x
);
3482 return scm_from_bool (xx
< 0.0);
3483 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3487 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3489 else if (SCM_FRACTIONP (x
))
3491 if (SCM_I_INUMP (y
))
3493 else if (SCM_BIGP (y
))
3495 else if (SCM_REALP (y
))
3497 double yy
= SCM_REAL_VALUE (y
);
3501 return scm_from_bool (0.0 < yy
);
3502 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3505 else if (SCM_COMPLEXP (y
))
3508 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3510 yy
= SCM_COMPLEX_REAL (y
);
3514 return scm_from_bool (0.0 < yy
);
3515 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3518 else if (SCM_FRACTIONP (y
))
3519 return scm_i_fraction_equalp (x
, y
);
3521 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3524 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3528 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3529 done are good for inums, but for bignums an answer can almost always be
3530 had by just examining a few high bits of the operands, as done by GMP in
3531 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3532 of the float exponent to take into account. */
3534 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3535 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3536 (SCM x
, SCM y
, SCM rest
),
3537 "Return @code{#t} if the list of parameters is monotonically\n"
3539 #define FUNC_NAME s_scm_i_num_less_p
3541 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3543 while (!scm_is_null (rest
))
3545 if (scm_is_false (scm_less_p (x
, y
)))
3549 rest
= scm_cdr (rest
);
3551 return scm_less_p (x
, y
);
3555 scm_less_p (SCM x
, SCM y
)
3558 if (SCM_I_INUMP (x
))
3560 long xx
= SCM_I_INUM (x
);
3561 if (SCM_I_INUMP (y
))
3563 long yy
= SCM_I_INUM (y
);
3564 return scm_from_bool (xx
< yy
);
3566 else if (SCM_BIGP (y
))
3568 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3569 scm_remember_upto_here_1 (y
);
3570 return scm_from_bool (sgn
> 0);
3572 else if (SCM_REALP (y
))
3573 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3574 else if (SCM_FRACTIONP (y
))
3576 /* "x < a/b" becomes "x*b < a" */
3578 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3579 y
= SCM_FRACTION_NUMERATOR (y
);
3583 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3585 else if (SCM_BIGP (x
))
3587 if (SCM_I_INUMP (y
))
3589 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3590 scm_remember_upto_here_1 (x
);
3591 return scm_from_bool (sgn
< 0);
3593 else if (SCM_BIGP (y
))
3595 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3596 scm_remember_upto_here_2 (x
, y
);
3597 return scm_from_bool (cmp
< 0);
3599 else if (SCM_REALP (y
))
3602 if (isnan (SCM_REAL_VALUE (y
)))
3604 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3605 scm_remember_upto_here_1 (x
);
3606 return scm_from_bool (cmp
< 0);
3608 else if (SCM_FRACTIONP (y
))
3611 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3613 else if (SCM_REALP (x
))
3615 if (SCM_I_INUMP (y
))
3616 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3617 else if (SCM_BIGP (y
))
3620 if (isnan (SCM_REAL_VALUE (x
)))
3622 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3623 scm_remember_upto_here_1 (y
);
3624 return scm_from_bool (cmp
> 0);
3626 else if (SCM_REALP (y
))
3627 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3628 else if (SCM_FRACTIONP (y
))
3630 double xx
= SCM_REAL_VALUE (x
);
3634 return scm_from_bool (xx
< 0.0);
3635 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3639 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3641 else if (SCM_FRACTIONP (x
))
3643 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3645 /* "a/b < y" becomes "a < y*b" */
3646 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3647 x
= SCM_FRACTION_NUMERATOR (x
);
3650 else if (SCM_REALP (y
))
3652 double yy
= SCM_REAL_VALUE (y
);
3656 return scm_from_bool (0.0 < yy
);
3657 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3660 else if (SCM_FRACTIONP (y
))
3662 /* "a/b < c/d" becomes "a*d < c*b" */
3663 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3664 SCM_FRACTION_DENOMINATOR (y
));
3665 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3666 SCM_FRACTION_DENOMINATOR (x
));
3672 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3675 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3679 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3680 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3681 (SCM x
, SCM y
, SCM rest
),
3682 "Return @code{#t} if the list of parameters is monotonically\n"
3684 #define FUNC_NAME s_scm_i_num_gr_p
3686 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3688 while (!scm_is_null (rest
))
3690 if (scm_is_false (scm_gr_p (x
, y
)))
3694 rest
= scm_cdr (rest
);
3696 return scm_gr_p (x
, y
);
3699 #define FUNC_NAME s_scm_i_num_gr_p
3701 scm_gr_p (SCM x
, SCM y
)
3703 if (!SCM_NUMBERP (x
))
3704 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3705 else if (!SCM_NUMBERP (y
))
3706 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3708 return scm_less_p (y
, x
);
3713 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3714 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3715 (SCM x
, SCM y
, SCM rest
),
3716 "Return @code{#t} if the list of parameters is monotonically\n"
3718 #define FUNC_NAME s_scm_i_num_leq_p
3720 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3722 while (!scm_is_null (rest
))
3724 if (scm_is_false (scm_leq_p (x
, y
)))
3728 rest
= scm_cdr (rest
);
3730 return scm_leq_p (x
, y
);
3733 #define FUNC_NAME s_scm_i_num_leq_p
3735 scm_leq_p (SCM x
, SCM y
)
3737 if (!SCM_NUMBERP (x
))
3738 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3739 else if (!SCM_NUMBERP (y
))
3740 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3741 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3744 return scm_not (scm_less_p (y
, x
));
3749 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3750 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3751 (SCM x
, SCM y
, SCM rest
),
3752 "Return @code{#t} if the list of parameters is monotonically\n"
3754 #define FUNC_NAME s_scm_i_num_geq_p
3756 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3758 while (!scm_is_null (rest
))
3760 if (scm_is_false (scm_geq_p (x
, y
)))
3764 rest
= scm_cdr (rest
);
3766 return scm_geq_p (x
, y
);
3769 #define FUNC_NAME s_scm_i_num_geq_p
3771 scm_geq_p (SCM x
, SCM y
)
3773 if (!SCM_NUMBERP (x
))
3774 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3775 else if (!SCM_NUMBERP (y
))
3776 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3777 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3780 return scm_not (scm_less_p (x
, y
));
3785 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3786 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3792 if (SCM_I_INUMP (z
))
3793 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3794 else if (SCM_BIGP (z
))
3796 else if (SCM_REALP (z
))
3797 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3798 else if (SCM_COMPLEXP (z
))
3799 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3800 && SCM_COMPLEX_IMAG (z
) == 0.0);
3801 else if (SCM_FRACTIONP (z
))
3804 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3808 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3809 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3813 scm_positive_p (SCM x
)
3815 if (SCM_I_INUMP (x
))
3816 return scm_from_bool (SCM_I_INUM (x
) > 0);
3817 else if (SCM_BIGP (x
))
3819 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3820 scm_remember_upto_here_1 (x
);
3821 return scm_from_bool (sgn
> 0);
3823 else if (SCM_REALP (x
))
3824 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3825 else if (SCM_FRACTIONP (x
))
3826 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3828 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3832 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3833 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3837 scm_negative_p (SCM x
)
3839 if (SCM_I_INUMP (x
))
3840 return scm_from_bool (SCM_I_INUM (x
) < 0);
3841 else if (SCM_BIGP (x
))
3843 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3844 scm_remember_upto_here_1 (x
);
3845 return scm_from_bool (sgn
< 0);
3847 else if (SCM_REALP (x
))
3848 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3849 else if (SCM_FRACTIONP (x
))
3850 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3852 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3856 /* scm_min and scm_max return an inexact when either argument is inexact, as
3857 required by r5rs. On that basis, for exact/inexact combinations the
3858 exact is converted to inexact to compare and possibly return. This is
3859 unlike scm_less_p above which takes some trouble to preserve all bits in
3860 its test, such trouble is not required for min and max. */
3862 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3863 (SCM x
, SCM y
, SCM rest
),
3864 "Return the maximum of all parameter values.")
3865 #define FUNC_NAME s_scm_i_max
3867 while (!scm_is_null (rest
))
3868 { x
= scm_max (x
, y
);
3870 rest
= scm_cdr (rest
);
3872 return scm_max (x
, y
);
3876 #define s_max s_scm_i_max
3877 #define g_max g_scm_i_max
3880 scm_max (SCM x
, SCM y
)
3885 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3886 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3889 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3892 if (SCM_I_INUMP (x
))
3894 long xx
= SCM_I_INUM (x
);
3895 if (SCM_I_INUMP (y
))
3897 long yy
= SCM_I_INUM (y
);
3898 return (xx
< yy
) ? y
: x
;
3900 else if (SCM_BIGP (y
))
3902 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3903 scm_remember_upto_here_1 (y
);
3904 return (sgn
< 0) ? x
: y
;
3906 else if (SCM_REALP (y
))
3909 /* if y==NaN then ">" is false and we return NaN */
3910 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3912 else if (SCM_FRACTIONP (y
))
3915 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3918 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3920 else if (SCM_BIGP (x
))
3922 if (SCM_I_INUMP (y
))
3924 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3925 scm_remember_upto_here_1 (x
);
3926 return (sgn
< 0) ? y
: x
;
3928 else if (SCM_BIGP (y
))
3930 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3931 scm_remember_upto_here_2 (x
, y
);
3932 return (cmp
> 0) ? x
: y
;
3934 else if (SCM_REALP (y
))
3936 /* if y==NaN then xx>yy is false, so we return the NaN y */
3939 xx
= scm_i_big2dbl (x
);
3940 yy
= SCM_REAL_VALUE (y
);
3941 return (xx
> yy
? scm_from_double (xx
) : y
);
3943 else if (SCM_FRACTIONP (y
))
3948 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3950 else if (SCM_REALP (x
))
3952 if (SCM_I_INUMP (y
))
3954 double z
= SCM_I_INUM (y
);
3955 /* if x==NaN then "<" is false and we return NaN */
3956 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3958 else if (SCM_BIGP (y
))
3963 else if (SCM_REALP (y
))
3965 /* if x==NaN then our explicit check means we return NaN
3966 if y==NaN then ">" is false and we return NaN
3967 calling isnan is unavoidable, since it's the only way to know
3968 which of x or y causes any compares to be false */
3969 double xx
= SCM_REAL_VALUE (x
);
3970 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3972 else if (SCM_FRACTIONP (y
))
3974 double yy
= scm_i_fraction2double (y
);
3975 double xx
= SCM_REAL_VALUE (x
);
3976 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3979 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3981 else if (SCM_FRACTIONP (x
))
3983 if (SCM_I_INUMP (y
))
3987 else if (SCM_BIGP (y
))
3991 else if (SCM_REALP (y
))
3993 double xx
= scm_i_fraction2double (x
);
3994 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3996 else if (SCM_FRACTIONP (y
))
4001 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4004 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4008 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4009 (SCM x
, SCM y
, SCM rest
),
4010 "Return the minimum of all parameter values.")
4011 #define FUNC_NAME s_scm_i_min
4013 while (!scm_is_null (rest
))
4014 { x
= scm_min (x
, y
);
4016 rest
= scm_cdr (rest
);
4018 return scm_min (x
, y
);
4022 #define s_min s_scm_i_min
4023 #define g_min g_scm_i_min
4026 scm_min (SCM x
, SCM y
)
4031 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4032 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4035 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4038 if (SCM_I_INUMP (x
))
4040 long xx
= SCM_I_INUM (x
);
4041 if (SCM_I_INUMP (y
))
4043 long yy
= SCM_I_INUM (y
);
4044 return (xx
< yy
) ? x
: y
;
4046 else if (SCM_BIGP (y
))
4048 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4049 scm_remember_upto_here_1 (y
);
4050 return (sgn
< 0) ? y
: x
;
4052 else if (SCM_REALP (y
))
4055 /* if y==NaN then "<" is false and we return NaN */
4056 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4058 else if (SCM_FRACTIONP (y
))
4061 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4064 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4066 else if (SCM_BIGP (x
))
4068 if (SCM_I_INUMP (y
))
4070 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4071 scm_remember_upto_here_1 (x
);
4072 return (sgn
< 0) ? x
: y
;
4074 else if (SCM_BIGP (y
))
4076 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4077 scm_remember_upto_here_2 (x
, y
);
4078 return (cmp
> 0) ? y
: x
;
4080 else if (SCM_REALP (y
))
4082 /* if y==NaN then xx<yy is false, so we return the NaN y */
4085 xx
= scm_i_big2dbl (x
);
4086 yy
= SCM_REAL_VALUE (y
);
4087 return (xx
< yy
? scm_from_double (xx
) : y
);
4089 else if (SCM_FRACTIONP (y
))
4094 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4096 else if (SCM_REALP (x
))
4098 if (SCM_I_INUMP (y
))
4100 double z
= SCM_I_INUM (y
);
4101 /* if x==NaN then "<" is false and we return NaN */
4102 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4104 else if (SCM_BIGP (y
))
4109 else if (SCM_REALP (y
))
4111 /* if x==NaN then our explicit check means we return NaN
4112 if y==NaN then "<" is false and we return NaN
4113 calling isnan is unavoidable, since it's the only way to know
4114 which of x or y causes any compares to be false */
4115 double xx
= SCM_REAL_VALUE (x
);
4116 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4118 else if (SCM_FRACTIONP (y
))
4120 double yy
= scm_i_fraction2double (y
);
4121 double xx
= SCM_REAL_VALUE (x
);
4122 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4125 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4127 else if (SCM_FRACTIONP (x
))
4129 if (SCM_I_INUMP (y
))
4133 else if (SCM_BIGP (y
))
4137 else if (SCM_REALP (y
))
4139 double xx
= scm_i_fraction2double (x
);
4140 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4142 else if (SCM_FRACTIONP (y
))
4147 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4150 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4154 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4155 (SCM x
, SCM y
, SCM rest
),
4156 "Return the sum of all parameter values. Return 0 if called without\n"
4158 #define FUNC_NAME s_scm_i_sum
4160 while (!scm_is_null (rest
))
4161 { x
= scm_sum (x
, y
);
4163 rest
= scm_cdr (rest
);
4165 return scm_sum (x
, y
);
4169 #define s_sum s_scm_i_sum
4170 #define g_sum g_scm_i_sum
4173 scm_sum (SCM x
, SCM y
)
4175 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4177 if (SCM_NUMBERP (x
)) return x
;
4178 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4179 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4182 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4184 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4186 long xx
= SCM_I_INUM (x
);
4187 long yy
= SCM_I_INUM (y
);
4188 long int z
= xx
+ yy
;
4189 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4191 else if (SCM_BIGP (y
))
4196 else if (SCM_REALP (y
))
4198 long int xx
= SCM_I_INUM (x
);
4199 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4201 else if (SCM_COMPLEXP (y
))
4203 long int xx
= SCM_I_INUM (x
);
4204 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4205 SCM_COMPLEX_IMAG (y
));
4207 else if (SCM_FRACTIONP (y
))
4208 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4209 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4210 SCM_FRACTION_DENOMINATOR (y
));
4212 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4213 } else if (SCM_BIGP (x
))
4215 if (SCM_I_INUMP (y
))
4220 inum
= SCM_I_INUM (y
);
4223 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4226 SCM result
= scm_i_mkbig ();
4227 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4228 scm_remember_upto_here_1 (x
);
4229 /* we know the result will have to be a bignum */
4232 return scm_i_normbig (result
);
4236 SCM result
= scm_i_mkbig ();
4237 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4238 scm_remember_upto_here_1 (x
);
4239 /* we know the result will have to be a bignum */
4242 return scm_i_normbig (result
);
4245 else if (SCM_BIGP (y
))
4247 SCM result
= scm_i_mkbig ();
4248 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4249 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4250 mpz_add (SCM_I_BIG_MPZ (result
),
4253 scm_remember_upto_here_2 (x
, y
);
4254 /* we know the result will have to be a bignum */
4257 return scm_i_normbig (result
);
4259 else if (SCM_REALP (y
))
4261 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4262 scm_remember_upto_here_1 (x
);
4263 return scm_from_double (result
);
4265 else if (SCM_COMPLEXP (y
))
4267 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4268 + SCM_COMPLEX_REAL (y
));
4269 scm_remember_upto_here_1 (x
);
4270 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4272 else if (SCM_FRACTIONP (y
))
4273 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4274 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4275 SCM_FRACTION_DENOMINATOR (y
));
4277 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4279 else if (SCM_REALP (x
))
4281 if (SCM_I_INUMP (y
))
4282 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4283 else if (SCM_BIGP (y
))
4285 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4286 scm_remember_upto_here_1 (y
);
4287 return scm_from_double (result
);
4289 else if (SCM_REALP (y
))
4290 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4291 else if (SCM_COMPLEXP (y
))
4292 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4293 SCM_COMPLEX_IMAG (y
));
4294 else if (SCM_FRACTIONP (y
))
4295 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4297 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4299 else if (SCM_COMPLEXP (x
))
4301 if (SCM_I_INUMP (y
))
4302 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4303 SCM_COMPLEX_IMAG (x
));
4304 else if (SCM_BIGP (y
))
4306 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4307 + SCM_COMPLEX_REAL (x
));
4308 scm_remember_upto_here_1 (y
);
4309 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4311 else if (SCM_REALP (y
))
4312 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4313 SCM_COMPLEX_IMAG (x
));
4314 else if (SCM_COMPLEXP (y
))
4315 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4316 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4317 else if (SCM_FRACTIONP (y
))
4318 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4319 SCM_COMPLEX_IMAG (x
));
4321 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4323 else if (SCM_FRACTIONP (x
))
4325 if (SCM_I_INUMP (y
))
4326 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4327 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4328 SCM_FRACTION_DENOMINATOR (x
));
4329 else if (SCM_BIGP (y
))
4330 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4331 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4332 SCM_FRACTION_DENOMINATOR (x
));
4333 else if (SCM_REALP (y
))
4334 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4335 else if (SCM_COMPLEXP (y
))
4336 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4337 SCM_COMPLEX_IMAG (y
));
4338 else if (SCM_FRACTIONP (y
))
4339 /* a/b + c/d = (ad + bc) / bd */
4340 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4341 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4342 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4344 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4347 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4351 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4353 "Return @math{@var{x}+1}.")
4354 #define FUNC_NAME s_scm_oneplus
4356 return scm_sum (x
, SCM_I_MAKINUM (1));
4361 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4362 (SCM x
, SCM y
, SCM rest
),
4363 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4364 "the sum of all but the first argument are subtracted from the first\n"
4366 #define FUNC_NAME s_scm_i_difference
4368 while (!scm_is_null (rest
))
4369 { x
= scm_difference (x
, y
);
4371 rest
= scm_cdr (rest
);
4373 return scm_difference (x
, y
);
4377 #define s_difference s_scm_i_difference
4378 #define g_difference g_scm_i_difference
4381 scm_difference (SCM x
, SCM y
)
4382 #define FUNC_NAME s_difference
4384 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4387 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4389 if (SCM_I_INUMP (x
))
4391 long xx
= -SCM_I_INUM (x
);
4392 if (SCM_FIXABLE (xx
))
4393 return SCM_I_MAKINUM (xx
);
4395 return scm_i_long2big (xx
);
4397 else if (SCM_BIGP (x
))
4398 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4399 bignum, but negating that gives a fixnum. */
4400 return scm_i_normbig (scm_i_clonebig (x
, 0));
4401 else if (SCM_REALP (x
))
4402 return scm_from_double (-SCM_REAL_VALUE (x
));
4403 else if (SCM_COMPLEXP (x
))
4404 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4405 -SCM_COMPLEX_IMAG (x
));
4406 else if (SCM_FRACTIONP (x
))
4407 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4408 SCM_FRACTION_DENOMINATOR (x
));
4410 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4413 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4415 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4417 long int xx
= SCM_I_INUM (x
);
4418 long int yy
= SCM_I_INUM (y
);
4419 long int z
= xx
- yy
;
4420 if (SCM_FIXABLE (z
))
4421 return SCM_I_MAKINUM (z
);
4423 return scm_i_long2big (z
);
4425 else if (SCM_BIGP (y
))
4427 /* inum-x - big-y */
4428 long xx
= SCM_I_INUM (x
);
4431 return scm_i_clonebig (y
, 0);
4434 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4435 SCM result
= scm_i_mkbig ();
4438 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4441 /* x - y == -(y + -x) */
4442 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4443 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4445 scm_remember_upto_here_1 (y
);
4447 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4448 /* we know the result will have to be a bignum */
4451 return scm_i_normbig (result
);
4454 else if (SCM_REALP (y
))
4456 long int xx
= SCM_I_INUM (x
);
4457 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4459 else if (SCM_COMPLEXP (y
))
4461 long int xx
= SCM_I_INUM (x
);
4462 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4463 - SCM_COMPLEX_IMAG (y
));
4465 else if (SCM_FRACTIONP (y
))
4466 /* a - b/c = (ac - b) / c */
4467 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4468 SCM_FRACTION_NUMERATOR (y
)),
4469 SCM_FRACTION_DENOMINATOR (y
));
4471 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4473 else if (SCM_BIGP (x
))
4475 if (SCM_I_INUMP (y
))
4477 /* big-x - inum-y */
4478 long yy
= SCM_I_INUM (y
);
4479 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4481 scm_remember_upto_here_1 (x
);
4483 return (SCM_FIXABLE (-yy
) ?
4484 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4487 SCM result
= scm_i_mkbig ();
4490 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4492 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4493 scm_remember_upto_here_1 (x
);
4495 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4496 /* we know the result will have to be a bignum */
4499 return scm_i_normbig (result
);
4502 else if (SCM_BIGP (y
))
4504 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4505 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4506 SCM result
= scm_i_mkbig ();
4507 mpz_sub (SCM_I_BIG_MPZ (result
),
4510 scm_remember_upto_here_2 (x
, y
);
4511 /* we know the result will have to be a bignum */
4512 if ((sgn_x
== 1) && (sgn_y
== -1))
4514 if ((sgn_x
== -1) && (sgn_y
== 1))
4516 return scm_i_normbig (result
);
4518 else if (SCM_REALP (y
))
4520 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4521 scm_remember_upto_here_1 (x
);
4522 return scm_from_double (result
);
4524 else if (SCM_COMPLEXP (y
))
4526 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4527 - SCM_COMPLEX_REAL (y
));
4528 scm_remember_upto_here_1 (x
);
4529 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4531 else if (SCM_FRACTIONP (y
))
4532 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4533 SCM_FRACTION_NUMERATOR (y
)),
4534 SCM_FRACTION_DENOMINATOR (y
));
4535 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4537 else if (SCM_REALP (x
))
4539 if (SCM_I_INUMP (y
))
4540 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4541 else if (SCM_BIGP (y
))
4543 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4544 scm_remember_upto_here_1 (x
);
4545 return scm_from_double (result
);
4547 else if (SCM_REALP (y
))
4548 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4549 else if (SCM_COMPLEXP (y
))
4550 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4551 -SCM_COMPLEX_IMAG (y
));
4552 else if (SCM_FRACTIONP (y
))
4553 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4555 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4557 else if (SCM_COMPLEXP (x
))
4559 if (SCM_I_INUMP (y
))
4560 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4561 SCM_COMPLEX_IMAG (x
));
4562 else if (SCM_BIGP (y
))
4564 double real_part
= (SCM_COMPLEX_REAL (x
)
4565 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4566 scm_remember_upto_here_1 (x
);
4567 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4569 else if (SCM_REALP (y
))
4570 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4571 SCM_COMPLEX_IMAG (x
));
4572 else if (SCM_COMPLEXP (y
))
4573 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4574 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4575 else if (SCM_FRACTIONP (y
))
4576 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4577 SCM_COMPLEX_IMAG (x
));
4579 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4581 else if (SCM_FRACTIONP (x
))
4583 if (SCM_I_INUMP (y
))
4584 /* a/b - c = (a - cb) / b */
4585 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4586 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4587 SCM_FRACTION_DENOMINATOR (x
));
4588 else if (SCM_BIGP (y
))
4589 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4590 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4591 SCM_FRACTION_DENOMINATOR (x
));
4592 else if (SCM_REALP (y
))
4593 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4594 else if (SCM_COMPLEXP (y
))
4595 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4596 -SCM_COMPLEX_IMAG (y
));
4597 else if (SCM_FRACTIONP (y
))
4598 /* a/b - c/d = (ad - bc) / bd */
4599 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4600 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4601 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4603 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4606 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4611 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4613 "Return @math{@var{x}-1}.")
4614 #define FUNC_NAME s_scm_oneminus
4616 return scm_difference (x
, SCM_I_MAKINUM (1));
4621 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4622 (SCM x
, SCM y
, SCM rest
),
4623 "Return the product of all arguments. If called without arguments,\n"
4625 #define FUNC_NAME s_scm_i_product
4627 while (!scm_is_null (rest
))
4628 { x
= scm_product (x
, y
);
4630 rest
= scm_cdr (rest
);
4632 return scm_product (x
, y
);
4636 #define s_product s_scm_i_product
4637 #define g_product g_scm_i_product
4640 scm_product (SCM x
, SCM y
)
4642 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4645 return SCM_I_MAKINUM (1L);
4646 else if (SCM_NUMBERP (x
))
4649 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4652 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4657 xx
= SCM_I_INUM (x
);
4661 case 0: return x
; break;
4662 case 1: return y
; break;
4665 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4667 long yy
= SCM_I_INUM (y
);
4669 SCM k
= SCM_I_MAKINUM (kk
);
4670 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4674 SCM result
= scm_i_long2big (xx
);
4675 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4676 return scm_i_normbig (result
);
4679 else if (SCM_BIGP (y
))
4681 SCM result
= scm_i_mkbig ();
4682 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4683 scm_remember_upto_here_1 (y
);
4686 else if (SCM_REALP (y
))
4687 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4688 else if (SCM_COMPLEXP (y
))
4689 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4690 xx
* SCM_COMPLEX_IMAG (y
));
4691 else if (SCM_FRACTIONP (y
))
4692 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4693 SCM_FRACTION_DENOMINATOR (y
));
4695 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4697 else if (SCM_BIGP (x
))
4699 if (SCM_I_INUMP (y
))
4704 else if (SCM_BIGP (y
))
4706 SCM result
= scm_i_mkbig ();
4707 mpz_mul (SCM_I_BIG_MPZ (result
),
4710 scm_remember_upto_here_2 (x
, y
);
4713 else if (SCM_REALP (y
))
4715 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4716 scm_remember_upto_here_1 (x
);
4717 return scm_from_double (result
);
4719 else if (SCM_COMPLEXP (y
))
4721 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4722 scm_remember_upto_here_1 (x
);
4723 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4724 z
* SCM_COMPLEX_IMAG (y
));
4726 else if (SCM_FRACTIONP (y
))
4727 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4728 SCM_FRACTION_DENOMINATOR (y
));
4730 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4732 else if (SCM_REALP (x
))
4734 if (SCM_I_INUMP (y
))
4736 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4737 if (scm_is_eq (y
, SCM_INUM0
))
4739 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4741 else if (SCM_BIGP (y
))
4743 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4744 scm_remember_upto_here_1 (y
);
4745 return scm_from_double (result
);
4747 else if (SCM_REALP (y
))
4748 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4749 else if (SCM_COMPLEXP (y
))
4750 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4751 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4752 else if (SCM_FRACTIONP (y
))
4753 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4755 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4757 else if (SCM_COMPLEXP (x
))
4759 if (SCM_I_INUMP (y
))
4761 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4762 if (scm_is_eq (y
, SCM_INUM0
))
4764 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4765 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4767 else if (SCM_BIGP (y
))
4769 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4770 scm_remember_upto_here_1 (y
);
4771 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4772 z
* SCM_COMPLEX_IMAG (x
));
4774 else if (SCM_REALP (y
))
4775 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4776 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4777 else if (SCM_COMPLEXP (y
))
4779 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4780 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4781 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4782 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4784 else if (SCM_FRACTIONP (y
))
4786 double yy
= scm_i_fraction2double (y
);
4787 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4788 yy
* SCM_COMPLEX_IMAG (x
));
4791 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4793 else if (SCM_FRACTIONP (x
))
4795 if (SCM_I_INUMP (y
))
4796 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4797 SCM_FRACTION_DENOMINATOR (x
));
4798 else if (SCM_BIGP (y
))
4799 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4800 SCM_FRACTION_DENOMINATOR (x
));
4801 else if (SCM_REALP (y
))
4802 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4803 else if (SCM_COMPLEXP (y
))
4805 double xx
= scm_i_fraction2double (x
);
4806 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4807 xx
* SCM_COMPLEX_IMAG (y
));
4809 else if (SCM_FRACTIONP (y
))
4810 /* a/b * c/d = ac / bd */
4811 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4812 SCM_FRACTION_NUMERATOR (y
)),
4813 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4814 SCM_FRACTION_DENOMINATOR (y
)));
4816 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4819 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4822 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4823 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4824 #define ALLOW_DIVIDE_BY_ZERO
4825 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4828 /* The code below for complex division is adapted from the GNU
4829 libstdc++, which adapted it from f2c's libF77, and is subject to
4832 /****************************************************************
4833 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4835 Permission to use, copy, modify, and distribute this software
4836 and its documentation for any purpose and without fee is hereby
4837 granted, provided that the above copyright notice appear in all
4838 copies and that both that the copyright notice and this
4839 permission notice and warranty disclaimer appear in supporting
4840 documentation, and that the names of AT&T Bell Laboratories or
4841 Bellcore or any of their entities not be used in advertising or
4842 publicity pertaining to distribution of the software without
4843 specific, written prior permission.
4845 AT&T and Bellcore disclaim all warranties with regard to this
4846 software, including all implied warranties of merchantability
4847 and fitness. In no event shall AT&T or Bellcore be liable for
4848 any special, indirect or consequential damages or any damages
4849 whatsoever resulting from loss of use, data or profits, whether
4850 in an action of contract, negligence or other tortious action,
4851 arising out of or in connection with the use or performance of
4853 ****************************************************************/
4855 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4856 (SCM x
, SCM y
, SCM rest
),
4857 "Divide the first argument by the product of the remaining\n"
4858 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4860 #define FUNC_NAME s_scm_i_divide
4862 while (!scm_is_null (rest
))
4863 { x
= scm_divide (x
, y
);
4865 rest
= scm_cdr (rest
);
4867 return scm_divide (x
, y
);
4871 #define s_divide s_scm_i_divide
4872 #define g_divide g_scm_i_divide
4875 do_divide (SCM x
, SCM y
, int inexact
)
4876 #define FUNC_NAME s_divide
4880 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4883 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4884 else if (SCM_I_INUMP (x
))
4886 long xx
= SCM_I_INUM (x
);
4887 if (xx
== 1 || xx
== -1)
4889 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4891 scm_num_overflow (s_divide
);
4896 return scm_from_double (1.0 / (double) xx
);
4897 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4900 else if (SCM_BIGP (x
))
4903 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4904 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4906 else if (SCM_REALP (x
))
4908 double xx
= SCM_REAL_VALUE (x
);
4909 #ifndef ALLOW_DIVIDE_BY_ZERO
4911 scm_num_overflow (s_divide
);
4914 return scm_from_double (1.0 / xx
);
4916 else if (SCM_COMPLEXP (x
))
4918 double r
= SCM_COMPLEX_REAL (x
);
4919 double i
= SCM_COMPLEX_IMAG (x
);
4920 if (fabs(r
) <= fabs(i
))
4923 double d
= i
* (1.0 + t
* t
);
4924 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4929 double d
= r
* (1.0 + t
* t
);
4930 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4933 else if (SCM_FRACTIONP (x
))
4934 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4935 SCM_FRACTION_NUMERATOR (x
));
4937 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4940 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4942 long xx
= SCM_I_INUM (x
);
4943 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4945 long yy
= SCM_I_INUM (y
);
4948 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4949 scm_num_overflow (s_divide
);
4951 return scm_from_double ((double) xx
/ (double) yy
);
4954 else if (xx
% yy
!= 0)
4957 return scm_from_double ((double) xx
/ (double) yy
);
4958 else return scm_i_make_ratio (x
, y
);
4963 if (SCM_FIXABLE (z
))
4964 return SCM_I_MAKINUM (z
);
4966 return scm_i_long2big (z
);
4969 else if (SCM_BIGP (y
))
4972 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4973 else return scm_i_make_ratio (x
, y
);
4975 else if (SCM_REALP (y
))
4977 double yy
= SCM_REAL_VALUE (y
);
4978 #ifndef ALLOW_DIVIDE_BY_ZERO
4980 scm_num_overflow (s_divide
);
4983 return scm_from_double ((double) xx
/ yy
);
4985 else if (SCM_COMPLEXP (y
))
4988 complex_div
: /* y _must_ be a complex number */
4990 double r
= SCM_COMPLEX_REAL (y
);
4991 double i
= SCM_COMPLEX_IMAG (y
);
4992 if (fabs(r
) <= fabs(i
))
4995 double d
= i
* (1.0 + t
* t
);
4996 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5001 double d
= r
* (1.0 + t
* t
);
5002 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5006 else if (SCM_FRACTIONP (y
))
5007 /* a / b/c = ac / b */
5008 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5009 SCM_FRACTION_NUMERATOR (y
));
5011 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5013 else if (SCM_BIGP (x
))
5015 if (SCM_I_INUMP (y
))
5017 long int yy
= SCM_I_INUM (y
);
5020 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5021 scm_num_overflow (s_divide
);
5023 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5024 scm_remember_upto_here_1 (x
);
5025 return (sgn
== 0) ? scm_nan () : scm_inf ();
5032 /* FIXME: HMM, what are the relative performance issues here?
5033 We need to test. Is it faster on average to test
5034 divisible_p, then perform whichever operation, or is it
5035 faster to perform the integer div opportunistically and
5036 switch to real if there's a remainder? For now we take the
5037 middle ground: test, then if divisible, use the faster div
5040 long abs_yy
= yy
< 0 ? -yy
: yy
;
5041 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5045 SCM result
= scm_i_mkbig ();
5046 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5047 scm_remember_upto_here_1 (x
);
5049 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5050 return scm_i_normbig (result
);
5055 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5056 else return scm_i_make_ratio (x
, y
);
5060 else if (SCM_BIGP (y
))
5062 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5065 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5066 scm_num_overflow (s_divide
);
5068 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5069 scm_remember_upto_here_1 (x
);
5070 return (sgn
== 0) ? scm_nan () : scm_inf ();
5078 /* It's easily possible for the ratio x/y to fit a double
5079 but one or both x and y be too big to fit a double,
5080 hence the use of mpq_get_d rather than converting and
5083 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5084 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5085 return scm_from_double (mpq_get_d (q
));
5089 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5093 SCM result
= scm_i_mkbig ();
5094 mpz_divexact (SCM_I_BIG_MPZ (result
),
5097 scm_remember_upto_here_2 (x
, y
);
5098 return scm_i_normbig (result
);
5101 return scm_i_make_ratio (x
, y
);
5105 else if (SCM_REALP (y
))
5107 double yy
= SCM_REAL_VALUE (y
);
5108 #ifndef ALLOW_DIVIDE_BY_ZERO
5110 scm_num_overflow (s_divide
);
5113 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5115 else if (SCM_COMPLEXP (y
))
5117 a
= scm_i_big2dbl (x
);
5120 else if (SCM_FRACTIONP (y
))
5121 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5122 SCM_FRACTION_NUMERATOR (y
));
5124 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5126 else if (SCM_REALP (x
))
5128 double rx
= SCM_REAL_VALUE (x
);
5129 if (SCM_I_INUMP (y
))
5131 long int yy
= SCM_I_INUM (y
);
5132 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5134 scm_num_overflow (s_divide
);
5137 return scm_from_double (rx
/ (double) yy
);
5139 else if (SCM_BIGP (y
))
5141 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5142 scm_remember_upto_here_1 (y
);
5143 return scm_from_double (rx
/ dby
);
5145 else if (SCM_REALP (y
))
5147 double yy
= SCM_REAL_VALUE (y
);
5148 #ifndef ALLOW_DIVIDE_BY_ZERO
5150 scm_num_overflow (s_divide
);
5153 return scm_from_double (rx
/ yy
);
5155 else if (SCM_COMPLEXP (y
))
5160 else if (SCM_FRACTIONP (y
))
5161 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5163 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5165 else if (SCM_COMPLEXP (x
))
5167 double rx
= SCM_COMPLEX_REAL (x
);
5168 double ix
= SCM_COMPLEX_IMAG (x
);
5169 if (SCM_I_INUMP (y
))
5171 long int yy
= SCM_I_INUM (y
);
5172 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5174 scm_num_overflow (s_divide
);
5179 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5182 else if (SCM_BIGP (y
))
5184 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5185 scm_remember_upto_here_1 (y
);
5186 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5188 else if (SCM_REALP (y
))
5190 double yy
= SCM_REAL_VALUE (y
);
5191 #ifndef ALLOW_DIVIDE_BY_ZERO
5193 scm_num_overflow (s_divide
);
5196 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5198 else if (SCM_COMPLEXP (y
))
5200 double ry
= SCM_COMPLEX_REAL (y
);
5201 double iy
= SCM_COMPLEX_IMAG (y
);
5202 if (fabs(ry
) <= fabs(iy
))
5205 double d
= iy
* (1.0 + t
* t
);
5206 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5211 double d
= ry
* (1.0 + t
* t
);
5212 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5215 else if (SCM_FRACTIONP (y
))
5217 double yy
= scm_i_fraction2double (y
);
5218 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5221 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5223 else if (SCM_FRACTIONP (x
))
5225 if (SCM_I_INUMP (y
))
5227 long int yy
= SCM_I_INUM (y
);
5228 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5230 scm_num_overflow (s_divide
);
5233 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5234 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5236 else if (SCM_BIGP (y
))
5238 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5239 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5241 else if (SCM_REALP (y
))
5243 double yy
= SCM_REAL_VALUE (y
);
5244 #ifndef ALLOW_DIVIDE_BY_ZERO
5246 scm_num_overflow (s_divide
);
5249 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5251 else if (SCM_COMPLEXP (y
))
5253 a
= scm_i_fraction2double (x
);
5256 else if (SCM_FRACTIONP (y
))
5257 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5258 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5260 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5263 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5267 scm_divide (SCM x
, SCM y
)
5269 return do_divide (x
, y
, 0);
5272 static SCM
scm_divide2real (SCM x
, SCM y
)
5274 return do_divide (x
, y
, 1);
5280 scm_c_truncate (double x
)
5291 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5292 half-way case (ie. when x is an integer plus 0.5) going upwards.
5293 Then half-way cases are identified and adjusted down if the
5294 round-upwards didn't give the desired even integer.
5296 "plus_half == result" identifies a half-way case. If plus_half, which is
5297 x + 0.5, is an integer then x must be an integer plus 0.5.
5299 An odd "result" value is identified with result/2 != floor(result/2).
5300 This is done with plus_half, since that value is ready for use sooner in
5301 a pipelined cpu, and we're already requiring plus_half == result.
5303 Note however that we need to be careful when x is big and already an
5304 integer. In that case "x+0.5" may round to an adjacent integer, causing
5305 us to return such a value, incorrectly. For instance if the hardware is
5306 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5307 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5308 returned. Or if the hardware is in round-upwards mode, then other bigger
5309 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5310 representable value, 2^128+2^76 (or whatever), again incorrect.
5312 These bad roundings of x+0.5 are avoided by testing at the start whether
5313 x is already an integer. If it is then clearly that's the desired result
5314 already. And if it's not then the exponent must be small enough to allow
5315 an 0.5 to be represented, and hence added without a bad rounding. */
5318 scm_c_round (double x
)
5320 double plus_half
, result
;
5325 plus_half
= x
+ 0.5;
5326 result
= floor (plus_half
);
5327 /* Adjust so that the rounding is towards even. */
5328 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5333 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5335 "Round the number @var{x} towards zero.")
5336 #define FUNC_NAME s_scm_truncate_number
5338 if (scm_is_false (scm_negative_p (x
)))
5339 return scm_floor (x
);
5341 return scm_ceiling (x
);
5345 static SCM exactly_one_half
;
5347 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5349 "Round the number @var{x} towards the nearest integer. "
5350 "When it is exactly halfway between two integers, "
5351 "round towards the even one.")
5352 #define FUNC_NAME s_scm_round_number
5354 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5356 else if (SCM_REALP (x
))
5357 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5360 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5361 single quotient+remainder division then examining to see which way
5362 the rounding should go. */
5363 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5364 SCM result
= scm_floor (plus_half
);
5365 /* Adjust so that the rounding is towards even. */
5366 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5367 && scm_is_true (scm_odd_p (result
)))
5368 return scm_difference (result
, SCM_I_MAKINUM (1));
5375 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5377 "Round the number @var{x} towards minus infinity.")
5378 #define FUNC_NAME s_scm_floor
5380 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5382 else if (SCM_REALP (x
))
5383 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5384 else if (SCM_FRACTIONP (x
))
5386 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5387 SCM_FRACTION_DENOMINATOR (x
));
5388 if (scm_is_false (scm_negative_p (x
)))
5390 /* For positive x, rounding towards zero is correct. */
5395 /* For negative x, we need to return q-1 unless x is an
5396 integer. But fractions are never integer, per our
5398 return scm_difference (q
, SCM_I_MAKINUM (1));
5402 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5406 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5408 "Round the number @var{x} towards infinity.")
5409 #define FUNC_NAME s_scm_ceiling
5411 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5413 else if (SCM_REALP (x
))
5414 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5415 else if (SCM_FRACTIONP (x
))
5417 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5418 SCM_FRACTION_DENOMINATOR (x
));
5419 if (scm_is_false (scm_positive_p (x
)))
5421 /* For negative x, rounding towards zero is correct. */
5426 /* For positive x, we need to return q+1 unless x is an
5427 integer. But fractions are never integer, per our
5429 return scm_sum (q
, SCM_I_MAKINUM (1));
5433 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5437 /* sin/cos/tan/asin/acos/atan
5438 sinh/cosh/tanh/asinh/acosh/atanh
5439 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5440 Written by Jerry D. Hedden, (C) FSF.
5441 See the file `COPYING' for terms applying to this program. */
5443 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5445 "Return @var{x} raised to the power of @var{y}.")
5446 #define FUNC_NAME s_scm_expt
5448 if (scm_is_true (scm_exact_p (x
)) && scm_is_integer (y
))
5449 return scm_integer_expt (x
, y
);
5450 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5452 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5455 return scm_exp (scm_product (scm_log (x
), y
));
5459 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5461 "Compute the sine of @var{z}.")
5462 #define FUNC_NAME s_scm_sin
5464 if (scm_is_real (z
))
5465 return scm_from_double (sin (scm_to_double (z
)));
5466 else if (SCM_COMPLEXP (z
))
5468 x
= SCM_COMPLEX_REAL (z
);
5469 y
= SCM_COMPLEX_IMAG (z
);
5470 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5471 cos (x
) * sinh (y
));
5474 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5478 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5480 "Compute the cosine of @var{z}.")
5481 #define FUNC_NAME s_scm_cos
5483 if (scm_is_real (z
))
5484 return scm_from_double (cos (scm_to_double (z
)));
5485 else if (SCM_COMPLEXP (z
))
5487 x
= SCM_COMPLEX_REAL (z
);
5488 y
= SCM_COMPLEX_IMAG (z
);
5489 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5490 -sin (x
) * sinh (y
));
5493 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5497 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5499 "Compute the tangent of @var{z}.")
5500 #define FUNC_NAME s_scm_tan
5502 if (scm_is_real (z
))
5503 return scm_from_double (tan (scm_to_double (z
)));
5504 else if (SCM_COMPLEXP (z
))
5506 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5507 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5508 w
= cos (x
) + cosh (y
);
5509 #ifndef ALLOW_DIVIDE_BY_ZERO
5511 scm_num_overflow (s_scm_tan
);
5513 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5516 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5520 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5522 "Compute the hyperbolic sine of @var{z}.")
5523 #define FUNC_NAME s_scm_sinh
5525 if (scm_is_real (z
))
5526 return scm_from_double (sinh (scm_to_double (z
)));
5527 else if (SCM_COMPLEXP (z
))
5529 x
= SCM_COMPLEX_REAL (z
);
5530 y
= SCM_COMPLEX_IMAG (z
);
5531 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5532 cosh (x
) * sin (y
));
5535 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5539 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5541 "Compute the hyperbolic cosine of @var{z}.")
5542 #define FUNC_NAME s_scm_cosh
5544 if (scm_is_real (z
))
5545 return scm_from_double (cosh (scm_to_double (z
)));
5546 else if (SCM_COMPLEXP (z
))
5548 x
= SCM_COMPLEX_REAL (z
);
5549 y
= SCM_COMPLEX_IMAG (z
);
5550 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5551 sinh (x
) * sin (y
));
5554 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5558 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5560 "Compute the hyperbolic tangent of @var{z}.")
5561 #define FUNC_NAME s_scm_tanh
5563 if (scm_is_real (z
))
5564 return scm_from_double (tanh (scm_to_double (z
)));
5565 else if (SCM_COMPLEXP (z
))
5567 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5568 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5569 w
= cosh (x
) + cos (y
);
5570 #ifndef ALLOW_DIVIDE_BY_ZERO
5572 scm_num_overflow (s_scm_tanh
);
5574 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5577 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5581 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5583 "Compute the arc sine of @var{z}.")
5584 #define FUNC_NAME s_scm_asin
5586 if (scm_is_real (z
))
5588 double w
= scm_to_double (z
);
5589 if (w
>= -1.0 && w
<= 1.0)
5590 return scm_from_double (asin (w
));
5592 return scm_product (scm_c_make_rectangular (0, -1),
5593 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5595 else if (SCM_COMPLEXP (z
))
5597 x
= SCM_COMPLEX_REAL (z
);
5598 y
= SCM_COMPLEX_IMAG (z
);
5599 return scm_product (scm_c_make_rectangular (0, -1),
5600 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5603 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5607 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5609 "Compute the arc cosine of @var{z}.")
5610 #define FUNC_NAME s_scm_acos
5612 if (scm_is_real (z
))
5614 double w
= scm_to_double (z
);
5615 if (w
>= -1.0 && w
<= 1.0)
5616 return scm_from_double (acos (w
));
5618 return scm_sum (scm_from_double (acos (0.0)),
5619 scm_product (scm_c_make_rectangular (0, 1),
5620 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5622 else if (SCM_COMPLEXP (z
))
5624 x
= SCM_COMPLEX_REAL (z
);
5625 y
= SCM_COMPLEX_IMAG (z
);
5626 return scm_sum (scm_from_double (acos (0.0)),
5627 scm_product (scm_c_make_rectangular (0, 1),
5628 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5631 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5635 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5637 "With one argument, compute the arc tangent of @var{z}.\n"
5638 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5639 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5640 #define FUNC_NAME s_scm_atan
5644 if (scm_is_real (z
))
5645 return scm_from_double (atan (scm_to_double (z
)));
5646 else if (SCM_COMPLEXP (z
))
5649 v
= SCM_COMPLEX_REAL (z
);
5650 w
= SCM_COMPLEX_IMAG (z
);
5651 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5652 scm_c_make_rectangular (v
, w
+ 1.0))),
5653 scm_c_make_rectangular (0, 2));
5656 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5658 else if (scm_is_real (z
))
5660 if (scm_is_real (y
))
5661 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5663 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5666 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5670 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5672 "Compute the inverse hyperbolic sine of @var{z}.")
5673 #define FUNC_NAME s_scm_sys_asinh
5675 if (scm_is_real (z
))
5676 return scm_from_double (asinh (scm_to_double (z
)));
5677 else if (scm_is_number (z
))
5678 return scm_log (scm_sum (z
,
5679 scm_sqrt (scm_sum (scm_product (z
, z
),
5680 SCM_I_MAKINUM (1)))));
5682 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5686 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5688 "Compute the inverse hyperbolic cosine of @var{z}.")
5689 #define FUNC_NAME s_scm_sys_acosh
5691 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5692 return scm_from_double (acosh (scm_to_double (z
)));
5693 else if (scm_is_number (z
))
5694 return scm_log (scm_sum (z
,
5695 scm_sqrt (scm_difference (scm_product (z
, z
),
5696 SCM_I_MAKINUM (1)))));
5698 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5702 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5704 "Compute the inverse hyperbolic tangent of @var{z}.")
5705 #define FUNC_NAME s_scm_sys_atanh
5707 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5708 return scm_from_double (atanh (scm_to_double (z
)));
5709 else if (scm_is_number (z
))
5710 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5711 scm_difference (SCM_I_MAKINUM (1), z
))),
5714 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5719 scm_c_make_rectangular (double re
, double im
)
5722 return scm_from_double (re
);
5726 SCM_NEWSMOB (z
, scm_tc16_complex
,
5727 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5729 SCM_COMPLEX_REAL (z
) = re
;
5730 SCM_COMPLEX_IMAG (z
) = im
;
5735 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5736 (SCM real_part
, SCM imaginary_part
),
5737 "Return a complex number constructed of the given @var{real-part} "
5738 "and @var{imaginary-part} parts.")
5739 #define FUNC_NAME s_scm_make_rectangular
5741 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5742 SCM_ARG1
, FUNC_NAME
, "real");
5743 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5744 SCM_ARG2
, FUNC_NAME
, "real");
5745 return scm_c_make_rectangular (scm_to_double (real_part
),
5746 scm_to_double (imaginary_part
));
5751 scm_c_make_polar (double mag
, double ang
)
5755 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5756 use it on Glibc-based systems that have it (it's a GNU extension). See
5757 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5759 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5760 sincos (ang
, &s
, &c
);
5765 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5768 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5770 "Return the complex number @var{x} * e^(i * @var{y}).")
5771 #define FUNC_NAME s_scm_make_polar
5773 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5774 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5775 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5780 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5781 /* "Return the real part of the number @var{z}."
5784 scm_real_part (SCM z
)
5786 if (SCM_I_INUMP (z
))
5788 else if (SCM_BIGP (z
))
5790 else if (SCM_REALP (z
))
5792 else if (SCM_COMPLEXP (z
))
5793 return scm_from_double (SCM_COMPLEX_REAL (z
));
5794 else if (SCM_FRACTIONP (z
))
5797 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5801 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5802 /* "Return the imaginary part of the number @var{z}."
5805 scm_imag_part (SCM z
)
5807 if (SCM_I_INUMP (z
))
5809 else if (SCM_BIGP (z
))
5811 else if (SCM_REALP (z
))
5813 else if (SCM_COMPLEXP (z
))
5814 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5815 else if (SCM_FRACTIONP (z
))
5818 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5821 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5822 /* "Return the numerator of the number @var{z}."
5825 scm_numerator (SCM z
)
5827 if (SCM_I_INUMP (z
))
5829 else if (SCM_BIGP (z
))
5831 else if (SCM_FRACTIONP (z
))
5832 return SCM_FRACTION_NUMERATOR (z
);
5833 else if (SCM_REALP (z
))
5834 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5836 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5840 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5841 /* "Return the denominator of the number @var{z}."
5844 scm_denominator (SCM z
)
5846 if (SCM_I_INUMP (z
))
5847 return SCM_I_MAKINUM (1);
5848 else if (SCM_BIGP (z
))
5849 return SCM_I_MAKINUM (1);
5850 else if (SCM_FRACTIONP (z
))
5851 return SCM_FRACTION_DENOMINATOR (z
);
5852 else if (SCM_REALP (z
))
5853 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5855 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5858 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5859 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5860 * "@code{abs} for real arguments, but also allows complex numbers."
5863 scm_magnitude (SCM z
)
5865 if (SCM_I_INUMP (z
))
5867 long int zz
= SCM_I_INUM (z
);
5870 else if (SCM_POSFIXABLE (-zz
))
5871 return SCM_I_MAKINUM (-zz
);
5873 return scm_i_long2big (-zz
);
5875 else if (SCM_BIGP (z
))
5877 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5878 scm_remember_upto_here_1 (z
);
5880 return scm_i_clonebig (z
, 0);
5884 else if (SCM_REALP (z
))
5885 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5886 else if (SCM_COMPLEXP (z
))
5887 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5888 else if (SCM_FRACTIONP (z
))
5890 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5892 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5893 SCM_FRACTION_DENOMINATOR (z
));
5896 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5900 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5901 /* "Return the angle of the complex number @var{z}."
5906 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5907 flo0 to save allocating a new flonum with scm_from_double each time.
5908 But if atan2 follows the floating point rounding mode, then the value
5909 is not a constant. Maybe it'd be close enough though. */
5910 if (SCM_I_INUMP (z
))
5912 if (SCM_I_INUM (z
) >= 0)
5915 return scm_from_double (atan2 (0.0, -1.0));
5917 else if (SCM_BIGP (z
))
5919 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5920 scm_remember_upto_here_1 (z
);
5922 return scm_from_double (atan2 (0.0, -1.0));
5926 else if (SCM_REALP (z
))
5928 if (SCM_REAL_VALUE (z
) >= 0)
5931 return scm_from_double (atan2 (0.0, -1.0));
5933 else if (SCM_COMPLEXP (z
))
5934 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5935 else if (SCM_FRACTIONP (z
))
5937 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5939 else return scm_from_double (atan2 (0.0, -1.0));
5942 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5946 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5947 /* Convert the number @var{x} to its inexact representation.\n"
5950 scm_exact_to_inexact (SCM z
)
5952 if (SCM_I_INUMP (z
))
5953 return scm_from_double ((double) SCM_I_INUM (z
));
5954 else if (SCM_BIGP (z
))
5955 return scm_from_double (scm_i_big2dbl (z
));
5956 else if (SCM_FRACTIONP (z
))
5957 return scm_from_double (scm_i_fraction2double (z
));
5958 else if (SCM_INEXACTP (z
))
5961 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5965 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5967 "Return an exact number that is numerically closest to @var{z}.")
5968 #define FUNC_NAME s_scm_inexact_to_exact
5970 if (SCM_I_INUMP (z
))
5972 else if (SCM_BIGP (z
))
5974 else if (SCM_REALP (z
))
5976 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
5977 SCM_OUT_OF_RANGE (1, z
);
5984 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5985 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5986 scm_i_mpz2num (mpq_denref (frac
)));
5988 /* When scm_i_make_ratio throws, we leak the memory allocated
5995 else if (SCM_FRACTIONP (z
))
5998 SCM_WRONG_TYPE_ARG (1, z
);
6002 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6004 "Returns the @emph{simplest} rational number differing\n"
6005 "from @var{x} by no more than @var{eps}.\n"
6007 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6008 "exact result when both its arguments are exact. Thus, you might need\n"
6009 "to use @code{inexact->exact} on the arguments.\n"
6012 "(rationalize (inexact->exact 1.2) 1/100)\n"
6015 #define FUNC_NAME s_scm_rationalize
6017 if (SCM_I_INUMP (x
))
6019 else if (SCM_BIGP (x
))
6021 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6023 /* Use continued fractions to find closest ratio. All
6024 arithmetic is done with exact numbers.
6027 SCM ex
= scm_inexact_to_exact (x
);
6028 SCM int_part
= scm_floor (ex
);
6029 SCM tt
= SCM_I_MAKINUM (1);
6030 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
6031 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
6035 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6038 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6039 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6041 /* We stop after a million iterations just to be absolutely sure
6042 that we don't go into an infinite loop. The process normally
6043 converges after less than a dozen iterations.
6046 eps
= scm_abs (eps
);
6047 while (++i
< 1000000)
6049 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6050 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6051 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6053 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6054 eps
))) /* abs(x-a/b) <= eps */
6056 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6057 if (scm_is_false (scm_exact_p (x
))
6058 || scm_is_false (scm_exact_p (eps
)))
6059 return scm_exact_to_inexact (res
);
6063 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6065 tt
= scm_floor (rx
); /* tt = floor (rx) */
6071 scm_num_overflow (s_scm_rationalize
);
6074 SCM_WRONG_TYPE_ARG (1, x
);
6078 /* conversion functions */
6081 scm_is_integer (SCM val
)
6083 return scm_is_true (scm_integer_p (val
));
6087 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6089 if (SCM_I_INUMP (val
))
6091 scm_t_signed_bits n
= SCM_I_INUM (val
);
6092 return n
>= min
&& n
<= max
;
6094 else if (SCM_BIGP (val
))
6096 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6098 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6100 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6102 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6103 return n
>= min
&& n
<= max
;
6113 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6114 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6117 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6118 SCM_I_BIG_MPZ (val
));
6120 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6132 return n
>= min
&& n
<= max
;
6140 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6142 if (SCM_I_INUMP (val
))
6144 scm_t_signed_bits n
= SCM_I_INUM (val
);
6145 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6147 else if (SCM_BIGP (val
))
6149 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6151 else if (max
<= ULONG_MAX
)
6153 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6155 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6156 return n
>= min
&& n
<= max
;
6166 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6169 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6170 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6173 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6174 SCM_I_BIG_MPZ (val
));
6176 return n
>= min
&& n
<= max
;
6184 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6186 scm_error (scm_out_of_range_key
,
6188 "Value out of range ~S to ~S: ~S",
6189 scm_list_3 (min
, max
, bad_val
),
6190 scm_list_1 (bad_val
));
6193 #define TYPE scm_t_intmax
6194 #define TYPE_MIN min
6195 #define TYPE_MAX max
6196 #define SIZEOF_TYPE 0
6197 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6198 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6199 #include "libguile/conv-integer.i.c"
6201 #define TYPE scm_t_uintmax
6202 #define TYPE_MIN min
6203 #define TYPE_MAX max
6204 #define SIZEOF_TYPE 0
6205 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6206 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6207 #include "libguile/conv-uinteger.i.c"
6209 #define TYPE scm_t_int8
6210 #define TYPE_MIN SCM_T_INT8_MIN
6211 #define TYPE_MAX SCM_T_INT8_MAX
6212 #define SIZEOF_TYPE 1
6213 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6214 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6215 #include "libguile/conv-integer.i.c"
6217 #define TYPE scm_t_uint8
6219 #define TYPE_MAX SCM_T_UINT8_MAX
6220 #define SIZEOF_TYPE 1
6221 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6222 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6223 #include "libguile/conv-uinteger.i.c"
6225 #define TYPE scm_t_int16
6226 #define TYPE_MIN SCM_T_INT16_MIN
6227 #define TYPE_MAX SCM_T_INT16_MAX
6228 #define SIZEOF_TYPE 2
6229 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6230 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6231 #include "libguile/conv-integer.i.c"
6233 #define TYPE scm_t_uint16
6235 #define TYPE_MAX SCM_T_UINT16_MAX
6236 #define SIZEOF_TYPE 2
6237 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6238 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6239 #include "libguile/conv-uinteger.i.c"
6241 #define TYPE scm_t_int32
6242 #define TYPE_MIN SCM_T_INT32_MIN
6243 #define TYPE_MAX SCM_T_INT32_MAX
6244 #define SIZEOF_TYPE 4
6245 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6246 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6247 #include "libguile/conv-integer.i.c"
6249 #define TYPE scm_t_uint32
6251 #define TYPE_MAX SCM_T_UINT32_MAX
6252 #define SIZEOF_TYPE 4
6253 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6254 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6255 #include "libguile/conv-uinteger.i.c"
6257 #define TYPE scm_t_wchar
6258 #define TYPE_MIN (scm_t_int32)-1
6259 #define TYPE_MAX (scm_t_int32)0x10ffff
6260 #define SIZEOF_TYPE 4
6261 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6262 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6263 #include "libguile/conv-integer.i.c"
6265 #define TYPE scm_t_int64
6266 #define TYPE_MIN SCM_T_INT64_MIN
6267 #define TYPE_MAX SCM_T_INT64_MAX
6268 #define SIZEOF_TYPE 8
6269 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6270 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6271 #include "libguile/conv-integer.i.c"
6273 #define TYPE scm_t_uint64
6275 #define TYPE_MAX SCM_T_UINT64_MAX
6276 #define SIZEOF_TYPE 8
6277 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6278 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6279 #include "libguile/conv-uinteger.i.c"
6282 scm_to_mpz (SCM val
, mpz_t rop
)
6284 if (SCM_I_INUMP (val
))
6285 mpz_set_si (rop
, SCM_I_INUM (val
));
6286 else if (SCM_BIGP (val
))
6287 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6289 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6293 scm_from_mpz (mpz_t val
)
6295 return scm_i_mpz2num (val
);
6299 scm_is_real (SCM val
)
6301 return scm_is_true (scm_real_p (val
));
6305 scm_is_rational (SCM val
)
6307 return scm_is_true (scm_rational_p (val
));
6311 scm_to_double (SCM val
)
6313 if (SCM_I_INUMP (val
))
6314 return SCM_I_INUM (val
);
6315 else if (SCM_BIGP (val
))
6316 return scm_i_big2dbl (val
);
6317 else if (SCM_FRACTIONP (val
))
6318 return scm_i_fraction2double (val
);
6319 else if (SCM_REALP (val
))
6320 return SCM_REAL_VALUE (val
);
6322 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6326 scm_from_double (double val
)
6330 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
6332 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
6333 SCM_REAL_VALUE (z
) = val
;
6338 #if SCM_ENABLE_DEPRECATED == 1
6341 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6343 scm_c_issue_deprecation_warning
6344 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6348 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6352 scm_out_of_range (NULL
, num
);
6355 return scm_to_double (num
);
6359 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6361 scm_c_issue_deprecation_warning
6362 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6366 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6370 scm_out_of_range (NULL
, num
);
6373 return scm_to_double (num
);
6379 scm_is_complex (SCM val
)
6381 return scm_is_true (scm_complex_p (val
));
6385 scm_c_real_part (SCM z
)
6387 if (SCM_COMPLEXP (z
))
6388 return SCM_COMPLEX_REAL (z
);
6391 /* Use the scm_real_part to get proper error checking and
6394 return scm_to_double (scm_real_part (z
));
6399 scm_c_imag_part (SCM z
)
6401 if (SCM_COMPLEXP (z
))
6402 return SCM_COMPLEX_IMAG (z
);
6405 /* Use the scm_imag_part to get proper error checking and
6406 dispatching. The result will almost always be 0.0, but not
6409 return scm_to_double (scm_imag_part (z
));
6414 scm_c_magnitude (SCM z
)
6416 return scm_to_double (scm_magnitude (z
));
6422 return scm_to_double (scm_angle (z
));
6426 scm_is_number (SCM z
)
6428 return scm_is_true (scm_number_p (z
));
6432 /* In the following functions we dispatch to the real-arg funcs like log()
6433 when we know the arg is real, instead of just handing everything to
6434 clog() for instance. This is in case clog() doesn't optimize for a
6435 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6436 well use it to go straight to the applicable C func. */
6438 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6440 "Return the natural logarithm of @var{z}.")
6441 #define FUNC_NAME s_scm_log
6443 if (SCM_COMPLEXP (z
))
6445 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6446 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6448 double re
= SCM_COMPLEX_REAL (z
);
6449 double im
= SCM_COMPLEX_IMAG (z
);
6450 return scm_c_make_rectangular (log (hypot (re
, im
)),
6456 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6457 although the value itself overflows. */
6458 double re
= scm_to_double (z
);
6459 double l
= log (fabs (re
));
6461 return scm_from_double (l
);
6463 return scm_c_make_rectangular (l
, M_PI
);
6469 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6471 "Return the base 10 logarithm of @var{z}.")
6472 #define FUNC_NAME s_scm_log10
6474 if (SCM_COMPLEXP (z
))
6476 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6477 clog() and a multiply by M_LOG10E, rather than the fallback
6478 log10+hypot+atan2.) */
6479 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6480 && defined SCM_COMPLEX_VALUE
6481 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6483 double re
= SCM_COMPLEX_REAL (z
);
6484 double im
= SCM_COMPLEX_IMAG (z
);
6485 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6486 M_LOG10E
* atan2 (im
, re
));
6491 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6492 although the value itself overflows. */
6493 double re
= scm_to_double (z
);
6494 double l
= log10 (fabs (re
));
6496 return scm_from_double (l
);
6498 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6504 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6506 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6507 "base of natural logarithms (2.71828@dots{}).")
6508 #define FUNC_NAME s_scm_exp
6510 if (SCM_COMPLEXP (z
))
6512 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6513 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6515 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6516 SCM_COMPLEX_IMAG (z
));
6521 /* When z is a negative bignum the conversion to double overflows,
6522 giving -infinity, but that's ok, the exp is still 0.0. */
6523 return scm_from_double (exp (scm_to_double (z
)));
6529 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6531 "Return the square root of @var{z}. Of the two possible roots\n"
6532 "(positive and negative), the one with the a positive real part\n"
6533 "is returned, or if that's zero then a positive imaginary part.\n"
6537 "(sqrt 9.0) @result{} 3.0\n"
6538 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6539 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6540 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6542 #define FUNC_NAME s_scm_sqrt
6544 if (SCM_COMPLEXP (x
))
6546 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6547 && defined SCM_COMPLEX_VALUE
6548 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6550 double re
= SCM_COMPLEX_REAL (x
);
6551 double im
= SCM_COMPLEX_IMAG (x
);
6552 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6553 0.5 * atan2 (im
, re
));
6558 double xx
= scm_to_double (x
);
6560 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6562 return scm_from_double (sqrt (xx
));
6574 mpz_init_set_si (z_negative_one
, -1);
6576 /* It may be possible to tune the performance of some algorithms by using
6577 * the following constants to avoid the creation of bignums. Please, before
6578 * using these values, remember the two rules of program optimization:
6579 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6580 scm_c_define ("most-positive-fixnum",
6581 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6582 scm_c_define ("most-negative-fixnum",
6583 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6585 scm_add_feature ("complex");
6586 scm_add_feature ("inexact");
6587 flo0
= scm_from_double (0.0);
6589 /* determine floating point precision */
6590 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6592 init_dblprec(&scm_dblprec
[i
-2],i
);
6593 init_fx_radix(fx_per_radix
[i
-2],i
);
6596 /* hard code precision for base 10 if the preprocessor tells us to... */
6597 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6600 exactly_one_half
= scm_divide (SCM_I_MAKINUM (1), SCM_I_MAKINUM (2));
6601 #include "libguile/numbers.x"