1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 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
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc() */
54 #include "libguile/_scm.h"
55 #include "libguile/feature.h"
56 #include "libguile/ports.h"
57 #include "libguile/root.h"
58 #include "libguile/smob.h"
59 #include "libguile/strings.h"
61 #include "libguile/validate.h"
62 #include "libguile/numbers.h"
63 #include "libguile/deprecation.h"
65 #include "libguile/eq.h"
67 #include "libguile/discouraged.h"
72 Wonder if this might be faster for some of our code? A switch on
73 the numtag would jump directly to the right case, and the
74 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
76 #define SCM_I_NUMTAG_NOTNUM 0
77 #define SCM_I_NUMTAG_INUM 1
78 #define SCM_I_NUMTAG_BIG scm_tc16_big
79 #define SCM_I_NUMTAG_REAL scm_tc16_real
80 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
81 #define SCM_I_NUMTAG(x) \
82 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
83 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
84 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
85 : SCM_I_NUMTAG_NOTNUM)))
87 /* the macro above will not work as is with fractions */
90 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
92 /* FLOBUFLEN is the maximum number of characters neccessary for the
93 * printed or scm_string representation of an inexact number.
95 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
98 #if ! defined (HAVE_ISNAN)
103 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
106 #if ! defined (HAVE_ISINF)
111 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
118 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
119 an explicit check. In some future gmp (don't know what version number),
120 mpz_cmp_d is supposed to do this itself. */
122 #define xmpz_cmp_d(z, d) \
123 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
125 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
128 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
129 isinf. It does have finite and isnan though, hence the use of those.
130 fpclass would be a possibility on that system too. */
134 #if defined (HAVE_ISINF)
136 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
137 return (! (finite (x
) || isnan (x
)));
146 #if defined (HAVE_ISNAN)
155 static mpz_t z_negative_one
;
159 SCM_C_INLINE_KEYWORD SCM
162 /* Return a newly created bignum. */
163 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
164 mpz_init (SCM_I_BIG_MPZ (z
));
168 SCM_C_INLINE_KEYWORD SCM
169 scm_i_long2big (long x
)
171 /* Return a newly created bignum initialized to X. */
172 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
173 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
177 SCM_C_INLINE_KEYWORD SCM
178 scm_i_ulong2big (unsigned long x
)
180 /* Return a newly created bignum initialized to X. */
181 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
182 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
186 SCM_C_INLINE_KEYWORD
static SCM
187 scm_i_clonebig (SCM src_big
, int same_sign_p
)
189 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
190 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
191 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
193 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
197 SCM_C_INLINE_KEYWORD
int
198 scm_i_bigcmp (SCM x
, SCM y
)
200 /* Return neg if x < y, pos if x > y, and 0 if x == y */
201 /* presume we already know x and y are bignums */
202 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
203 scm_remember_upto_here_2 (x
, y
);
207 SCM_C_INLINE_KEYWORD SCM
208 scm_i_dbl2big (double d
)
210 /* results are only defined if d is an integer */
211 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
212 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
216 /* Convert a integer in double representation to a SCM number. */
218 SCM_C_INLINE_KEYWORD SCM
219 scm_i_dbl2num (double u
)
221 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
222 powers of 2, so there's no rounding when making "double" values
223 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
224 get rounded on a 64-bit machine, hence the "+1".
226 The use of floor() to force to an integer value ensures we get a
227 "numerically closest" value without depending on how a
228 double->long cast or how mpz_set_d will round. For reference,
229 double->long probably follows the hardware rounding mode,
230 mpz_set_d truncates towards zero. */
232 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
233 representable as a double? */
235 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
236 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
237 return SCM_I_MAKINUM ((long) u
);
239 return scm_i_dbl2big (u
);
242 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
243 with R5RS exact->inexact.
245 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
246 (ie. truncate towards zero), then adjust to get the closest double by
247 examining the next lower bit and adding 1 (to the absolute value) if
250 Bignums exactly half way between representable doubles are rounded to the
251 next higher absolute value (ie. away from zero). This seems like an
252 adequate interpretation of R5RS "numerically closest", and it's easier
253 and faster than a full "nearest-even" style.
255 The bit test must be done on the absolute value of the mpz_t, which means
256 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
257 negatives as twos complement.
259 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
260 following the hardware rounding mode, but applied to the absolute value
261 of the mpz_t operand. This is not what we want so we put the high
262 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
263 mpz_get_d is supposed to always truncate towards zero.
265 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
266 is a slowdown. It'd be faster to pick out the relevant high bits with
267 mpz_getlimbn if we could be bothered coding that, and if the new
268 truncating gmp doesn't come out. */
271 scm_i_big2dbl (SCM b
)
276 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
280 /* Current GMP, eg. 4.1.3, force truncation towards zero */
282 if (bits
> DBL_MANT_DIG
)
284 size_t shift
= bits
- DBL_MANT_DIG
;
285 mpz_init2 (tmp
, DBL_MANT_DIG
);
286 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
287 result
= ldexp (mpz_get_d (tmp
), shift
);
292 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
297 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
300 if (bits
> DBL_MANT_DIG
)
302 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
303 /* test bit number "pos" in absolute value */
304 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
305 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
307 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
311 scm_remember_upto_here_1 (b
);
315 SCM_C_INLINE_KEYWORD SCM
316 scm_i_normbig (SCM b
)
318 /* convert a big back to a fixnum if it'll fit */
319 /* presume b is a bignum */
320 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
322 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
323 if (SCM_FIXABLE (val
))
324 b
= SCM_I_MAKINUM (val
);
329 static SCM_C_INLINE_KEYWORD SCM
330 scm_i_mpz2num (mpz_t b
)
332 /* convert a mpz number to a SCM number. */
333 if (mpz_fits_slong_p (b
))
335 long val
= mpz_get_si (b
);
336 if (SCM_FIXABLE (val
))
337 return SCM_I_MAKINUM (val
);
341 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
342 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
347 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
348 static SCM
scm_divide2real (SCM x
, SCM y
);
351 scm_i_make_ratio (SCM numerator
, SCM denominator
)
352 #define FUNC_NAME "make-ratio"
354 /* First make sure the arguments are proper.
356 if (SCM_I_INUMP (denominator
))
358 if (scm_is_eq (denominator
, SCM_INUM0
))
359 scm_num_overflow ("make-ratio");
360 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
365 if (!(SCM_BIGP(denominator
)))
366 SCM_WRONG_TYPE_ARG (2, denominator
);
368 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
369 SCM_WRONG_TYPE_ARG (1, numerator
);
371 /* Then flip signs so that the denominator is positive.
373 if (scm_is_true (scm_negative_p (denominator
)))
375 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
376 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
379 /* Now consider for each of the four fixnum/bignum combinations
380 whether the rational number is really an integer.
382 if (SCM_I_INUMP (numerator
))
384 long x
= SCM_I_INUM (numerator
);
385 if (scm_is_eq (numerator
, SCM_INUM0
))
387 if (SCM_I_INUMP (denominator
))
390 y
= SCM_I_INUM (denominator
);
392 return SCM_I_MAKINUM(1);
394 return SCM_I_MAKINUM (x
/ y
);
398 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
399 of that value for the denominator, as a bignum. Apart from
400 that case, abs(bignum) > abs(inum) so inum/bignum is not an
402 if (x
== SCM_MOST_NEGATIVE_FIXNUM
403 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
404 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
405 return SCM_I_MAKINUM(-1);
408 else if (SCM_BIGP (numerator
))
410 if (SCM_I_INUMP (denominator
))
412 long yy
= SCM_I_INUM (denominator
);
413 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
414 return scm_divide (numerator
, denominator
);
418 if (scm_is_eq (numerator
, denominator
))
419 return SCM_I_MAKINUM(1);
420 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
421 SCM_I_BIG_MPZ (denominator
)))
422 return scm_divide(numerator
, denominator
);
426 /* No, it's a proper fraction.
428 return scm_double_cell (scm_tc16_fraction
,
429 SCM_UNPACK (numerator
),
430 SCM_UNPACK (denominator
), 0);
434 static void scm_i_fraction_reduce (SCM z
)
436 if (!(SCM_FRACTION_REDUCED (z
)))
439 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
440 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
443 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
444 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
446 SCM_FRACTION_REDUCED_SET (z
);
451 scm_i_fraction2double (SCM z
)
453 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
454 SCM_FRACTION_DENOMINATOR (z
)));
457 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
459 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
461 #define FUNC_NAME s_scm_exact_p
467 if (SCM_FRACTIONP (x
))
471 SCM_WRONG_TYPE_ARG (1, x
);
476 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
478 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
480 #define FUNC_NAME s_scm_odd_p
484 long val
= SCM_I_INUM (n
);
485 return scm_from_bool ((val
& 1L) != 0);
487 else if (SCM_BIGP (n
))
489 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
490 scm_remember_upto_here_1 (n
);
491 return scm_from_bool (odd_p
);
493 else if (scm_is_true (scm_inf_p (n
)))
495 else if (SCM_REALP (n
))
497 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
503 SCM_WRONG_TYPE_ARG (1, n
);
506 SCM_WRONG_TYPE_ARG (1, n
);
511 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
513 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
515 #define FUNC_NAME s_scm_even_p
519 long val
= SCM_I_INUM (n
);
520 return scm_from_bool ((val
& 1L) == 0);
522 else if (SCM_BIGP (n
))
524 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
525 scm_remember_upto_here_1 (n
);
526 return scm_from_bool (even_p
);
528 else if (scm_is_true (scm_inf_p (n
)))
530 else if (SCM_REALP (n
))
532 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
538 SCM_WRONG_TYPE_ARG (1, n
);
541 SCM_WRONG_TYPE_ARG (1, n
);
545 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
547 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
548 "or @samp{-inf.0}, @code{#f} otherwise.")
549 #define FUNC_NAME s_scm_inf_p
552 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
553 else if (SCM_COMPLEXP (x
))
554 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
555 || xisinf (SCM_COMPLEX_IMAG (x
)));
561 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
563 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
565 #define FUNC_NAME s_scm_nan_p
568 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
569 else if (SCM_COMPLEXP (n
))
570 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
571 || xisnan (SCM_COMPLEX_IMAG (n
)));
577 /* Guile's idea of infinity. */
578 static double guile_Inf
;
580 /* Guile's idea of not a number. */
581 static double guile_NaN
;
584 guile_ieee_init (void)
586 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
588 /* Some version of gcc on some old version of Linux used to crash when
589 trying to make Inf and NaN. */
592 /* C99 INFINITY, when available.
593 FIXME: The standard allows for INFINITY to be something that overflows
594 at compile time. We ought to have a configure test to check for that
595 before trying to use it. (But in practice we believe this is not a
596 problem on any system guile is likely to target.) */
597 guile_Inf
= INFINITY
;
600 extern unsigned int DINFINITY
[2];
601 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
608 if (guile_Inf
== tmp
)
616 #if defined (HAVE_ISNAN)
619 /* C99 NAN, when available */
623 extern unsigned int DQNAN
[2];
624 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
626 guile_NaN
= guile_Inf
/ guile_Inf
;
632 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
635 #define FUNC_NAME s_scm_inf
637 static int initialized
= 0;
643 return scm_from_double (guile_Inf
);
647 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
650 #define FUNC_NAME s_scm_nan
652 static int initialized
= 0;
658 return scm_from_double (guile_NaN
);
663 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
665 "Return the absolute value of @var{x}.")
670 long int xx
= SCM_I_INUM (x
);
673 else if (SCM_POSFIXABLE (-xx
))
674 return SCM_I_MAKINUM (-xx
);
676 return scm_i_long2big (-xx
);
678 else if (SCM_BIGP (x
))
680 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
682 return scm_i_clonebig (x
, 0);
686 else if (SCM_REALP (x
))
688 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
689 double xx
= SCM_REAL_VALUE (x
);
691 return scm_from_double (-xx
);
695 else if (SCM_FRACTIONP (x
))
697 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
699 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
700 SCM_FRACTION_DENOMINATOR (x
));
703 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
708 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
709 /* "Return the quotient of the numbers @var{x} and @var{y}."
712 scm_quotient (SCM x
, SCM y
)
716 long xx
= SCM_I_INUM (x
);
719 long yy
= SCM_I_INUM (y
);
721 scm_num_overflow (s_quotient
);
726 return SCM_I_MAKINUM (z
);
728 return scm_i_long2big (z
);
731 else if (SCM_BIGP (y
))
733 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
734 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
735 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
737 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
738 scm_remember_upto_here_1 (y
);
739 return SCM_I_MAKINUM (-1);
742 return SCM_I_MAKINUM (0);
745 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
747 else if (SCM_BIGP (x
))
751 long yy
= SCM_I_INUM (y
);
753 scm_num_overflow (s_quotient
);
758 SCM result
= scm_i_mkbig ();
761 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
764 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
767 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
768 scm_remember_upto_here_1 (x
);
769 return scm_i_normbig (result
);
772 else if (SCM_BIGP (y
))
774 SCM result
= scm_i_mkbig ();
775 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
778 scm_remember_upto_here_2 (x
, y
);
779 return scm_i_normbig (result
);
782 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
785 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
788 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
789 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
791 * "(remainder 13 4) @result{} 1\n"
792 * "(remainder -13 4) @result{} -1\n"
796 scm_remainder (SCM x
, SCM y
)
802 long yy
= SCM_I_INUM (y
);
804 scm_num_overflow (s_remainder
);
807 long z
= SCM_I_INUM (x
) % yy
;
808 return SCM_I_MAKINUM (z
);
811 else if (SCM_BIGP (y
))
813 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
814 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
815 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
817 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
818 scm_remember_upto_here_1 (y
);
819 return SCM_I_MAKINUM (0);
825 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
827 else if (SCM_BIGP (x
))
831 long yy
= SCM_I_INUM (y
);
833 scm_num_overflow (s_remainder
);
836 SCM result
= scm_i_mkbig ();
839 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
840 scm_remember_upto_here_1 (x
);
841 return scm_i_normbig (result
);
844 else if (SCM_BIGP (y
))
846 SCM result
= scm_i_mkbig ();
847 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
850 scm_remember_upto_here_2 (x
, y
);
851 return scm_i_normbig (result
);
854 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
857 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
861 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
862 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
864 * "(modulo 13 4) @result{} 1\n"
865 * "(modulo -13 4) @result{} 3\n"
869 scm_modulo (SCM x
, SCM y
)
873 long xx
= SCM_I_INUM (x
);
876 long yy
= SCM_I_INUM (y
);
878 scm_num_overflow (s_modulo
);
881 /* C99 specifies that "%" is the remainder corresponding to a
882 quotient rounded towards zero, and that's also traditional
883 for machine division, so z here should be well defined. */
901 return SCM_I_MAKINUM (result
);
904 else if (SCM_BIGP (y
))
906 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
913 SCM pos_y
= scm_i_clonebig (y
, 0);
914 /* do this after the last scm_op */
915 mpz_init_set_si (z_x
, xx
);
916 result
= pos_y
; /* re-use this bignum */
917 mpz_mod (SCM_I_BIG_MPZ (result
),
919 SCM_I_BIG_MPZ (pos_y
));
920 scm_remember_upto_here_1 (pos_y
);
924 result
= scm_i_mkbig ();
925 /* do this after the last scm_op */
926 mpz_init_set_si (z_x
, xx
);
927 mpz_mod (SCM_I_BIG_MPZ (result
),
930 scm_remember_upto_here_1 (y
);
933 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
934 mpz_add (SCM_I_BIG_MPZ (result
),
936 SCM_I_BIG_MPZ (result
));
937 scm_remember_upto_here_1 (y
);
938 /* and do this before the next one */
940 return scm_i_normbig (result
);
944 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
946 else if (SCM_BIGP (x
))
950 long yy
= SCM_I_INUM (y
);
952 scm_num_overflow (s_modulo
);
955 SCM result
= scm_i_mkbig ();
956 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
958 (yy
< 0) ? - yy
: yy
);
959 scm_remember_upto_here_1 (x
);
960 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
961 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
962 SCM_I_BIG_MPZ (result
),
964 return scm_i_normbig (result
);
967 else if (SCM_BIGP (y
))
970 SCM result
= scm_i_mkbig ();
971 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
972 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
973 mpz_mod (SCM_I_BIG_MPZ (result
),
975 SCM_I_BIG_MPZ (pos_y
));
977 scm_remember_upto_here_1 (x
);
978 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
979 mpz_add (SCM_I_BIG_MPZ (result
),
981 SCM_I_BIG_MPZ (result
));
982 scm_remember_upto_here_2 (y
, pos_y
);
983 return scm_i_normbig (result
);
987 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
990 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
993 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
994 /* "Return the greatest common divisor of all arguments.\n"
995 * "If called without arguments, 0 is returned."
998 scm_gcd (SCM x
, SCM y
)
1001 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1003 if (SCM_I_INUMP (x
))
1005 if (SCM_I_INUMP (y
))
1007 long xx
= SCM_I_INUM (x
);
1008 long yy
= SCM_I_INUM (y
);
1009 long u
= xx
< 0 ? -xx
: xx
;
1010 long v
= yy
< 0 ? -yy
: yy
;
1020 /* Determine a common factor 2^k */
1021 while (!(1 & (u
| v
)))
1027 /* Now, any factor 2^n can be eliminated */
1047 return (SCM_POSFIXABLE (result
)
1048 ? SCM_I_MAKINUM (result
)
1049 : scm_i_long2big (result
));
1051 else if (SCM_BIGP (y
))
1057 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1059 else if (SCM_BIGP (x
))
1061 if (SCM_I_INUMP (y
))
1063 unsigned long result
;
1066 yy
= SCM_I_INUM (y
);
1071 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1072 scm_remember_upto_here_1 (x
);
1073 return (SCM_POSFIXABLE (result
)
1074 ? SCM_I_MAKINUM (result
)
1075 : scm_from_ulong (result
));
1077 else if (SCM_BIGP (y
))
1079 SCM result
= scm_i_mkbig ();
1080 mpz_gcd (SCM_I_BIG_MPZ (result
),
1083 scm_remember_upto_here_2 (x
, y
);
1084 return scm_i_normbig (result
);
1087 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1090 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1093 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1094 /* "Return the least common multiple of the arguments.\n"
1095 * "If called without arguments, 1 is returned."
1098 scm_lcm (SCM n1
, SCM n2
)
1100 if (SCM_UNBNDP (n2
))
1102 if (SCM_UNBNDP (n1
))
1103 return SCM_I_MAKINUM (1L);
1104 n2
= SCM_I_MAKINUM (1L);
1107 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1108 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1109 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1110 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1112 if (SCM_I_INUMP (n1
))
1114 if (SCM_I_INUMP (n2
))
1116 SCM d
= scm_gcd (n1
, n2
);
1117 if (scm_is_eq (d
, SCM_INUM0
))
1120 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1124 /* inum n1, big n2 */
1127 SCM result
= scm_i_mkbig ();
1128 long nn1
= SCM_I_INUM (n1
);
1129 if (nn1
== 0) return SCM_INUM0
;
1130 if (nn1
< 0) nn1
= - nn1
;
1131 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1132 scm_remember_upto_here_1 (n2
);
1140 if (SCM_I_INUMP (n2
))
1147 SCM result
= scm_i_mkbig ();
1148 mpz_lcm(SCM_I_BIG_MPZ (result
),
1150 SCM_I_BIG_MPZ (n2
));
1151 scm_remember_upto_here_2(n1
, n2
);
1152 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1158 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1163 + + + x (map digit:logand X Y)
1164 + - + x (map digit:logand X (lognot (+ -1 Y)))
1165 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1166 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1171 + + + (map digit:logior X Y)
1172 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1173 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1174 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1179 + + + (map digit:logxor X Y)
1180 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1181 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1182 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1187 + + (any digit:logand X Y)
1188 + - (any digit:logand X (lognot (+ -1 Y)))
1189 - + (any digit:logand (lognot (+ -1 X)) Y)
1194 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1196 "Return the bitwise AND of the integer arguments.\n\n"
1198 "(logand) @result{} -1\n"
1199 "(logand 7) @result{} 7\n"
1200 "(logand #b111 #b011 #b001) @result{} 1\n"
1202 #define FUNC_NAME s_scm_logand
1206 if (SCM_UNBNDP (n2
))
1208 if (SCM_UNBNDP (n1
))
1209 return SCM_I_MAKINUM (-1);
1210 else if (!SCM_NUMBERP (n1
))
1211 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1212 else if (SCM_NUMBERP (n1
))
1215 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1218 if (SCM_I_INUMP (n1
))
1220 nn1
= SCM_I_INUM (n1
);
1221 if (SCM_I_INUMP (n2
))
1223 long nn2
= SCM_I_INUM (n2
);
1224 return SCM_I_MAKINUM (nn1
& nn2
);
1226 else if SCM_BIGP (n2
)
1232 SCM result_z
= scm_i_mkbig ();
1234 mpz_init_set_si (nn1_z
, nn1
);
1235 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1236 scm_remember_upto_here_1 (n2
);
1238 return scm_i_normbig (result_z
);
1242 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1244 else if (SCM_BIGP (n1
))
1246 if (SCM_I_INUMP (n2
))
1249 nn1
= SCM_I_INUM (n1
);
1252 else if (SCM_BIGP (n2
))
1254 SCM result_z
= scm_i_mkbig ();
1255 mpz_and (SCM_I_BIG_MPZ (result_z
),
1257 SCM_I_BIG_MPZ (n2
));
1258 scm_remember_upto_here_2 (n1
, n2
);
1259 return scm_i_normbig (result_z
);
1262 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1265 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1270 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1272 "Return the bitwise OR of the integer arguments.\n\n"
1274 "(logior) @result{} 0\n"
1275 "(logior 7) @result{} 7\n"
1276 "(logior #b000 #b001 #b011) @result{} 3\n"
1278 #define FUNC_NAME s_scm_logior
1282 if (SCM_UNBNDP (n2
))
1284 if (SCM_UNBNDP (n1
))
1286 else if (SCM_NUMBERP (n1
))
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1292 if (SCM_I_INUMP (n1
))
1294 nn1
= SCM_I_INUM (n1
);
1295 if (SCM_I_INUMP (n2
))
1297 long nn2
= SCM_I_INUM (n2
);
1298 return SCM_I_MAKINUM (nn1
| nn2
);
1300 else if (SCM_BIGP (n2
))
1306 SCM result_z
= scm_i_mkbig ();
1308 mpz_init_set_si (nn1_z
, nn1
);
1309 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1310 scm_remember_upto_here_1 (n2
);
1316 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1318 else if (SCM_BIGP (n1
))
1320 if (SCM_I_INUMP (n2
))
1323 nn1
= SCM_I_INUM (n1
);
1326 else if (SCM_BIGP (n2
))
1328 SCM result_z
= scm_i_mkbig ();
1329 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1331 SCM_I_BIG_MPZ (n2
));
1332 scm_remember_upto_here_2 (n1
, n2
);
1336 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1339 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1344 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1346 "Return the bitwise XOR of the integer arguments. A bit is\n"
1347 "set in the result if it is set in an odd number of arguments.\n"
1349 "(logxor) @result{} 0\n"
1350 "(logxor 7) @result{} 7\n"
1351 "(logxor #b000 #b001 #b011) @result{} 2\n"
1352 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1354 #define FUNC_NAME s_scm_logxor
1358 if (SCM_UNBNDP (n2
))
1360 if (SCM_UNBNDP (n1
))
1362 else if (SCM_NUMBERP (n1
))
1365 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1368 if (SCM_I_INUMP (n1
))
1370 nn1
= SCM_I_INUM (n1
);
1371 if (SCM_I_INUMP (n2
))
1373 long nn2
= SCM_I_INUM (n2
);
1374 return SCM_I_MAKINUM (nn1
^ nn2
);
1376 else if (SCM_BIGP (n2
))
1380 SCM result_z
= scm_i_mkbig ();
1382 mpz_init_set_si (nn1_z
, nn1
);
1383 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1384 scm_remember_upto_here_1 (n2
);
1386 return scm_i_normbig (result_z
);
1390 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1392 else if (SCM_BIGP (n1
))
1394 if (SCM_I_INUMP (n2
))
1397 nn1
= SCM_I_INUM (n1
);
1400 else if (SCM_BIGP (n2
))
1402 SCM result_z
= scm_i_mkbig ();
1403 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1405 SCM_I_BIG_MPZ (n2
));
1406 scm_remember_upto_here_2 (n1
, n2
);
1407 return scm_i_normbig (result_z
);
1410 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1413 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1418 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1421 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1422 "(logtest #b0100 #b1011) @result{} #f\n"
1423 "(logtest #b0100 #b0111) @result{} #t\n"
1425 #define FUNC_NAME s_scm_logtest
1429 if (SCM_I_INUMP (j
))
1431 nj
= SCM_I_INUM (j
);
1432 if (SCM_I_INUMP (k
))
1434 long nk
= SCM_I_INUM (k
);
1435 return scm_from_bool (nj
& nk
);
1437 else if (SCM_BIGP (k
))
1445 mpz_init_set_si (nj_z
, nj
);
1446 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1447 scm_remember_upto_here_1 (k
);
1448 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1454 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1456 else if (SCM_BIGP (j
))
1458 if (SCM_I_INUMP (k
))
1461 nj
= SCM_I_INUM (j
);
1464 else if (SCM_BIGP (k
))
1468 mpz_init (result_z
);
1472 scm_remember_upto_here_2 (j
, k
);
1473 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1474 mpz_clear (result_z
);
1478 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1481 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1486 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1489 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1490 "(logbit? 0 #b1101) @result{} #t\n"
1491 "(logbit? 1 #b1101) @result{} #f\n"
1492 "(logbit? 2 #b1101) @result{} #t\n"
1493 "(logbit? 3 #b1101) @result{} #t\n"
1494 "(logbit? 4 #b1101) @result{} #f\n"
1496 #define FUNC_NAME s_scm_logbit_p
1498 unsigned long int iindex
;
1499 iindex
= scm_to_ulong (index
);
1501 if (SCM_I_INUMP (j
))
1503 /* bits above what's in an inum follow the sign bit */
1504 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1505 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1507 else if (SCM_BIGP (j
))
1509 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1510 scm_remember_upto_here_1 (j
);
1511 return scm_from_bool (val
);
1514 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1519 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1521 "Return the integer which is the ones-complement of the integer\n"
1525 "(number->string (lognot #b10000000) 2)\n"
1526 " @result{} \"-10000001\"\n"
1527 "(number->string (lognot #b0) 2)\n"
1528 " @result{} \"-1\"\n"
1530 #define FUNC_NAME s_scm_lognot
1532 if (SCM_I_INUMP (n
)) {
1533 /* No overflow here, just need to toggle all the bits making up the inum.
1534 Enhancement: No need to strip the tag and add it back, could just xor
1535 a block of 1 bits, if that worked with the various debug versions of
1537 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1539 } else if (SCM_BIGP (n
)) {
1540 SCM result
= scm_i_mkbig ();
1541 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1542 scm_remember_upto_here_1 (n
);
1546 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1551 /* returns 0 if IN is not an integer. OUT must already be
1554 coerce_to_big (SCM in
, mpz_t out
)
1557 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1558 else if (SCM_I_INUMP (in
))
1559 mpz_set_si (out
, SCM_I_INUM (in
));
1566 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1567 (SCM n
, SCM k
, SCM m
),
1568 "Return @var{n} raised to the integer exponent\n"
1569 "@var{k}, modulo @var{m}.\n"
1572 "(modulo-expt 2 3 5)\n"
1575 #define FUNC_NAME s_scm_modulo_expt
1581 /* There are two classes of error we might encounter --
1582 1) Math errors, which we'll report by calling scm_num_overflow,
1584 2) wrong-type errors, which of course we'll report by calling
1586 We don't report those errors immediately, however; instead we do
1587 some cleanup first. These variables tell us which error (if
1588 any) we should report after cleaning up.
1590 int report_overflow
= 0;
1592 int position_of_wrong_type
= 0;
1593 SCM value_of_wrong_type
= SCM_INUM0
;
1595 SCM result
= SCM_UNDEFINED
;
1601 if (scm_is_eq (m
, SCM_INUM0
))
1603 report_overflow
= 1;
1607 if (!coerce_to_big (n
, n_tmp
))
1609 value_of_wrong_type
= n
;
1610 position_of_wrong_type
= 1;
1614 if (!coerce_to_big (k
, k_tmp
))
1616 value_of_wrong_type
= k
;
1617 position_of_wrong_type
= 2;
1621 if (!coerce_to_big (m
, m_tmp
))
1623 value_of_wrong_type
= m
;
1624 position_of_wrong_type
= 3;
1628 /* if the exponent K is negative, and we simply call mpz_powm, we
1629 will get a divide-by-zero exception when an inverse 1/n mod m
1630 doesn't exist (or is not unique). Since exceptions are hard to
1631 handle, we'll attempt the inversion "by hand" -- that way, we get
1632 a simple failure code, which is easy to handle. */
1634 if (-1 == mpz_sgn (k_tmp
))
1636 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1638 report_overflow
= 1;
1641 mpz_neg (k_tmp
, k_tmp
);
1644 result
= scm_i_mkbig ();
1645 mpz_powm (SCM_I_BIG_MPZ (result
),
1650 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1651 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1658 if (report_overflow
)
1659 scm_num_overflow (FUNC_NAME
);
1661 if (position_of_wrong_type
)
1662 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1663 value_of_wrong_type
);
1665 return scm_i_normbig (result
);
1669 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1671 "Return @var{n} raised to the exact integer exponent\n"
1675 "(integer-expt 2 5)\n"
1677 "(integer-expt -3 3)\n"
1680 #define FUNC_NAME s_scm_integer_expt
1683 SCM z_i2
= SCM_BOOL_F
;
1685 SCM acc
= SCM_I_MAKINUM (1L);
1687 /* 0^0 == 1 according to R5RS */
1688 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1689 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1690 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1691 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1693 if (SCM_I_INUMP (k
))
1694 i2
= SCM_I_INUM (k
);
1695 else if (SCM_BIGP (k
))
1697 z_i2
= scm_i_clonebig (k
, 1);
1698 scm_remember_upto_here_1 (k
);
1702 SCM_WRONG_TYPE_ARG (2, k
);
1706 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1708 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1709 n
= scm_divide (n
, SCM_UNDEFINED
);
1713 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1717 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1719 return scm_product (acc
, n
);
1721 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1722 acc
= scm_product (acc
, n
);
1723 n
= scm_product (n
, n
);
1724 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1732 n
= scm_divide (n
, SCM_UNDEFINED
);
1739 return scm_product (acc
, n
);
1741 acc
= scm_product (acc
, n
);
1742 n
= scm_product (n
, n
);
1749 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1751 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1752 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1754 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1755 "@var{cnt} is negative it's a division, rounded towards negative\n"
1756 "infinity. (Note that this is not the same rounding as\n"
1757 "@code{quotient} does.)\n"
1759 "With @var{n} viewed as an infinite precision twos complement,\n"
1760 "@code{ash} means a left shift introducing zero bits, or a right\n"
1761 "shift dropping bits.\n"
1764 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1765 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1767 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1768 "(ash -23 -2) @result{} -6\n"
1770 #define FUNC_NAME s_scm_ash
1773 bits_to_shift
= scm_to_long (cnt
);
1775 if (bits_to_shift
< 0)
1777 /* Shift right by abs(cnt) bits. This is realized as a division
1778 by div:=2^abs(cnt). However, to guarantee the floor
1779 rounding, negative values require some special treatment.
1781 SCM div
= scm_integer_expt (SCM_I_MAKINUM (2),
1782 scm_from_long (-bits_to_shift
));
1784 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1785 if (scm_is_false (scm_negative_p (n
)))
1786 return scm_quotient (n
, div
);
1788 return scm_sum (SCM_I_MAKINUM (-1L),
1789 scm_quotient (scm_sum (SCM_I_MAKINUM (1L), n
), div
));
1792 /* Shift left is done by multiplication with 2^CNT */
1793 return scm_product (n
, scm_integer_expt (SCM_I_MAKINUM (2), cnt
));
1798 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1799 (SCM n
, SCM start
, SCM end
),
1800 "Return the integer composed of the @var{start} (inclusive)\n"
1801 "through @var{end} (exclusive) bits of @var{n}. The\n"
1802 "@var{start}th bit becomes the 0-th bit in the result.\n"
1805 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1806 " @result{} \"1010\"\n"
1807 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1808 " @result{} \"10110\"\n"
1810 #define FUNC_NAME s_scm_bit_extract
1812 unsigned long int istart
, iend
, bits
;
1813 istart
= scm_to_ulong (start
);
1814 iend
= scm_to_ulong (end
);
1815 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1817 /* how many bits to keep */
1818 bits
= iend
- istart
;
1820 if (SCM_I_INUMP (n
))
1822 long int in
= SCM_I_INUM (n
);
1824 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1825 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1826 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1828 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1830 /* Since we emulate two's complement encoded numbers, this
1831 * special case requires us to produce a result that has
1832 * more bits than can be stored in a fixnum.
1834 SCM result
= scm_i_long2big (in
);
1835 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1840 /* mask down to requisite bits */
1841 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1842 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1844 else if (SCM_BIGP (n
))
1849 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1853 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1854 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1855 such bits into a ulong. */
1856 result
= scm_i_mkbig ();
1857 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1858 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1859 result
= scm_i_normbig (result
);
1861 scm_remember_upto_here_1 (n
);
1865 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1870 static const char scm_logtab
[] = {
1871 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1874 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1876 "Return the number of bits in integer @var{n}. If integer is\n"
1877 "positive, the 1-bits in its binary representation are counted.\n"
1878 "If negative, the 0-bits in its two's-complement binary\n"
1879 "representation are counted. If 0, 0 is returned.\n"
1882 "(logcount #b10101010)\n"
1889 #define FUNC_NAME s_scm_logcount
1891 if (SCM_I_INUMP (n
))
1893 unsigned long int c
= 0;
1894 long int nn
= SCM_I_INUM (n
);
1899 c
+= scm_logtab
[15 & nn
];
1902 return SCM_I_MAKINUM (c
);
1904 else if (SCM_BIGP (n
))
1906 unsigned long count
;
1907 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1908 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1910 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1911 scm_remember_upto_here_1 (n
);
1912 return SCM_I_MAKINUM (count
);
1915 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1920 static const char scm_ilentab
[] = {
1921 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1925 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1927 "Return the number of bits necessary to represent @var{n}.\n"
1930 "(integer-length #b10101010)\n"
1932 "(integer-length 0)\n"
1934 "(integer-length #b1111)\n"
1937 #define FUNC_NAME s_scm_integer_length
1939 if (SCM_I_INUMP (n
))
1941 unsigned long int c
= 0;
1943 long int nn
= SCM_I_INUM (n
);
1949 l
= scm_ilentab
[15 & nn
];
1952 return SCM_I_MAKINUM (c
- 4 + l
);
1954 else if (SCM_BIGP (n
))
1956 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1957 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1958 1 too big, so check for that and adjust. */
1959 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1960 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1961 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1962 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1964 scm_remember_upto_here_1 (n
);
1965 return SCM_I_MAKINUM (size
);
1968 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1972 /*** NUMBERS -> STRINGS ***/
1973 #define SCM_MAX_DBL_PREC 60
1974 #define SCM_MAX_DBL_RADIX 36
1976 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1977 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1978 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1981 void init_dblprec(int *prec
, int radix
) {
1982 /* determine floating point precision by adding successively
1983 smaller increments to 1.0 until it is considered == 1.0 */
1984 double f
= ((double)1.0)/radix
;
1985 double fsum
= 1.0 + f
;
1990 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2002 void init_fx_radix(double *fx_list
, int radix
)
2004 /* initialize a per-radix list of tolerances. When added
2005 to a number < 1.0, we can determine if we should raund
2006 up and quit converting a number to a string. */
2010 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2011 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2014 /* use this array as a way to generate a single digit */
2015 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2018 idbl2str (double f
, char *a
, int radix
)
2020 int efmt
, dpt
, d
, i
, wp
;
2022 #ifdef DBL_MIN_10_EXP
2025 #endif /* DBL_MIN_10_EXP */
2030 radix
> SCM_MAX_DBL_RADIX
)
2032 /* revert to existing behavior */
2036 wp
= scm_dblprec
[radix
-2];
2037 fx
= fx_per_radix
[radix
-2];
2041 #ifdef HAVE_COPYSIGN
2042 double sgn
= copysign (1.0, f
);
2047 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2053 strcpy (a
, "-inf.0");
2055 strcpy (a
, "+inf.0");
2058 else if (xisnan (f
))
2060 strcpy (a
, "+nan.0");
2070 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2071 make-uniform-vector, from causing infinite loops. */
2072 /* just do the checking...if it passes, we do the conversion for our
2073 radix again below */
2080 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2088 while (f_cpy
> 10.0)
2091 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2112 if (f
+ fx
[wp
] >= radix
)
2119 /* adding 9999 makes this equivalent to abs(x) % 3 */
2120 dpt
= (exp
+ 9999) % 3;
2124 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2146 a
[ch
++] = number_chars
[d
];
2149 if (f
+ fx
[wp
] >= 1.0)
2151 a
[ch
- 1] = number_chars
[d
+1];
2163 if ((dpt
> 4) && (exp
> 6))
2165 d
= (a
[0] == '-' ? 2 : 1);
2166 for (i
= ch
++; i
> d
; i
--)
2179 if (a
[ch
- 1] == '.')
2180 a
[ch
++] = '0'; /* trailing zero */
2189 for (i
= radix
; i
<= exp
; i
*= radix
);
2190 for (i
/= radix
; i
; i
/= radix
)
2192 a
[ch
++] = number_chars
[exp
/ i
];
2201 icmplx2str (double real
, double imag
, char *str
, int radix
)
2205 i
= idbl2str (real
, str
, radix
);
2208 /* Don't output a '+' for negative numbers or for Inf and
2209 NaN. They will provide their own sign. */
2210 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2212 i
+= idbl2str (imag
, &str
[i
], radix
);
2219 iflo2str (SCM flt
, char *str
, int radix
)
2222 if (SCM_REALP (flt
))
2223 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2225 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2230 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2231 characters in the result.
2233 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2235 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2240 return scm_iuint2str (-num
, rad
, p
) + 1;
2243 return scm_iuint2str (num
, rad
, p
);
2246 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2247 characters in the result.
2249 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2251 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2255 scm_t_uintmax n
= num
;
2257 for (n
/= rad
; n
> 0; n
/= rad
)
2267 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2272 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2274 "Return a string holding the external representation of the\n"
2275 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2276 "inexact, a radix of 10 will be used.")
2277 #define FUNC_NAME s_scm_number_to_string
2281 if (SCM_UNBNDP (radix
))
2284 base
= scm_to_signed_integer (radix
, 2, 36);
2286 if (SCM_I_INUMP (n
))
2288 char num_buf
[SCM_INTBUFLEN
];
2289 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2290 return scm_from_locale_stringn (num_buf
, length
);
2292 else if (SCM_BIGP (n
))
2294 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2295 scm_remember_upto_here_1 (n
);
2296 return scm_take_locale_string (str
);
2298 else if (SCM_FRACTIONP (n
))
2300 scm_i_fraction_reduce (n
);
2301 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2302 scm_from_locale_string ("/"),
2303 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2305 else if (SCM_INEXACTP (n
))
2307 char num_buf
[FLOBUFLEN
];
2308 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2311 SCM_WRONG_TYPE_ARG (1, n
);
2316 /* These print routines used to be stubbed here so that scm_repl.c
2317 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2320 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2322 char num_buf
[FLOBUFLEN
];
2323 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2328 scm_i_print_double (double val
, SCM port
)
2330 char num_buf
[FLOBUFLEN
];
2331 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2335 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2338 char num_buf
[FLOBUFLEN
];
2339 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2344 scm_i_print_complex (double real
, double imag
, SCM port
)
2346 char num_buf
[FLOBUFLEN
];
2347 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2351 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2354 scm_i_fraction_reduce (sexp
);
2355 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2356 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2357 scm_remember_upto_here_1 (str
);
2362 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2364 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2365 scm_remember_upto_here_1 (exp
);
2366 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2370 /*** END nums->strs ***/
2373 /*** STRINGS -> NUMBERS ***/
2375 /* The following functions implement the conversion from strings to numbers.
2376 * The implementation somehow follows the grammar for numbers as it is given
2377 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2378 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2379 * points should be noted about the implementation:
2380 * * Each function keeps a local index variable 'idx' that points at the
2381 * current position within the parsed string. The global index is only
2382 * updated if the function could parse the corresponding syntactic unit
2384 * * Similarly, the functions keep track of indicators of inexactness ('#',
2385 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2386 * global exactness information is only updated after each part has been
2387 * successfully parsed.
2388 * * Sequences of digits are parsed into temporary variables holding fixnums.
2389 * Only if these fixnums would overflow, the result variables are updated
2390 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2391 * the temporary variables holding the fixnums are cleared, and the process
2392 * starts over again. If for example fixnums were able to store five decimal
2393 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2394 * and the result was computed as 12345 * 100000 + 67890. In other words,
2395 * only every five digits two bignum operations were performed.
2398 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2400 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2402 /* In non ASCII-style encodings the following macro might not work. */
2403 #define XDIGIT2UINT(d) \
2404 (isdigit ((int) (unsigned char) d) \
2406 : tolower ((int) (unsigned char) d) - 'a' + 10)
2409 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2410 unsigned int radix
, enum t_exactness
*p_exactness
)
2412 unsigned int idx
= *p_idx
;
2413 unsigned int hash_seen
= 0;
2414 scm_t_bits shift
= 1;
2416 unsigned int digit_value
;
2424 if (!isxdigit ((int) (unsigned char) c
))
2426 digit_value
= XDIGIT2UINT (c
);
2427 if (digit_value
>= radix
)
2431 result
= SCM_I_MAKINUM (digit_value
);
2435 if (isxdigit ((int) (unsigned char) c
))
2439 digit_value
= XDIGIT2UINT (c
);
2440 if (digit_value
>= radix
)
2452 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2454 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2456 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2463 shift
= shift
* radix
;
2464 add
= add
* radix
+ digit_value
;
2469 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2471 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2475 *p_exactness
= INEXACT
;
2481 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2482 * covers the parts of the rules that start at a potential point. The value
2483 * of the digits up to the point have been parsed by the caller and are given
2484 * in variable result. The content of *p_exactness indicates, whether a hash
2485 * has already been seen in the digits before the point.
2488 /* In non ASCII-style encodings the following macro might not work. */
2489 #define DIGIT2UINT(d) ((d) - '0')
2492 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2493 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2495 unsigned int idx
= *p_idx
;
2496 enum t_exactness x
= *p_exactness
;
2501 if (mem
[idx
] == '.')
2503 scm_t_bits shift
= 1;
2505 unsigned int digit_value
;
2506 SCM big_shift
= SCM_I_MAKINUM (1);
2512 if (isdigit ((int) (unsigned char) c
))
2517 digit_value
= DIGIT2UINT (c
);
2528 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2530 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2531 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2533 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2541 add
= add
* 10 + digit_value
;
2547 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2548 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2549 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2552 result
= scm_divide (result
, big_shift
);
2554 /* We've seen a decimal point, thus the value is implicitly inexact. */
2566 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2593 if (!isdigit ((int) (unsigned char) c
))
2597 exponent
= DIGIT2UINT (c
);
2601 if (isdigit ((int) (unsigned char) c
))
2604 if (exponent
<= SCM_MAXEXP
)
2605 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2611 if (exponent
> SCM_MAXEXP
)
2613 size_t exp_len
= idx
- start
;
2614 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2615 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2616 scm_out_of_range ("string->number", exp_num
);
2619 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2621 result
= scm_product (result
, e
);
2623 result
= scm_divide2real (result
, e
);
2625 /* We've seen an exponent, thus the value is implicitly inexact. */
2643 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2646 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2647 unsigned int radix
, enum t_exactness
*p_exactness
)
2649 unsigned int idx
= *p_idx
;
2655 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2661 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2663 enum t_exactness x
= EXACT
;
2665 /* Cobble up the fractional part. We might want to set the
2666 NaN's mantissa from it. */
2668 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2673 if (mem
[idx
] == '.')
2677 else if (idx
+ 1 == len
)
2679 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2682 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2683 p_idx
, p_exactness
);
2687 enum t_exactness x
= EXACT
;
2690 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2691 if (scm_is_false (uinteger
))
2696 else if (mem
[idx
] == '/')
2702 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2703 if (scm_is_false (divisor
))
2706 /* both are int/big here, I assume */
2707 result
= scm_i_make_ratio (uinteger
, divisor
);
2709 else if (radix
== 10)
2711 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2712 if (scm_is_false (result
))
2723 /* When returning an inexact zero, make sure it is represented as a
2724 floating point value so that we can change its sign.
2726 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2727 result
= scm_from_double (0.0);
2733 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2736 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2737 unsigned int radix
, enum t_exactness
*p_exactness
)
2761 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2762 if (scm_is_false (ureal
))
2764 /* input must be either +i or -i */
2769 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2775 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2782 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2783 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2792 /* either +<ureal>i or -<ureal>i */
2799 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2802 /* polar input: <real>@<real>. */
2827 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2828 if (scm_is_false (angle
))
2833 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2834 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2836 result
= scm_make_polar (ureal
, angle
);
2841 /* expecting input matching <real>[+-]<ureal>?i */
2848 int sign
= (c
== '+') ? 1 : -1;
2849 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2851 if (scm_is_false (imag
))
2852 imag
= SCM_I_MAKINUM (sign
);
2853 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2854 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2858 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2865 return scm_make_rectangular (ureal
, imag
);
2874 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2876 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2879 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2881 unsigned int idx
= 0;
2882 unsigned int radix
= NO_RADIX
;
2883 enum t_exactness forced_x
= NO_EXACTNESS
;
2884 enum t_exactness implicit_x
= EXACT
;
2887 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2888 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2890 switch (mem
[idx
+ 1])
2893 if (radix
!= NO_RADIX
)
2898 if (radix
!= NO_RADIX
)
2903 if (forced_x
!= NO_EXACTNESS
)
2908 if (forced_x
!= NO_EXACTNESS
)
2913 if (radix
!= NO_RADIX
)
2918 if (radix
!= NO_RADIX
)
2928 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2929 if (radix
== NO_RADIX
)
2930 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2932 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2934 if (scm_is_false (result
))
2940 if (SCM_INEXACTP (result
))
2941 return scm_inexact_to_exact (result
);
2945 if (SCM_INEXACTP (result
))
2948 return scm_exact_to_inexact (result
);
2951 if (implicit_x
== INEXACT
)
2953 if (SCM_INEXACTP (result
))
2956 return scm_exact_to_inexact (result
);
2964 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2965 (SCM string
, SCM radix
),
2966 "Return a number of the maximally precise representation\n"
2967 "expressed by the given @var{string}. @var{radix} must be an\n"
2968 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2969 "is a default radix that may be overridden by an explicit radix\n"
2970 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2971 "supplied, then the default radix is 10. If string is not a\n"
2972 "syntactically valid notation for a number, then\n"
2973 "@code{string->number} returns @code{#f}.")
2974 #define FUNC_NAME s_scm_string_to_number
2978 SCM_VALIDATE_STRING (1, string
);
2980 if (SCM_UNBNDP (radix
))
2983 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
2985 answer
= scm_i_mem2number (scm_i_string_chars (string
),
2986 scm_i_string_length (string
),
2988 scm_remember_upto_here_1 (string
);
2994 /*** END strs->nums ***/
2998 scm_bigequal (SCM x
, SCM y
)
3000 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3001 scm_remember_upto_here_2 (x
, y
);
3002 return scm_from_bool (0 == result
);
3006 scm_real_equalp (SCM x
, SCM y
)
3008 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3012 scm_complex_equalp (SCM x
, SCM y
)
3014 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3015 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3019 scm_i_fraction_equalp (SCM x
, SCM y
)
3021 scm_i_fraction_reduce (x
);
3022 scm_i_fraction_reduce (y
);
3023 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3024 SCM_FRACTION_NUMERATOR (y
)))
3025 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3026 SCM_FRACTION_DENOMINATOR (y
))))
3033 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3035 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3037 #define FUNC_NAME s_scm_number_p
3039 return scm_from_bool (SCM_NUMBERP (x
));
3043 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3045 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3046 "otherwise. Note that the sets of real, rational and integer\n"
3047 "values form subsets of the set of complex numbers, i. e. the\n"
3048 "predicate will also be fulfilled if @var{x} is a real,\n"
3049 "rational or integer number.")
3050 #define FUNC_NAME s_scm_complex_p
3052 /* all numbers are complex. */
3053 return scm_number_p (x
);
3057 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3059 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3060 "otherwise. Note that the set of integer values forms a subset of\n"
3061 "the set of real numbers, i. e. the predicate will also be\n"
3062 "fulfilled if @var{x} is an integer number.")
3063 #define FUNC_NAME s_scm_real_p
3065 /* we can't represent irrational numbers. */
3066 return scm_rational_p (x
);
3070 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3072 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3073 "otherwise. Note that the set of integer values forms a subset of\n"
3074 "the set of rational numbers, i. e. the predicate will also be\n"
3075 "fulfilled if @var{x} is an integer number.")
3076 #define FUNC_NAME s_scm_rational_p
3078 if (SCM_I_INUMP (x
))
3080 else if (SCM_IMP (x
))
3082 else if (SCM_BIGP (x
))
3084 else if (SCM_FRACTIONP (x
))
3086 else if (SCM_REALP (x
))
3087 /* due to their limited precision, all floating point numbers are
3088 rational as well. */
3095 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3097 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3099 #define FUNC_NAME s_scm_integer_p
3102 if (SCM_I_INUMP (x
))
3108 if (!SCM_INEXACTP (x
))
3110 if (SCM_COMPLEXP (x
))
3112 r
= SCM_REAL_VALUE (x
);
3113 /* +/-inf passes r==floor(r), making those #t */
3121 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3123 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3125 #define FUNC_NAME s_scm_inexact_p
3127 if (SCM_INEXACTP (x
))
3129 if (SCM_NUMBERP (x
))
3131 SCM_WRONG_TYPE_ARG (1, x
);
3136 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3137 /* "Return @code{#t} if all parameters are numerically equal." */
3139 scm_num_eq_p (SCM x
, SCM y
)
3142 if (SCM_I_INUMP (x
))
3144 long xx
= SCM_I_INUM (x
);
3145 if (SCM_I_INUMP (y
))
3147 long yy
= SCM_I_INUM (y
);
3148 return scm_from_bool (xx
== yy
);
3150 else if (SCM_BIGP (y
))
3152 else if (SCM_REALP (y
))
3153 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3154 else if (SCM_COMPLEXP (y
))
3155 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3156 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3157 else if (SCM_FRACTIONP (y
))
3160 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3162 else if (SCM_BIGP (x
))
3164 if (SCM_I_INUMP (y
))
3166 else if (SCM_BIGP (y
))
3168 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3169 scm_remember_upto_here_2 (x
, y
);
3170 return scm_from_bool (0 == cmp
);
3172 else if (SCM_REALP (y
))
3175 if (xisnan (SCM_REAL_VALUE (y
)))
3177 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3178 scm_remember_upto_here_1 (x
);
3179 return scm_from_bool (0 == cmp
);
3181 else if (SCM_COMPLEXP (y
))
3184 if (0.0 != SCM_COMPLEX_IMAG (y
))
3186 if (xisnan (SCM_COMPLEX_REAL (y
)))
3188 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3189 scm_remember_upto_here_1 (x
);
3190 return scm_from_bool (0 == cmp
);
3192 else if (SCM_FRACTIONP (y
))
3195 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3197 else if (SCM_REALP (x
))
3199 if (SCM_I_INUMP (y
))
3200 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3201 else if (SCM_BIGP (y
))
3204 if (xisnan (SCM_REAL_VALUE (x
)))
3206 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3207 scm_remember_upto_here_1 (y
);
3208 return scm_from_bool (0 == cmp
);
3210 else if (SCM_REALP (y
))
3211 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3212 else if (SCM_COMPLEXP (y
))
3213 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3214 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3215 else if (SCM_FRACTIONP (y
))
3217 double xx
= SCM_REAL_VALUE (x
);
3221 return scm_from_bool (xx
< 0.0);
3222 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3226 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3228 else if (SCM_COMPLEXP (x
))
3230 if (SCM_I_INUMP (y
))
3231 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3232 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3233 else if (SCM_BIGP (y
))
3236 if (0.0 != SCM_COMPLEX_IMAG (x
))
3238 if (xisnan (SCM_COMPLEX_REAL (x
)))
3240 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3241 scm_remember_upto_here_1 (y
);
3242 return scm_from_bool (0 == cmp
);
3244 else if (SCM_REALP (y
))
3245 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3246 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3247 else if (SCM_COMPLEXP (y
))
3248 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3249 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3250 else if (SCM_FRACTIONP (y
))
3253 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3255 xx
= SCM_COMPLEX_REAL (x
);
3259 return scm_from_bool (xx
< 0.0);
3260 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3264 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3266 else if (SCM_FRACTIONP (x
))
3268 if (SCM_I_INUMP (y
))
3270 else if (SCM_BIGP (y
))
3272 else if (SCM_REALP (y
))
3274 double yy
= SCM_REAL_VALUE (y
);
3278 return scm_from_bool (0.0 < yy
);
3279 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3282 else if (SCM_COMPLEXP (y
))
3285 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3287 yy
= SCM_COMPLEX_REAL (y
);
3291 return scm_from_bool (0.0 < yy
);
3292 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3295 else if (SCM_FRACTIONP (y
))
3296 return scm_i_fraction_equalp (x
, y
);
3298 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3301 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3305 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3306 done are good for inums, but for bignums an answer can almost always be
3307 had by just examining a few high bits of the operands, as done by GMP in
3308 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3309 of the float exponent to take into account. */
3311 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3312 /* "Return @code{#t} if the list of parameters is monotonically\n"
3316 scm_less_p (SCM x
, SCM y
)
3319 if (SCM_I_INUMP (x
))
3321 long xx
= SCM_I_INUM (x
);
3322 if (SCM_I_INUMP (y
))
3324 long yy
= SCM_I_INUM (y
);
3325 return scm_from_bool (xx
< yy
);
3327 else if (SCM_BIGP (y
))
3329 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3330 scm_remember_upto_here_1 (y
);
3331 return scm_from_bool (sgn
> 0);
3333 else if (SCM_REALP (y
))
3334 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3335 else if (SCM_FRACTIONP (y
))
3337 /* "x < a/b" becomes "x*b < a" */
3339 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3340 y
= SCM_FRACTION_NUMERATOR (y
);
3344 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3346 else if (SCM_BIGP (x
))
3348 if (SCM_I_INUMP (y
))
3350 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3351 scm_remember_upto_here_1 (x
);
3352 return scm_from_bool (sgn
< 0);
3354 else if (SCM_BIGP (y
))
3356 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3357 scm_remember_upto_here_2 (x
, y
);
3358 return scm_from_bool (cmp
< 0);
3360 else if (SCM_REALP (y
))
3363 if (xisnan (SCM_REAL_VALUE (y
)))
3365 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3366 scm_remember_upto_here_1 (x
);
3367 return scm_from_bool (cmp
< 0);
3369 else if (SCM_FRACTIONP (y
))
3372 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3374 else if (SCM_REALP (x
))
3376 if (SCM_I_INUMP (y
))
3377 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3378 else if (SCM_BIGP (y
))
3381 if (xisnan (SCM_REAL_VALUE (x
)))
3383 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3384 scm_remember_upto_here_1 (y
);
3385 return scm_from_bool (cmp
> 0);
3387 else if (SCM_REALP (y
))
3388 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3389 else if (SCM_FRACTIONP (y
))
3391 double xx
= SCM_REAL_VALUE (x
);
3395 return scm_from_bool (xx
< 0.0);
3396 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3400 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3402 else if (SCM_FRACTIONP (x
))
3404 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3406 /* "a/b < y" becomes "a < y*b" */
3407 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3408 x
= SCM_FRACTION_NUMERATOR (x
);
3411 else if (SCM_REALP (y
))
3413 double yy
= SCM_REAL_VALUE (y
);
3417 return scm_from_bool (0.0 < yy
);
3418 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3421 else if (SCM_FRACTIONP (y
))
3423 /* "a/b < c/d" becomes "a*d < c*b" */
3424 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3425 SCM_FRACTION_DENOMINATOR (y
));
3426 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3427 SCM_FRACTION_DENOMINATOR (x
));
3433 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3436 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3440 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3441 /* "Return @code{#t} if the list of parameters is monotonically\n"
3444 #define FUNC_NAME s_scm_gr_p
3446 scm_gr_p (SCM x
, SCM y
)
3448 if (!SCM_NUMBERP (x
))
3449 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3450 else if (!SCM_NUMBERP (y
))
3451 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3453 return scm_less_p (y
, x
);
3458 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3459 /* "Return @code{#t} if the list of parameters is monotonically\n"
3462 #define FUNC_NAME s_scm_leq_p
3464 scm_leq_p (SCM x
, SCM y
)
3466 if (!SCM_NUMBERP (x
))
3467 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3468 else if (!SCM_NUMBERP (y
))
3469 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3470 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3473 return scm_not (scm_less_p (y
, x
));
3478 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3479 /* "Return @code{#t} if the list of parameters is monotonically\n"
3482 #define FUNC_NAME s_scm_geq_p
3484 scm_geq_p (SCM x
, SCM y
)
3486 if (!SCM_NUMBERP (x
))
3487 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3488 else if (!SCM_NUMBERP (y
))
3489 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3490 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3493 return scm_not (scm_less_p (x
, y
));
3498 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3499 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3505 if (SCM_I_INUMP (z
))
3506 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3507 else if (SCM_BIGP (z
))
3509 else if (SCM_REALP (z
))
3510 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3511 else if (SCM_COMPLEXP (z
))
3512 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3513 && SCM_COMPLEX_IMAG (z
) == 0.0);
3514 else if (SCM_FRACTIONP (z
))
3517 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3521 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3522 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3526 scm_positive_p (SCM x
)
3528 if (SCM_I_INUMP (x
))
3529 return scm_from_bool (SCM_I_INUM (x
) > 0);
3530 else if (SCM_BIGP (x
))
3532 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3533 scm_remember_upto_here_1 (x
);
3534 return scm_from_bool (sgn
> 0);
3536 else if (SCM_REALP (x
))
3537 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3538 else if (SCM_FRACTIONP (x
))
3539 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3541 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3545 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3546 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3550 scm_negative_p (SCM x
)
3552 if (SCM_I_INUMP (x
))
3553 return scm_from_bool (SCM_I_INUM (x
) < 0);
3554 else if (SCM_BIGP (x
))
3556 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3557 scm_remember_upto_here_1 (x
);
3558 return scm_from_bool (sgn
< 0);
3560 else if (SCM_REALP (x
))
3561 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3562 else if (SCM_FRACTIONP (x
))
3563 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3565 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3569 /* scm_min and scm_max return an inexact when either argument is inexact, as
3570 required by r5rs. On that basis, for exact/inexact combinations the
3571 exact is converted to inexact to compare and possibly return. This is
3572 unlike scm_less_p above which takes some trouble to preserve all bits in
3573 its test, such trouble is not required for min and max. */
3575 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3576 /* "Return the maximum of all parameter values."
3579 scm_max (SCM x
, SCM y
)
3584 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3585 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3588 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3591 if (SCM_I_INUMP (x
))
3593 long xx
= SCM_I_INUM (x
);
3594 if (SCM_I_INUMP (y
))
3596 long yy
= SCM_I_INUM (y
);
3597 return (xx
< yy
) ? y
: x
;
3599 else if (SCM_BIGP (y
))
3601 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3602 scm_remember_upto_here_1 (y
);
3603 return (sgn
< 0) ? x
: y
;
3605 else if (SCM_REALP (y
))
3608 /* if y==NaN then ">" is false and we return NaN */
3609 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3611 else if (SCM_FRACTIONP (y
))
3614 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3617 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3619 else if (SCM_BIGP (x
))
3621 if (SCM_I_INUMP (y
))
3623 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3624 scm_remember_upto_here_1 (x
);
3625 return (sgn
< 0) ? y
: x
;
3627 else if (SCM_BIGP (y
))
3629 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3630 scm_remember_upto_here_2 (x
, y
);
3631 return (cmp
> 0) ? x
: y
;
3633 else if (SCM_REALP (y
))
3635 /* if y==NaN then xx>yy is false, so we return the NaN y */
3638 xx
= scm_i_big2dbl (x
);
3639 yy
= SCM_REAL_VALUE (y
);
3640 return (xx
> yy
? scm_from_double (xx
) : y
);
3642 else if (SCM_FRACTIONP (y
))
3647 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3649 else if (SCM_REALP (x
))
3651 if (SCM_I_INUMP (y
))
3653 double z
= SCM_I_INUM (y
);
3654 /* if x==NaN then "<" is false and we return NaN */
3655 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3657 else if (SCM_BIGP (y
))
3662 else if (SCM_REALP (y
))
3664 /* if x==NaN then our explicit check means we return NaN
3665 if y==NaN then ">" is false and we return NaN
3666 calling isnan is unavoidable, since it's the only way to know
3667 which of x or y causes any compares to be false */
3668 double xx
= SCM_REAL_VALUE (x
);
3669 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3671 else if (SCM_FRACTIONP (y
))
3673 double yy
= scm_i_fraction2double (y
);
3674 double xx
= SCM_REAL_VALUE (x
);
3675 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3678 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3680 else if (SCM_FRACTIONP (x
))
3682 if (SCM_I_INUMP (y
))
3686 else if (SCM_BIGP (y
))
3690 else if (SCM_REALP (y
))
3692 double xx
= scm_i_fraction2double (x
);
3693 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3695 else if (SCM_FRACTIONP (y
))
3700 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3703 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3707 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3708 /* "Return the minium of all parameter values."
3711 scm_min (SCM x
, SCM y
)
3716 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3717 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3720 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3723 if (SCM_I_INUMP (x
))
3725 long xx
= SCM_I_INUM (x
);
3726 if (SCM_I_INUMP (y
))
3728 long yy
= SCM_I_INUM (y
);
3729 return (xx
< yy
) ? x
: y
;
3731 else if (SCM_BIGP (y
))
3733 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3734 scm_remember_upto_here_1 (y
);
3735 return (sgn
< 0) ? y
: x
;
3737 else if (SCM_REALP (y
))
3740 /* if y==NaN then "<" is false and we return NaN */
3741 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3743 else if (SCM_FRACTIONP (y
))
3746 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3749 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3751 else if (SCM_BIGP (x
))
3753 if (SCM_I_INUMP (y
))
3755 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3756 scm_remember_upto_here_1 (x
);
3757 return (sgn
< 0) ? x
: y
;
3759 else if (SCM_BIGP (y
))
3761 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3762 scm_remember_upto_here_2 (x
, y
);
3763 return (cmp
> 0) ? y
: x
;
3765 else if (SCM_REALP (y
))
3767 /* if y==NaN then xx<yy is false, so we return the NaN y */
3770 xx
= scm_i_big2dbl (x
);
3771 yy
= SCM_REAL_VALUE (y
);
3772 return (xx
< yy
? scm_from_double (xx
) : y
);
3774 else if (SCM_FRACTIONP (y
))
3779 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3781 else if (SCM_REALP (x
))
3783 if (SCM_I_INUMP (y
))
3785 double z
= SCM_I_INUM (y
);
3786 /* if x==NaN then "<" is false and we return NaN */
3787 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3789 else if (SCM_BIGP (y
))
3794 else if (SCM_REALP (y
))
3796 /* if x==NaN then our explicit check means we return NaN
3797 if y==NaN then "<" is false and we return NaN
3798 calling isnan is unavoidable, since it's the only way to know
3799 which of x or y causes any compares to be false */
3800 double xx
= SCM_REAL_VALUE (x
);
3801 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3803 else if (SCM_FRACTIONP (y
))
3805 double yy
= scm_i_fraction2double (y
);
3806 double xx
= SCM_REAL_VALUE (x
);
3807 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3810 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3812 else if (SCM_FRACTIONP (x
))
3814 if (SCM_I_INUMP (y
))
3818 else if (SCM_BIGP (y
))
3822 else if (SCM_REALP (y
))
3824 double xx
= scm_i_fraction2double (x
);
3825 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3827 else if (SCM_FRACTIONP (y
))
3832 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3835 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3839 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3840 /* "Return the sum of all parameter values. Return 0 if called without\n"
3844 scm_sum (SCM x
, SCM y
)
3848 if (SCM_NUMBERP (x
)) return x
;
3849 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3850 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3853 if (SCM_I_INUMP (x
))
3855 if (SCM_I_INUMP (y
))
3857 long xx
= SCM_I_INUM (x
);
3858 long yy
= SCM_I_INUM (y
);
3859 long int z
= xx
+ yy
;
3860 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3862 else if (SCM_BIGP (y
))
3867 else if (SCM_REALP (y
))
3869 long int xx
= SCM_I_INUM (x
);
3870 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3872 else if (SCM_COMPLEXP (y
))
3874 long int xx
= SCM_I_INUM (x
);
3875 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3876 SCM_COMPLEX_IMAG (y
));
3878 else if (SCM_FRACTIONP (y
))
3879 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3880 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3881 SCM_FRACTION_DENOMINATOR (y
));
3883 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3884 } else if (SCM_BIGP (x
))
3886 if (SCM_I_INUMP (y
))
3891 inum
= SCM_I_INUM (y
);
3894 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3897 SCM result
= scm_i_mkbig ();
3898 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3899 scm_remember_upto_here_1 (x
);
3900 /* we know the result will have to be a bignum */
3903 return scm_i_normbig (result
);
3907 SCM result
= scm_i_mkbig ();
3908 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3909 scm_remember_upto_here_1 (x
);
3910 /* we know the result will have to be a bignum */
3913 return scm_i_normbig (result
);
3916 else if (SCM_BIGP (y
))
3918 SCM result
= scm_i_mkbig ();
3919 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3920 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3921 mpz_add (SCM_I_BIG_MPZ (result
),
3924 scm_remember_upto_here_2 (x
, y
);
3925 /* we know the result will have to be a bignum */
3928 return scm_i_normbig (result
);
3930 else if (SCM_REALP (y
))
3932 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3933 scm_remember_upto_here_1 (x
);
3934 return scm_from_double (result
);
3936 else if (SCM_COMPLEXP (y
))
3938 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3939 + SCM_COMPLEX_REAL (y
));
3940 scm_remember_upto_here_1 (x
);
3941 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
3943 else if (SCM_FRACTIONP (y
))
3944 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3945 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3946 SCM_FRACTION_DENOMINATOR (y
));
3948 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3950 else if (SCM_REALP (x
))
3952 if (SCM_I_INUMP (y
))
3953 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3954 else if (SCM_BIGP (y
))
3956 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3957 scm_remember_upto_here_1 (y
);
3958 return scm_from_double (result
);
3960 else if (SCM_REALP (y
))
3961 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3962 else if (SCM_COMPLEXP (y
))
3963 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3964 SCM_COMPLEX_IMAG (y
));
3965 else if (SCM_FRACTIONP (y
))
3966 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3968 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3970 else if (SCM_COMPLEXP (x
))
3972 if (SCM_I_INUMP (y
))
3973 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3974 SCM_COMPLEX_IMAG (x
));
3975 else if (SCM_BIGP (y
))
3977 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3978 + SCM_COMPLEX_REAL (x
));
3979 scm_remember_upto_here_1 (y
);
3980 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
3982 else if (SCM_REALP (y
))
3983 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3984 SCM_COMPLEX_IMAG (x
));
3985 else if (SCM_COMPLEXP (y
))
3986 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3987 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3988 else if (SCM_FRACTIONP (y
))
3989 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3990 SCM_COMPLEX_IMAG (x
));
3992 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3994 else if (SCM_FRACTIONP (x
))
3996 if (SCM_I_INUMP (y
))
3997 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3998 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3999 SCM_FRACTION_DENOMINATOR (x
));
4000 else if (SCM_BIGP (y
))
4001 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4002 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4003 SCM_FRACTION_DENOMINATOR (x
));
4004 else if (SCM_REALP (y
))
4005 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4006 else if (SCM_COMPLEXP (y
))
4007 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4008 SCM_COMPLEX_IMAG (y
));
4009 else if (SCM_FRACTIONP (y
))
4010 /* a/b + c/d = (ad + bc) / bd */
4011 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4012 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4013 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4015 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4018 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4022 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4023 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4024 * the sum of all but the first argument are subtracted from the first
4026 #define FUNC_NAME s_difference
4028 scm_difference (SCM x
, SCM y
)
4033 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4035 if (SCM_I_INUMP (x
))
4037 long xx
= -SCM_I_INUM (x
);
4038 if (SCM_FIXABLE (xx
))
4039 return SCM_I_MAKINUM (xx
);
4041 return scm_i_long2big (xx
);
4043 else if (SCM_BIGP (x
))
4044 /* FIXME: do we really need to normalize here? */
4045 return scm_i_normbig (scm_i_clonebig (x
, 0));
4046 else if (SCM_REALP (x
))
4047 return scm_from_double (-SCM_REAL_VALUE (x
));
4048 else if (SCM_COMPLEXP (x
))
4049 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4050 -SCM_COMPLEX_IMAG (x
));
4051 else if (SCM_FRACTIONP (x
))
4052 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4053 SCM_FRACTION_DENOMINATOR (x
));
4055 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4058 if (SCM_I_INUMP (x
))
4060 if (SCM_I_INUMP (y
))
4062 long int xx
= SCM_I_INUM (x
);
4063 long int yy
= SCM_I_INUM (y
);
4064 long int z
= xx
- yy
;
4065 if (SCM_FIXABLE (z
))
4066 return SCM_I_MAKINUM (z
);
4068 return scm_i_long2big (z
);
4070 else if (SCM_BIGP (y
))
4072 /* inum-x - big-y */
4073 long xx
= SCM_I_INUM (x
);
4076 return scm_i_clonebig (y
, 0);
4079 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4080 SCM result
= scm_i_mkbig ();
4083 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4086 /* x - y == -(y + -x) */
4087 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4088 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4090 scm_remember_upto_here_1 (y
);
4092 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4093 /* we know the result will have to be a bignum */
4096 return scm_i_normbig (result
);
4099 else if (SCM_REALP (y
))
4101 long int xx
= SCM_I_INUM (x
);
4102 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4104 else if (SCM_COMPLEXP (y
))
4106 long int xx
= SCM_I_INUM (x
);
4107 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4108 - SCM_COMPLEX_IMAG (y
));
4110 else if (SCM_FRACTIONP (y
))
4111 /* a - b/c = (ac - b) / c */
4112 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4113 SCM_FRACTION_NUMERATOR (y
)),
4114 SCM_FRACTION_DENOMINATOR (y
));
4116 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4118 else if (SCM_BIGP (x
))
4120 if (SCM_I_INUMP (y
))
4122 /* big-x - inum-y */
4123 long yy
= SCM_I_INUM (y
);
4124 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4126 scm_remember_upto_here_1 (x
);
4128 return (SCM_FIXABLE (-yy
) ?
4129 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4132 SCM result
= scm_i_mkbig ();
4135 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4137 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4138 scm_remember_upto_here_1 (x
);
4140 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4141 /* we know the result will have to be a bignum */
4144 return scm_i_normbig (result
);
4147 else if (SCM_BIGP (y
))
4149 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4150 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4151 SCM result
= scm_i_mkbig ();
4152 mpz_sub (SCM_I_BIG_MPZ (result
),
4155 scm_remember_upto_here_2 (x
, y
);
4156 /* we know the result will have to be a bignum */
4157 if ((sgn_x
== 1) && (sgn_y
== -1))
4159 if ((sgn_x
== -1) && (sgn_y
== 1))
4161 return scm_i_normbig (result
);
4163 else if (SCM_REALP (y
))
4165 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4166 scm_remember_upto_here_1 (x
);
4167 return scm_from_double (result
);
4169 else if (SCM_COMPLEXP (y
))
4171 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4172 - SCM_COMPLEX_REAL (y
));
4173 scm_remember_upto_here_1 (x
);
4174 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4176 else if (SCM_FRACTIONP (y
))
4177 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4178 SCM_FRACTION_NUMERATOR (y
)),
4179 SCM_FRACTION_DENOMINATOR (y
));
4180 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4182 else if (SCM_REALP (x
))
4184 if (SCM_I_INUMP (y
))
4185 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4186 else if (SCM_BIGP (y
))
4188 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4189 scm_remember_upto_here_1 (x
);
4190 return scm_from_double (result
);
4192 else if (SCM_REALP (y
))
4193 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4194 else if (SCM_COMPLEXP (y
))
4195 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4196 -SCM_COMPLEX_IMAG (y
));
4197 else if (SCM_FRACTIONP (y
))
4198 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4200 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4202 else if (SCM_COMPLEXP (x
))
4204 if (SCM_I_INUMP (y
))
4205 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4206 SCM_COMPLEX_IMAG (x
));
4207 else if (SCM_BIGP (y
))
4209 double real_part
= (SCM_COMPLEX_REAL (x
)
4210 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4211 scm_remember_upto_here_1 (x
);
4212 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4214 else if (SCM_REALP (y
))
4215 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4216 SCM_COMPLEX_IMAG (x
));
4217 else if (SCM_COMPLEXP (y
))
4218 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4219 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4220 else if (SCM_FRACTIONP (y
))
4221 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4222 SCM_COMPLEX_IMAG (x
));
4224 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4226 else if (SCM_FRACTIONP (x
))
4228 if (SCM_I_INUMP (y
))
4229 /* a/b - c = (a - cb) / b */
4230 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4231 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4232 SCM_FRACTION_DENOMINATOR (x
));
4233 else if (SCM_BIGP (y
))
4234 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4235 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4236 SCM_FRACTION_DENOMINATOR (x
));
4237 else if (SCM_REALP (y
))
4238 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4239 else if (SCM_COMPLEXP (y
))
4240 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4241 -SCM_COMPLEX_IMAG (y
));
4242 else if (SCM_FRACTIONP (y
))
4243 /* a/b - c/d = (ad - bc) / bd */
4244 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4245 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4246 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4248 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4251 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4256 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4257 /* "Return the product of all arguments. If called without arguments,\n"
4261 scm_product (SCM x
, SCM y
)
4266 return SCM_I_MAKINUM (1L);
4267 else if (SCM_NUMBERP (x
))
4270 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4273 if (SCM_I_INUMP (x
))
4278 xx
= SCM_I_INUM (x
);
4282 case 0: return x
; break;
4283 case 1: return y
; break;
4286 if (SCM_I_INUMP (y
))
4288 long yy
= SCM_I_INUM (y
);
4290 SCM k
= SCM_I_MAKINUM (kk
);
4291 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4295 SCM result
= scm_i_long2big (xx
);
4296 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4297 return scm_i_normbig (result
);
4300 else if (SCM_BIGP (y
))
4302 SCM result
= scm_i_mkbig ();
4303 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4304 scm_remember_upto_here_1 (y
);
4307 else if (SCM_REALP (y
))
4308 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4309 else if (SCM_COMPLEXP (y
))
4310 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4311 xx
* SCM_COMPLEX_IMAG (y
));
4312 else if (SCM_FRACTIONP (y
))
4313 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4314 SCM_FRACTION_DENOMINATOR (y
));
4316 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4318 else if (SCM_BIGP (x
))
4320 if (SCM_I_INUMP (y
))
4325 else if (SCM_BIGP (y
))
4327 SCM result
= scm_i_mkbig ();
4328 mpz_mul (SCM_I_BIG_MPZ (result
),
4331 scm_remember_upto_here_2 (x
, y
);
4334 else if (SCM_REALP (y
))
4336 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4337 scm_remember_upto_here_1 (x
);
4338 return scm_from_double (result
);
4340 else if (SCM_COMPLEXP (y
))
4342 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4343 scm_remember_upto_here_1 (x
);
4344 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4345 z
* SCM_COMPLEX_IMAG (y
));
4347 else if (SCM_FRACTIONP (y
))
4348 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4349 SCM_FRACTION_DENOMINATOR (y
));
4351 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4353 else if (SCM_REALP (x
))
4355 if (SCM_I_INUMP (y
))
4356 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4357 else if (SCM_BIGP (y
))
4359 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4360 scm_remember_upto_here_1 (y
);
4361 return scm_from_double (result
);
4363 else if (SCM_REALP (y
))
4364 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4365 else if (SCM_COMPLEXP (y
))
4366 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4367 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4368 else if (SCM_FRACTIONP (y
))
4369 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4371 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4373 else if (SCM_COMPLEXP (x
))
4375 if (SCM_I_INUMP (y
))
4376 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4377 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4378 else if (SCM_BIGP (y
))
4380 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4381 scm_remember_upto_here_1 (y
);
4382 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4383 z
* SCM_COMPLEX_IMAG (x
));
4385 else if (SCM_REALP (y
))
4386 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4387 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4388 else if (SCM_COMPLEXP (y
))
4390 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4391 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4392 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4393 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4395 else if (SCM_FRACTIONP (y
))
4397 double yy
= scm_i_fraction2double (y
);
4398 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4399 yy
* SCM_COMPLEX_IMAG (x
));
4402 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4404 else if (SCM_FRACTIONP (x
))
4406 if (SCM_I_INUMP (y
))
4407 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4408 SCM_FRACTION_DENOMINATOR (x
));
4409 else if (SCM_BIGP (y
))
4410 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4411 SCM_FRACTION_DENOMINATOR (x
));
4412 else if (SCM_REALP (y
))
4413 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4414 else if (SCM_COMPLEXP (y
))
4416 double xx
= scm_i_fraction2double (x
);
4417 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4418 xx
* SCM_COMPLEX_IMAG (y
));
4420 else if (SCM_FRACTIONP (y
))
4421 /* a/b * c/d = ac / bd */
4422 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4423 SCM_FRACTION_NUMERATOR (y
)),
4424 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4425 SCM_FRACTION_DENOMINATOR (y
)));
4427 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4430 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4433 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4434 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4435 #define ALLOW_DIVIDE_BY_ZERO
4436 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4439 /* The code below for complex division is adapted from the GNU
4440 libstdc++, which adapted it from f2c's libF77, and is subject to
4443 /****************************************************************
4444 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4446 Permission to use, copy, modify, and distribute this software
4447 and its documentation for any purpose and without fee is hereby
4448 granted, provided that the above copyright notice appear in all
4449 copies and that both that the copyright notice and this
4450 permission notice and warranty disclaimer appear in supporting
4451 documentation, and that the names of AT&T Bell Laboratories or
4452 Bellcore or any of their entities not be used in advertising or
4453 publicity pertaining to distribution of the software without
4454 specific, written prior permission.
4456 AT&T and Bellcore disclaim all warranties with regard to this
4457 software, including all implied warranties of merchantability
4458 and fitness. In no event shall AT&T or Bellcore be liable for
4459 any special, indirect or consequential damages or any damages
4460 whatsoever resulting from loss of use, data or profits, whether
4461 in an action of contract, negligence or other tortious action,
4462 arising out of or in connection with the use or performance of
4464 ****************************************************************/
4466 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4467 /* Divide the first argument by the product of the remaining
4468 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4470 #define FUNC_NAME s_divide
4472 scm_i_divide (SCM x
, SCM y
, int inexact
)
4479 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4480 else if (SCM_I_INUMP (x
))
4482 long xx
= SCM_I_INUM (x
);
4483 if (xx
== 1 || xx
== -1)
4485 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4487 scm_num_overflow (s_divide
);
4492 return scm_from_double (1.0 / (double) xx
);
4493 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4496 else if (SCM_BIGP (x
))
4499 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4500 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4502 else if (SCM_REALP (x
))
4504 double xx
= SCM_REAL_VALUE (x
);
4505 #ifndef ALLOW_DIVIDE_BY_ZERO
4507 scm_num_overflow (s_divide
);
4510 return scm_from_double (1.0 / xx
);
4512 else if (SCM_COMPLEXP (x
))
4514 double r
= SCM_COMPLEX_REAL (x
);
4515 double i
= SCM_COMPLEX_IMAG (x
);
4519 double d
= i
* (1.0 + t
* t
);
4520 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4525 double d
= r
* (1.0 + t
* t
);
4526 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4529 else if (SCM_FRACTIONP (x
))
4530 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4531 SCM_FRACTION_NUMERATOR (x
));
4533 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4536 if (SCM_I_INUMP (x
))
4538 long xx
= SCM_I_INUM (x
);
4539 if (SCM_I_INUMP (y
))
4541 long yy
= SCM_I_INUM (y
);
4544 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4545 scm_num_overflow (s_divide
);
4547 return scm_from_double ((double) xx
/ (double) yy
);
4550 else if (xx
% yy
!= 0)
4553 return scm_from_double ((double) xx
/ (double) yy
);
4554 else return scm_i_make_ratio (x
, y
);
4559 if (SCM_FIXABLE (z
))
4560 return SCM_I_MAKINUM (z
);
4562 return scm_i_long2big (z
);
4565 else if (SCM_BIGP (y
))
4568 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4569 else return scm_i_make_ratio (x
, y
);
4571 else if (SCM_REALP (y
))
4573 double yy
= SCM_REAL_VALUE (y
);
4574 #ifndef ALLOW_DIVIDE_BY_ZERO
4576 scm_num_overflow (s_divide
);
4579 return scm_from_double ((double) xx
/ yy
);
4581 else if (SCM_COMPLEXP (y
))
4584 complex_div
: /* y _must_ be a complex number */
4586 double r
= SCM_COMPLEX_REAL (y
);
4587 double i
= SCM_COMPLEX_IMAG (y
);
4591 double d
= i
* (1.0 + t
* t
);
4592 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4597 double d
= r
* (1.0 + t
* t
);
4598 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4602 else if (SCM_FRACTIONP (y
))
4603 /* a / b/c = ac / b */
4604 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4605 SCM_FRACTION_NUMERATOR (y
));
4607 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4609 else if (SCM_BIGP (x
))
4611 if (SCM_I_INUMP (y
))
4613 long int yy
= SCM_I_INUM (y
);
4616 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4617 scm_num_overflow (s_divide
);
4619 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4620 scm_remember_upto_here_1 (x
);
4621 return (sgn
== 0) ? scm_nan () : scm_inf ();
4628 /* FIXME: HMM, what are the relative performance issues here?
4629 We need to test. Is it faster on average to test
4630 divisible_p, then perform whichever operation, or is it
4631 faster to perform the integer div opportunistically and
4632 switch to real if there's a remainder? For now we take the
4633 middle ground: test, then if divisible, use the faster div
4636 long abs_yy
= yy
< 0 ? -yy
: yy
;
4637 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4641 SCM result
= scm_i_mkbig ();
4642 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4643 scm_remember_upto_here_1 (x
);
4645 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4646 return scm_i_normbig (result
);
4651 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4652 else return scm_i_make_ratio (x
, y
);
4656 else if (SCM_BIGP (y
))
4658 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4661 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4662 scm_num_overflow (s_divide
);
4664 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4665 scm_remember_upto_here_1 (x
);
4666 return (sgn
== 0) ? scm_nan () : scm_inf ();
4672 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4676 SCM result
= scm_i_mkbig ();
4677 mpz_divexact (SCM_I_BIG_MPZ (result
),
4680 scm_remember_upto_here_2 (x
, y
);
4681 return scm_i_normbig (result
);
4687 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4688 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4689 scm_remember_upto_here_2 (x
, y
);
4690 return scm_from_double (dbx
/ dby
);
4692 else return scm_i_make_ratio (x
, y
);
4696 else if (SCM_REALP (y
))
4698 double yy
= SCM_REAL_VALUE (y
);
4699 #ifndef ALLOW_DIVIDE_BY_ZERO
4701 scm_num_overflow (s_divide
);
4704 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4706 else if (SCM_COMPLEXP (y
))
4708 a
= scm_i_big2dbl (x
);
4711 else if (SCM_FRACTIONP (y
))
4712 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4713 SCM_FRACTION_NUMERATOR (y
));
4715 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4717 else if (SCM_REALP (x
))
4719 double rx
= SCM_REAL_VALUE (x
);
4720 if (SCM_I_INUMP (y
))
4722 long int yy
= SCM_I_INUM (y
);
4723 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4725 scm_num_overflow (s_divide
);
4728 return scm_from_double (rx
/ (double) yy
);
4730 else if (SCM_BIGP (y
))
4732 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4733 scm_remember_upto_here_1 (y
);
4734 return scm_from_double (rx
/ dby
);
4736 else if (SCM_REALP (y
))
4738 double yy
= SCM_REAL_VALUE (y
);
4739 #ifndef ALLOW_DIVIDE_BY_ZERO
4741 scm_num_overflow (s_divide
);
4744 return scm_from_double (rx
/ yy
);
4746 else if (SCM_COMPLEXP (y
))
4751 else if (SCM_FRACTIONP (y
))
4752 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4754 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4756 else if (SCM_COMPLEXP (x
))
4758 double rx
= SCM_COMPLEX_REAL (x
);
4759 double ix
= SCM_COMPLEX_IMAG (x
);
4760 if (SCM_I_INUMP (y
))
4762 long int yy
= SCM_I_INUM (y
);
4763 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4765 scm_num_overflow (s_divide
);
4770 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4773 else if (SCM_BIGP (y
))
4775 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4776 scm_remember_upto_here_1 (y
);
4777 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4779 else if (SCM_REALP (y
))
4781 double yy
= SCM_REAL_VALUE (y
);
4782 #ifndef ALLOW_DIVIDE_BY_ZERO
4784 scm_num_overflow (s_divide
);
4787 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4789 else if (SCM_COMPLEXP (y
))
4791 double ry
= SCM_COMPLEX_REAL (y
);
4792 double iy
= SCM_COMPLEX_IMAG (y
);
4796 double d
= iy
* (1.0 + t
* t
);
4797 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4802 double d
= ry
* (1.0 + t
* t
);
4803 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4806 else if (SCM_FRACTIONP (y
))
4808 double yy
= scm_i_fraction2double (y
);
4809 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4812 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4814 else if (SCM_FRACTIONP (x
))
4816 if (SCM_I_INUMP (y
))
4818 long int yy
= SCM_I_INUM (y
);
4819 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4821 scm_num_overflow (s_divide
);
4824 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4825 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4827 else if (SCM_BIGP (y
))
4829 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4830 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4832 else if (SCM_REALP (y
))
4834 double yy
= SCM_REAL_VALUE (y
);
4835 #ifndef ALLOW_DIVIDE_BY_ZERO
4837 scm_num_overflow (s_divide
);
4840 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4842 else if (SCM_COMPLEXP (y
))
4844 a
= scm_i_fraction2double (x
);
4847 else if (SCM_FRACTIONP (y
))
4848 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4849 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4851 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4854 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4858 scm_divide (SCM x
, SCM y
)
4860 return scm_i_divide (x
, y
, 0);
4863 static SCM
scm_divide2real (SCM x
, SCM y
)
4865 return scm_i_divide (x
, y
, 1);
4871 scm_asinh (double x
)
4876 #define asinh scm_asinh
4877 return log (x
+ sqrt (x
* x
+ 1));
4880 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4881 /* "Return the inverse hyperbolic sine of @var{x}."
4886 scm_acosh (double x
)
4891 #define acosh scm_acosh
4892 return log (x
+ sqrt (x
* x
- 1));
4895 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4896 /* "Return the inverse hyperbolic cosine of @var{x}."
4901 scm_atanh (double x
)
4906 #define atanh scm_atanh
4907 return 0.5 * log ((1 + x
) / (1 - x
));
4910 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4911 /* "Return the inverse hyperbolic tangent of @var{x}."
4916 scm_c_truncate (double x
)
4927 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4928 half-way case (ie. when x is an integer plus 0.5) going upwards.
4929 Then half-way cases are identified and adjusted down if the
4930 round-upwards didn't give the desired even integer.
4932 "plus_half == result" identifies a half-way case. If plus_half, which is
4933 x + 0.5, is an integer then x must be an integer plus 0.5.
4935 An odd "result" value is identified with result/2 != floor(result/2).
4936 This is done with plus_half, since that value is ready for use sooner in
4937 a pipelined cpu, and we're already requiring plus_half == result.
4939 Note however that we need to be careful when x is big and already an
4940 integer. In that case "x+0.5" may round to an adjacent integer, causing
4941 us to return such a value, incorrectly. For instance if the hardware is
4942 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4943 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4944 returned. Or if the hardware is in round-upwards mode, then other bigger
4945 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4946 representable value, 2^128+2^76 (or whatever), again incorrect.
4948 These bad roundings of x+0.5 are avoided by testing at the start whether
4949 x is already an integer. If it is then clearly that's the desired result
4950 already. And if it's not then the exponent must be small enough to allow
4951 an 0.5 to be represented, and hence added without a bad rounding. */
4954 scm_c_round (double x
)
4956 double plus_half
, result
;
4961 plus_half
= x
+ 0.5;
4962 result
= floor (plus_half
);
4963 /* Adjust so that the rounding is towards even. */
4964 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4969 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4971 "Round the number @var{x} towards zero.")
4972 #define FUNC_NAME s_scm_truncate_number
4974 if (scm_is_false (scm_negative_p (x
)))
4975 return scm_floor (x
);
4977 return scm_ceiling (x
);
4981 static SCM exactly_one_half
;
4983 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4985 "Round the number @var{x} towards the nearest integer. "
4986 "When it is exactly halfway between two integers, "
4987 "round towards the even one.")
4988 #define FUNC_NAME s_scm_round_number
4990 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4992 else if (SCM_REALP (x
))
4993 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
4996 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
4997 single quotient+remainder division then examining to see which way
4998 the rounding should go. */
4999 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5000 SCM result
= scm_floor (plus_half
);
5001 /* Adjust so that the rounding is towards even. */
5002 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5003 && scm_is_true (scm_odd_p (result
)))
5004 return scm_difference (result
, SCM_I_MAKINUM (1));
5011 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5013 "Round the number @var{x} towards minus infinity.")
5014 #define FUNC_NAME s_scm_floor
5016 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5018 else if (SCM_REALP (x
))
5019 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5020 else if (SCM_FRACTIONP (x
))
5022 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5023 SCM_FRACTION_DENOMINATOR (x
));
5024 if (scm_is_false (scm_negative_p (x
)))
5026 /* For positive x, rounding towards zero is correct. */
5031 /* For negative x, we need to return q-1 unless x is an
5032 integer. But fractions are never integer, per our
5034 return scm_difference (q
, SCM_I_MAKINUM (1));
5038 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5042 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5044 "Round the number @var{x} towards infinity.")
5045 #define FUNC_NAME s_scm_ceiling
5047 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5049 else if (SCM_REALP (x
))
5050 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5051 else if (SCM_FRACTIONP (x
))
5053 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5054 SCM_FRACTION_DENOMINATOR (x
));
5055 if (scm_is_false (scm_positive_p (x
)))
5057 /* For negative x, rounding towards zero is correct. */
5062 /* For positive x, we need to return q+1 unless x is an
5063 integer. But fractions are never integer, per our
5065 return scm_sum (q
, SCM_I_MAKINUM (1));
5069 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5073 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5074 /* "Return the square root of the real number @var{x}."
5076 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5077 /* "Return the absolute value of the real number @var{x}."
5079 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5080 /* "Return the @var{x}th power of e."
5082 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5083 /* "Return the natural logarithm of the real number @var{x}."
5085 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5086 /* "Return the sine of the real number @var{x}."
5088 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5089 /* "Return the cosine of the real number @var{x}."
5091 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5092 /* "Return the tangent of the real number @var{x}."
5094 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5095 /* "Return the arc sine of the real number @var{x}."
5097 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5098 /* "Return the arc cosine of the real number @var{x}."
5100 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5101 /* "Return the arc tangent of the real number @var{x}."
5103 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5104 /* "Return the hyperbolic sine of the real number @var{x}."
5106 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5107 /* "Return the hyperbolic cosine of the real number @var{x}."
5109 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5110 /* "Return the hyperbolic tangent of the real number @var{x}."
5118 static void scm_two_doubles (SCM x
,
5120 const char *sstring
,
5124 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5126 if (SCM_I_INUMP (x
))
5127 xy
->x
= SCM_I_INUM (x
);
5128 else if (SCM_BIGP (x
))
5129 xy
->x
= scm_i_big2dbl (x
);
5130 else if (SCM_REALP (x
))
5131 xy
->x
= SCM_REAL_VALUE (x
);
5132 else if (SCM_FRACTIONP (x
))
5133 xy
->x
= scm_i_fraction2double (x
);
5135 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5137 if (SCM_I_INUMP (y
))
5138 xy
->y
= SCM_I_INUM (y
);
5139 else if (SCM_BIGP (y
))
5140 xy
->y
= scm_i_big2dbl (y
);
5141 else if (SCM_REALP (y
))
5142 xy
->y
= SCM_REAL_VALUE (y
);
5143 else if (SCM_FRACTIONP (y
))
5144 xy
->y
= scm_i_fraction2double (y
);
5146 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5150 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5152 "Return @var{x} raised to the power of @var{y}. This\n"
5153 "procedure does not accept complex arguments.")
5154 #define FUNC_NAME s_scm_sys_expt
5157 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5158 return scm_from_double (pow (xy
.x
, xy
.y
));
5163 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5165 "Return the arc tangent of the two arguments @var{x} and\n"
5166 "@var{y}. This is similar to calculating the arc tangent of\n"
5167 "@var{x} / @var{y}, except that the signs of both arguments\n"
5168 "are used to determine the quadrant of the result. This\n"
5169 "procedure does not accept complex arguments.")
5170 #define FUNC_NAME s_scm_sys_atan2
5173 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5174 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5179 scm_c_make_rectangular (double re
, double im
)
5182 return scm_from_double (re
);
5186 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5188 SCM_COMPLEX_REAL (z
) = re
;
5189 SCM_COMPLEX_IMAG (z
) = im
;
5194 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5195 (SCM real
, SCM imaginary
),
5196 "Return a complex number constructed of the given @var{real} and\n"
5197 "@var{imaginary} parts.")
5198 #define FUNC_NAME s_scm_make_rectangular
5201 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5202 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5207 scm_c_make_polar (double mag
, double ang
)
5211 sincos (ang
, &s
, &c
);
5216 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5219 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5221 "Return the complex number @var{x} * e^(i * @var{y}).")
5222 #define FUNC_NAME s_scm_make_polar
5225 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5226 return scm_c_make_polar (xy
.x
, xy
.y
);
5231 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5232 /* "Return the real part of the number @var{z}."
5235 scm_real_part (SCM z
)
5237 if (SCM_I_INUMP (z
))
5239 else if (SCM_BIGP (z
))
5241 else if (SCM_REALP (z
))
5243 else if (SCM_COMPLEXP (z
))
5244 return scm_from_double (SCM_COMPLEX_REAL (z
));
5245 else if (SCM_FRACTIONP (z
))
5248 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5252 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5253 /* "Return the imaginary part of the number @var{z}."
5256 scm_imag_part (SCM z
)
5258 if (SCM_I_INUMP (z
))
5260 else if (SCM_BIGP (z
))
5262 else if (SCM_REALP (z
))
5264 else if (SCM_COMPLEXP (z
))
5265 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5266 else if (SCM_FRACTIONP (z
))
5269 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5272 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5273 /* "Return the numerator of the number @var{z}."
5276 scm_numerator (SCM z
)
5278 if (SCM_I_INUMP (z
))
5280 else if (SCM_BIGP (z
))
5282 else if (SCM_FRACTIONP (z
))
5284 scm_i_fraction_reduce (z
);
5285 return SCM_FRACTION_NUMERATOR (z
);
5287 else if (SCM_REALP (z
))
5288 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5290 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5294 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5295 /* "Return the denominator of the number @var{z}."
5298 scm_denominator (SCM z
)
5300 if (SCM_I_INUMP (z
))
5301 return SCM_I_MAKINUM (1);
5302 else if (SCM_BIGP (z
))
5303 return SCM_I_MAKINUM (1);
5304 else if (SCM_FRACTIONP (z
))
5306 scm_i_fraction_reduce (z
);
5307 return SCM_FRACTION_DENOMINATOR (z
);
5309 else if (SCM_REALP (z
))
5310 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5312 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5315 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5316 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5317 * "@code{abs} for real arguments, but also allows complex numbers."
5320 scm_magnitude (SCM z
)
5322 if (SCM_I_INUMP (z
))
5324 long int zz
= SCM_I_INUM (z
);
5327 else if (SCM_POSFIXABLE (-zz
))
5328 return SCM_I_MAKINUM (-zz
);
5330 return scm_i_long2big (-zz
);
5332 else if (SCM_BIGP (z
))
5334 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5335 scm_remember_upto_here_1 (z
);
5337 return scm_i_clonebig (z
, 0);
5341 else if (SCM_REALP (z
))
5342 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5343 else if (SCM_COMPLEXP (z
))
5344 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5345 else if (SCM_FRACTIONP (z
))
5347 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5349 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5350 SCM_FRACTION_DENOMINATOR (z
));
5353 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5357 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5358 /* "Return the angle of the complex number @var{z}."
5363 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5364 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5365 But if atan2 follows the floating point rounding mode, then the value
5366 is not a constant. Maybe it'd be close enough though. */
5367 if (SCM_I_INUMP (z
))
5369 if (SCM_I_INUM (z
) >= 0)
5372 return scm_from_double (atan2 (0.0, -1.0));
5374 else if (SCM_BIGP (z
))
5376 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5377 scm_remember_upto_here_1 (z
);
5379 return scm_from_double (atan2 (0.0, -1.0));
5383 else if (SCM_REALP (z
))
5385 if (SCM_REAL_VALUE (z
) >= 0)
5388 return scm_from_double (atan2 (0.0, -1.0));
5390 else if (SCM_COMPLEXP (z
))
5391 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5392 else if (SCM_FRACTIONP (z
))
5394 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5396 else return scm_from_double (atan2 (0.0, -1.0));
5399 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5403 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5404 /* Convert the number @var{x} to its inexact representation.\n"
5407 scm_exact_to_inexact (SCM z
)
5409 if (SCM_I_INUMP (z
))
5410 return scm_from_double ((double) SCM_I_INUM (z
));
5411 else if (SCM_BIGP (z
))
5412 return scm_from_double (scm_i_big2dbl (z
));
5413 else if (SCM_FRACTIONP (z
))
5414 return scm_from_double (scm_i_fraction2double (z
));
5415 else if (SCM_INEXACTP (z
))
5418 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5422 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5424 "Return an exact number that is numerically closest to @var{z}.")
5425 #define FUNC_NAME s_scm_inexact_to_exact
5427 if (SCM_I_INUMP (z
))
5429 else if (SCM_BIGP (z
))
5431 else if (SCM_REALP (z
))
5433 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5434 SCM_OUT_OF_RANGE (1, z
);
5441 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5442 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5443 scm_i_mpz2num (mpq_denref (frac
)));
5445 /* When scm_i_make_ratio throws, we leak the memory allocated
5452 else if (SCM_FRACTIONP (z
))
5455 SCM_WRONG_TYPE_ARG (1, z
);
5459 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5461 "Return an exact number that is within @var{err} of @var{x}.")
5462 #define FUNC_NAME s_scm_rationalize
5464 if (SCM_I_INUMP (x
))
5466 else if (SCM_BIGP (x
))
5468 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5470 /* Use continued fractions to find closest ratio. All
5471 arithmetic is done with exact numbers.
5474 SCM ex
= scm_inexact_to_exact (x
);
5475 SCM int_part
= scm_floor (ex
);
5476 SCM tt
= SCM_I_MAKINUM (1);
5477 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5478 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5482 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5485 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5486 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5488 /* We stop after a million iterations just to be absolutely sure
5489 that we don't go into an infinite loop. The process normally
5490 converges after less than a dozen iterations.
5493 err
= scm_abs (err
);
5494 while (++i
< 1000000)
5496 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5497 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5498 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5500 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5501 err
))) /* abs(x-a/b) <= err */
5503 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5504 if (scm_is_false (scm_exact_p (x
))
5505 || scm_is_false (scm_exact_p (err
)))
5506 return scm_exact_to_inexact (res
);
5510 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5512 tt
= scm_floor (rx
); /* tt = floor (rx) */
5518 scm_num_overflow (s_scm_rationalize
);
5521 SCM_WRONG_TYPE_ARG (1, x
);
5525 /* conversion functions */
5528 scm_is_integer (SCM val
)
5530 return scm_is_true (scm_integer_p (val
));
5534 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5536 if (SCM_I_INUMP (val
))
5538 scm_t_signed_bits n
= SCM_I_INUM (val
);
5539 return n
>= min
&& n
<= max
;
5541 else if (SCM_BIGP (val
))
5543 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5545 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5547 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5549 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5550 return n
>= min
&& n
<= max
;
5560 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5561 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5564 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5565 SCM_I_BIG_MPZ (val
));
5567 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5579 return n
>= min
&& n
<= max
;
5587 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5589 if (SCM_I_INUMP (val
))
5591 scm_t_signed_bits n
= SCM_I_INUM (val
);
5592 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5594 else if (SCM_BIGP (val
))
5596 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5598 else if (max
<= ULONG_MAX
)
5600 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5602 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5603 return n
>= min
&& n
<= max
;
5613 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5616 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5617 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5620 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5621 SCM_I_BIG_MPZ (val
));
5623 return n
>= min
&& n
<= max
;
5631 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5633 scm_error (scm_out_of_range_key
,
5635 "Value out of range ~S to ~S: ~S",
5636 scm_list_3 (min
, max
, bad_val
),
5637 scm_list_1 (bad_val
));
5640 #define TYPE scm_t_intmax
5641 #define TYPE_MIN min
5642 #define TYPE_MAX max
5643 #define SIZEOF_TYPE 0
5644 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5645 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5646 #include "libguile/conv-integer.i.c"
5648 #define TYPE scm_t_uintmax
5649 #define TYPE_MIN min
5650 #define TYPE_MAX max
5651 #define SIZEOF_TYPE 0
5652 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5653 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5654 #include "libguile/conv-uinteger.i.c"
5656 #define TYPE scm_t_int8
5657 #define TYPE_MIN SCM_T_INT8_MIN
5658 #define TYPE_MAX SCM_T_INT8_MAX
5659 #define SIZEOF_TYPE 1
5660 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5661 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5662 #include "libguile/conv-integer.i.c"
5664 #define TYPE scm_t_uint8
5666 #define TYPE_MAX SCM_T_UINT8_MAX
5667 #define SIZEOF_TYPE 1
5668 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5669 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5670 #include "libguile/conv-uinteger.i.c"
5672 #define TYPE scm_t_int16
5673 #define TYPE_MIN SCM_T_INT16_MIN
5674 #define TYPE_MAX SCM_T_INT16_MAX
5675 #define SIZEOF_TYPE 2
5676 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5677 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5678 #include "libguile/conv-integer.i.c"
5680 #define TYPE scm_t_uint16
5682 #define TYPE_MAX SCM_T_UINT16_MAX
5683 #define SIZEOF_TYPE 2
5684 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5685 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5686 #include "libguile/conv-uinteger.i.c"
5688 #define TYPE scm_t_int32
5689 #define TYPE_MIN SCM_T_INT32_MIN
5690 #define TYPE_MAX SCM_T_INT32_MAX
5691 #define SIZEOF_TYPE 4
5692 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5693 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5694 #include "libguile/conv-integer.i.c"
5696 #define TYPE scm_t_uint32
5698 #define TYPE_MAX SCM_T_UINT32_MAX
5699 #define SIZEOF_TYPE 4
5700 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5701 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5702 #include "libguile/conv-uinteger.i.c"
5704 #if SCM_HAVE_T_INT64
5706 #define TYPE scm_t_int64
5707 #define TYPE_MIN SCM_T_INT64_MIN
5708 #define TYPE_MAX SCM_T_INT64_MAX
5709 #define SIZEOF_TYPE 8
5710 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5711 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5712 #include "libguile/conv-integer.i.c"
5714 #define TYPE scm_t_uint64
5716 #define TYPE_MAX SCM_T_UINT64_MAX
5717 #define SIZEOF_TYPE 8
5718 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5719 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5720 #include "libguile/conv-uinteger.i.c"
5725 scm_to_mpz (SCM val
, mpz_t rop
)
5727 if (SCM_I_INUMP (val
))
5728 mpz_set_si (rop
, SCM_I_INUM (val
));
5729 else if (SCM_BIGP (val
))
5730 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5732 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5736 scm_from_mpz (mpz_t val
)
5738 return scm_i_mpz2num (val
);
5742 scm_is_real (SCM val
)
5744 return scm_is_true (scm_real_p (val
));
5748 scm_is_rational (SCM val
)
5750 return scm_is_true (scm_rational_p (val
));
5754 scm_to_double (SCM val
)
5756 if (SCM_I_INUMP (val
))
5757 return SCM_I_INUM (val
);
5758 else if (SCM_BIGP (val
))
5759 return scm_i_big2dbl (val
);
5760 else if (SCM_FRACTIONP (val
))
5761 return scm_i_fraction2double (val
);
5762 else if (SCM_REALP (val
))
5763 return SCM_REAL_VALUE (val
);
5765 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5769 scm_from_double (double val
)
5771 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5772 SCM_REAL_VALUE (z
) = val
;
5776 #if SCM_ENABLE_DISCOURAGED == 1
5779 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5783 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5787 scm_out_of_range (NULL
, num
);
5790 return scm_to_double (num
);
5794 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5798 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5802 scm_out_of_range (NULL
, num
);
5805 return scm_to_double (num
);
5811 scm_is_complex (SCM val
)
5813 return scm_is_true (scm_complex_p (val
));
5817 scm_c_real_part (SCM z
)
5819 if (SCM_COMPLEXP (z
))
5820 return SCM_COMPLEX_REAL (z
);
5823 /* Use the scm_real_part to get proper error checking and
5826 return scm_to_double (scm_real_part (z
));
5831 scm_c_imag_part (SCM z
)
5833 if (SCM_COMPLEXP (z
))
5834 return SCM_COMPLEX_IMAG (z
);
5837 /* Use the scm_imag_part to get proper error checking and
5838 dispatching. The result will almost always be 0.0, but not
5841 return scm_to_double (scm_imag_part (z
));
5846 scm_c_magnitude (SCM z
)
5848 return scm_to_double (scm_magnitude (z
));
5854 return scm_to_double (scm_angle (z
));
5858 scm_is_number (SCM z
)
5860 return scm_is_true (scm_number_p (z
));
5868 mpz_init_set_si (z_negative_one
, -1);
5870 /* It may be possible to tune the performance of some algorithms by using
5871 * the following constants to avoid the creation of bignums. Please, before
5872 * using these values, remember the two rules of program optimization:
5873 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5874 scm_c_define ("most-positive-fixnum",
5875 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5876 scm_c_define ("most-negative-fixnum",
5877 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5879 scm_add_feature ("complex");
5880 scm_add_feature ("inexact");
5881 scm_flo0
= scm_from_double (0.0);
5883 /* determine floating point precision */
5884 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5886 init_dblprec(&scm_dblprec
[i
-2],i
);
5887 init_fx_radix(fx_per_radix
[i
-2],i
);
5890 /* hard code precision for base 10 if the preprocessor tells us to... */
5891 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5894 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5895 SCM_I_MAKINUM (2)));
5896 #include "libguile/numbers.x"