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 /* FIXME: I think this may be a bug on some arches -- results
882 of % with negative second arg are undefined... */
900 return SCM_I_MAKINUM (result
);
903 else if (SCM_BIGP (y
))
905 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
912 SCM pos_y
= scm_i_clonebig (y
, 0);
913 /* do this after the last scm_op */
914 mpz_init_set_si (z_x
, xx
);
915 result
= pos_y
; /* re-use this bignum */
916 mpz_mod (SCM_I_BIG_MPZ (result
),
918 SCM_I_BIG_MPZ (pos_y
));
919 scm_remember_upto_here_1 (pos_y
);
923 result
= scm_i_mkbig ();
924 /* do this after the last scm_op */
925 mpz_init_set_si (z_x
, xx
);
926 mpz_mod (SCM_I_BIG_MPZ (result
),
929 scm_remember_upto_here_1 (y
);
932 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
933 mpz_add (SCM_I_BIG_MPZ (result
),
935 SCM_I_BIG_MPZ (result
));
936 scm_remember_upto_here_1 (y
);
937 /* and do this before the next one */
939 return scm_i_normbig (result
);
943 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
945 else if (SCM_BIGP (x
))
949 long yy
= SCM_I_INUM (y
);
951 scm_num_overflow (s_modulo
);
954 SCM result
= scm_i_mkbig ();
955 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
957 (yy
< 0) ? - yy
: yy
);
958 scm_remember_upto_here_1 (x
);
959 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
960 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
961 SCM_I_BIG_MPZ (result
),
963 return scm_i_normbig (result
);
966 else if (SCM_BIGP (y
))
969 SCM result
= scm_i_mkbig ();
970 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
971 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
972 mpz_mod (SCM_I_BIG_MPZ (result
),
974 SCM_I_BIG_MPZ (pos_y
));
976 scm_remember_upto_here_1 (x
);
977 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
978 mpz_add (SCM_I_BIG_MPZ (result
),
980 SCM_I_BIG_MPZ (result
));
981 scm_remember_upto_here_2 (y
, pos_y
);
982 return scm_i_normbig (result
);
986 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
989 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
992 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
993 /* "Return the greatest common divisor of all arguments.\n"
994 * "If called without arguments, 0 is returned."
997 scm_gcd (SCM x
, SCM y
)
1000 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1002 if (SCM_I_INUMP (x
))
1004 if (SCM_I_INUMP (y
))
1006 long xx
= SCM_I_INUM (x
);
1007 long yy
= SCM_I_INUM (y
);
1008 long u
= xx
< 0 ? -xx
: xx
;
1009 long v
= yy
< 0 ? -yy
: yy
;
1019 /* Determine a common factor 2^k */
1020 while (!(1 & (u
| v
)))
1026 /* Now, any factor 2^n can be eliminated */
1046 return (SCM_POSFIXABLE (result
)
1047 ? SCM_I_MAKINUM (result
)
1048 : scm_i_long2big (result
));
1050 else if (SCM_BIGP (y
))
1056 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1058 else if (SCM_BIGP (x
))
1060 if (SCM_I_INUMP (y
))
1062 unsigned long result
;
1065 yy
= SCM_I_INUM (y
);
1070 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1071 scm_remember_upto_here_1 (x
);
1072 return (SCM_POSFIXABLE (result
)
1073 ? SCM_I_MAKINUM (result
)
1074 : scm_from_ulong (result
));
1076 else if (SCM_BIGP (y
))
1078 SCM result
= scm_i_mkbig ();
1079 mpz_gcd (SCM_I_BIG_MPZ (result
),
1082 scm_remember_upto_here_2 (x
, y
);
1083 return scm_i_normbig (result
);
1086 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1089 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1092 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1093 /* "Return the least common multiple of the arguments.\n"
1094 * "If called without arguments, 1 is returned."
1097 scm_lcm (SCM n1
, SCM n2
)
1099 if (SCM_UNBNDP (n2
))
1101 if (SCM_UNBNDP (n1
))
1102 return SCM_I_MAKINUM (1L);
1103 n2
= SCM_I_MAKINUM (1L);
1106 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1107 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1108 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1109 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1111 if (SCM_I_INUMP (n1
))
1113 if (SCM_I_INUMP (n2
))
1115 SCM d
= scm_gcd (n1
, n2
);
1116 if (scm_is_eq (d
, SCM_INUM0
))
1119 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1123 /* inum n1, big n2 */
1126 SCM result
= scm_i_mkbig ();
1127 long nn1
= SCM_I_INUM (n1
);
1128 if (nn1
== 0) return SCM_INUM0
;
1129 if (nn1
< 0) nn1
= - nn1
;
1130 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1131 scm_remember_upto_here_1 (n2
);
1139 if (SCM_I_INUMP (n2
))
1146 SCM result
= scm_i_mkbig ();
1147 mpz_lcm(SCM_I_BIG_MPZ (result
),
1149 SCM_I_BIG_MPZ (n2
));
1150 scm_remember_upto_here_2(n1
, n2
);
1151 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1157 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1162 + + + x (map digit:logand X Y)
1163 + - + x (map digit:logand X (lognot (+ -1 Y)))
1164 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1165 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1170 + + + (map digit:logior X Y)
1171 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1172 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1173 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1178 + + + (map digit:logxor X Y)
1179 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1180 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1181 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1186 + + (any digit:logand X Y)
1187 + - (any digit:logand X (lognot (+ -1 Y)))
1188 - + (any digit:logand (lognot (+ -1 X)) Y)
1193 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1195 "Return the bitwise AND of the integer arguments.\n\n"
1197 "(logand) @result{} -1\n"
1198 "(logand 7) @result{} 7\n"
1199 "(logand #b111 #b011 #b001) @result{} 1\n"
1201 #define FUNC_NAME s_scm_logand
1205 if (SCM_UNBNDP (n2
))
1207 if (SCM_UNBNDP (n1
))
1208 return SCM_I_MAKINUM (-1);
1209 else if (!SCM_NUMBERP (n1
))
1210 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1211 else if (SCM_NUMBERP (n1
))
1214 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1217 if (SCM_I_INUMP (n1
))
1219 nn1
= SCM_I_INUM (n1
);
1220 if (SCM_I_INUMP (n2
))
1222 long nn2
= SCM_I_INUM (n2
);
1223 return SCM_I_MAKINUM (nn1
& nn2
);
1225 else if SCM_BIGP (n2
)
1231 SCM result_z
= scm_i_mkbig ();
1233 mpz_init_set_si (nn1_z
, nn1
);
1234 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1235 scm_remember_upto_here_1 (n2
);
1237 return scm_i_normbig (result_z
);
1241 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1243 else if (SCM_BIGP (n1
))
1245 if (SCM_I_INUMP (n2
))
1248 nn1
= SCM_I_INUM (n1
);
1251 else if (SCM_BIGP (n2
))
1253 SCM result_z
= scm_i_mkbig ();
1254 mpz_and (SCM_I_BIG_MPZ (result_z
),
1256 SCM_I_BIG_MPZ (n2
));
1257 scm_remember_upto_here_2 (n1
, n2
);
1258 return scm_i_normbig (result_z
);
1261 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1264 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1269 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1271 "Return the bitwise OR of the integer arguments.\n\n"
1273 "(logior) @result{} 0\n"
1274 "(logior 7) @result{} 7\n"
1275 "(logior #b000 #b001 #b011) @result{} 3\n"
1277 #define FUNC_NAME s_scm_logior
1281 if (SCM_UNBNDP (n2
))
1283 if (SCM_UNBNDP (n1
))
1285 else if (SCM_NUMBERP (n1
))
1288 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1291 if (SCM_I_INUMP (n1
))
1293 nn1
= SCM_I_INUM (n1
);
1294 if (SCM_I_INUMP (n2
))
1296 long nn2
= SCM_I_INUM (n2
);
1297 return SCM_I_MAKINUM (nn1
| nn2
);
1299 else if (SCM_BIGP (n2
))
1305 SCM result_z
= scm_i_mkbig ();
1307 mpz_init_set_si (nn1_z
, nn1
);
1308 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1309 scm_remember_upto_here_1 (n2
);
1315 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1317 else if (SCM_BIGP (n1
))
1319 if (SCM_I_INUMP (n2
))
1322 nn1
= SCM_I_INUM (n1
);
1325 else if (SCM_BIGP (n2
))
1327 SCM result_z
= scm_i_mkbig ();
1328 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1330 SCM_I_BIG_MPZ (n2
));
1331 scm_remember_upto_here_2 (n1
, n2
);
1335 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1338 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1343 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1345 "Return the bitwise XOR of the integer arguments. A bit is\n"
1346 "set in the result if it is set in an odd number of arguments.\n"
1348 "(logxor) @result{} 0\n"
1349 "(logxor 7) @result{} 7\n"
1350 "(logxor #b000 #b001 #b011) @result{} 2\n"
1351 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1353 #define FUNC_NAME s_scm_logxor
1357 if (SCM_UNBNDP (n2
))
1359 if (SCM_UNBNDP (n1
))
1361 else if (SCM_NUMBERP (n1
))
1364 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1367 if (SCM_I_INUMP (n1
))
1369 nn1
= SCM_I_INUM (n1
);
1370 if (SCM_I_INUMP (n2
))
1372 long nn2
= SCM_I_INUM (n2
);
1373 return SCM_I_MAKINUM (nn1
^ nn2
);
1375 else if (SCM_BIGP (n2
))
1379 SCM result_z
= scm_i_mkbig ();
1381 mpz_init_set_si (nn1_z
, nn1
);
1382 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1383 scm_remember_upto_here_1 (n2
);
1385 return scm_i_normbig (result_z
);
1389 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1391 else if (SCM_BIGP (n1
))
1393 if (SCM_I_INUMP (n2
))
1396 nn1
= SCM_I_INUM (n1
);
1399 else if (SCM_BIGP (n2
))
1401 SCM result_z
= scm_i_mkbig ();
1402 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1404 SCM_I_BIG_MPZ (n2
));
1405 scm_remember_upto_here_2 (n1
, n2
);
1406 return scm_i_normbig (result_z
);
1409 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1412 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1417 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1420 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1421 "(logtest #b0100 #b1011) @result{} #f\n"
1422 "(logtest #b0100 #b0111) @result{} #t\n"
1424 #define FUNC_NAME s_scm_logtest
1428 if (SCM_I_INUMP (j
))
1430 nj
= SCM_I_INUM (j
);
1431 if (SCM_I_INUMP (k
))
1433 long nk
= SCM_I_INUM (k
);
1434 return scm_from_bool (nj
& nk
);
1436 else if (SCM_BIGP (k
))
1444 mpz_init_set_si (nj_z
, nj
);
1445 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1446 scm_remember_upto_here_1 (k
);
1447 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1453 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1455 else if (SCM_BIGP (j
))
1457 if (SCM_I_INUMP (k
))
1460 nj
= SCM_I_INUM (j
);
1463 else if (SCM_BIGP (k
))
1467 mpz_init (result_z
);
1471 scm_remember_upto_here_2 (j
, k
);
1472 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1473 mpz_clear (result_z
);
1477 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1480 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1485 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1488 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1489 "(logbit? 0 #b1101) @result{} #t\n"
1490 "(logbit? 1 #b1101) @result{} #f\n"
1491 "(logbit? 2 #b1101) @result{} #t\n"
1492 "(logbit? 3 #b1101) @result{} #t\n"
1493 "(logbit? 4 #b1101) @result{} #f\n"
1495 #define FUNC_NAME s_scm_logbit_p
1497 unsigned long int iindex
;
1498 iindex
= scm_to_ulong (index
);
1500 if (SCM_I_INUMP (j
))
1502 /* bits above what's in an inum follow the sign bit */
1503 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1504 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1506 else if (SCM_BIGP (j
))
1508 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1509 scm_remember_upto_here_1 (j
);
1510 return scm_from_bool (val
);
1513 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1518 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1520 "Return the integer which is the ones-complement of the integer\n"
1524 "(number->string (lognot #b10000000) 2)\n"
1525 " @result{} \"-10000001\"\n"
1526 "(number->string (lognot #b0) 2)\n"
1527 " @result{} \"-1\"\n"
1529 #define FUNC_NAME s_scm_lognot
1531 if (SCM_I_INUMP (n
)) {
1532 /* No overflow here, just need to toggle all the bits making up the inum.
1533 Enhancement: No need to strip the tag and add it back, could just xor
1534 a block of 1 bits, if that worked with the various debug versions of
1536 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1538 } else if (SCM_BIGP (n
)) {
1539 SCM result
= scm_i_mkbig ();
1540 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1541 scm_remember_upto_here_1 (n
);
1545 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1550 /* returns 0 if IN is not an integer. OUT must already be
1553 coerce_to_big (SCM in
, mpz_t out
)
1556 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1557 else if (SCM_I_INUMP (in
))
1558 mpz_set_si (out
, SCM_I_INUM (in
));
1565 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1566 (SCM n
, SCM k
, SCM m
),
1567 "Return @var{n} raised to the integer exponent\n"
1568 "@var{k}, modulo @var{m}.\n"
1571 "(modulo-expt 2 3 5)\n"
1574 #define FUNC_NAME s_scm_modulo_expt
1580 /* There are two classes of error we might encounter --
1581 1) Math errors, which we'll report by calling scm_num_overflow,
1583 2) wrong-type errors, which of course we'll report by calling
1585 We don't report those errors immediately, however; instead we do
1586 some cleanup first. These variables tell us which error (if
1587 any) we should report after cleaning up.
1589 int report_overflow
= 0;
1591 int position_of_wrong_type
= 0;
1592 SCM value_of_wrong_type
= SCM_INUM0
;
1594 SCM result
= SCM_UNDEFINED
;
1600 if (scm_is_eq (m
, SCM_INUM0
))
1602 report_overflow
= 1;
1606 if (!coerce_to_big (n
, n_tmp
))
1608 value_of_wrong_type
= n
;
1609 position_of_wrong_type
= 1;
1613 if (!coerce_to_big (k
, k_tmp
))
1615 value_of_wrong_type
= k
;
1616 position_of_wrong_type
= 2;
1620 if (!coerce_to_big (m
, m_tmp
))
1622 value_of_wrong_type
= m
;
1623 position_of_wrong_type
= 3;
1627 /* if the exponent K is negative, and we simply call mpz_powm, we
1628 will get a divide-by-zero exception when an inverse 1/n mod m
1629 doesn't exist (or is not unique). Since exceptions are hard to
1630 handle, we'll attempt the inversion "by hand" -- that way, we get
1631 a simple failure code, which is easy to handle. */
1633 if (-1 == mpz_sgn (k_tmp
))
1635 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1637 report_overflow
= 1;
1640 mpz_neg (k_tmp
, k_tmp
);
1643 result
= scm_i_mkbig ();
1644 mpz_powm (SCM_I_BIG_MPZ (result
),
1649 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1650 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1657 if (report_overflow
)
1658 scm_num_overflow (FUNC_NAME
);
1660 if (position_of_wrong_type
)
1661 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1662 value_of_wrong_type
);
1664 return scm_i_normbig (result
);
1668 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1670 "Return @var{n} raised to the exact integer exponent\n"
1674 "(integer-expt 2 5)\n"
1676 "(integer-expt -3 3)\n"
1679 #define FUNC_NAME s_scm_integer_expt
1682 SCM z_i2
= SCM_BOOL_F
;
1684 SCM acc
= SCM_I_MAKINUM (1L);
1686 /* 0^0 == 1 according to R5RS */
1687 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1688 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1689 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1690 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1692 if (SCM_I_INUMP (k
))
1693 i2
= SCM_I_INUM (k
);
1694 else if (SCM_BIGP (k
))
1696 z_i2
= scm_i_clonebig (k
, 1);
1697 scm_remember_upto_here_1 (k
);
1701 SCM_WRONG_TYPE_ARG (2, k
);
1705 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1707 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1708 n
= scm_divide (n
, SCM_UNDEFINED
);
1712 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1716 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1718 return scm_product (acc
, n
);
1720 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1721 acc
= scm_product (acc
, n
);
1722 n
= scm_product (n
, n
);
1723 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1731 n
= scm_divide (n
, SCM_UNDEFINED
);
1738 return scm_product (acc
, n
);
1740 acc
= scm_product (acc
, n
);
1741 n
= scm_product (n
, n
);
1748 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1750 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1751 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1753 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1754 "@var{cnt} is negative it's a division, rounded towards negative\n"
1755 "infinity. (Note that this is not the same rounding as\n"
1756 "@code{quotient} does.)\n"
1758 "With @var{n} viewed as an infinite precision twos complement,\n"
1759 "@code{ash} means a left shift introducing zero bits, or a right\n"
1760 "shift dropping bits.\n"
1763 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1764 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1766 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1767 "(ash -23 -2) @result{} -6\n"
1769 #define FUNC_NAME s_scm_ash
1772 bits_to_shift
= scm_to_long (cnt
);
1774 if (bits_to_shift
< 0)
1776 /* Shift right by abs(cnt) bits. This is realized as a division
1777 by div:=2^abs(cnt). However, to guarantee the floor
1778 rounding, negative values require some special treatment.
1780 SCM div
= scm_integer_expt (SCM_I_MAKINUM (2),
1781 scm_from_long (-bits_to_shift
));
1783 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1784 if (scm_is_false (scm_negative_p (n
)))
1785 return scm_quotient (n
, div
);
1787 return scm_sum (SCM_I_MAKINUM (-1L),
1788 scm_quotient (scm_sum (SCM_I_MAKINUM (1L), n
), div
));
1791 /* Shift left is done by multiplication with 2^CNT */
1792 return scm_product (n
, scm_integer_expt (SCM_I_MAKINUM (2), cnt
));
1797 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1798 (SCM n
, SCM start
, SCM end
),
1799 "Return the integer composed of the @var{start} (inclusive)\n"
1800 "through @var{end} (exclusive) bits of @var{n}. The\n"
1801 "@var{start}th bit becomes the 0-th bit in the result.\n"
1804 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1805 " @result{} \"1010\"\n"
1806 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1807 " @result{} \"10110\"\n"
1809 #define FUNC_NAME s_scm_bit_extract
1811 unsigned long int istart
, iend
, bits
;
1812 istart
= scm_to_ulong (start
);
1813 iend
= scm_to_ulong (end
);
1814 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1816 /* how many bits to keep */
1817 bits
= iend
- istart
;
1819 if (SCM_I_INUMP (n
))
1821 long int in
= SCM_I_INUM (n
);
1823 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1824 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1825 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1827 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1829 /* Since we emulate two's complement encoded numbers, this
1830 * special case requires us to produce a result that has
1831 * more bits than can be stored in a fixnum.
1833 SCM result
= scm_i_long2big (in
);
1834 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1839 /* mask down to requisite bits */
1840 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1841 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1843 else if (SCM_BIGP (n
))
1848 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1852 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1853 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1854 such bits into a ulong. */
1855 result
= scm_i_mkbig ();
1856 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1857 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1858 result
= scm_i_normbig (result
);
1860 scm_remember_upto_here_1 (n
);
1864 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1869 static const char scm_logtab
[] = {
1870 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1873 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1875 "Return the number of bits in integer @var{n}. If integer is\n"
1876 "positive, the 1-bits in its binary representation are counted.\n"
1877 "If negative, the 0-bits in its two's-complement binary\n"
1878 "representation are counted. If 0, 0 is returned.\n"
1881 "(logcount #b10101010)\n"
1888 #define FUNC_NAME s_scm_logcount
1890 if (SCM_I_INUMP (n
))
1892 unsigned long int c
= 0;
1893 long int nn
= SCM_I_INUM (n
);
1898 c
+= scm_logtab
[15 & nn
];
1901 return SCM_I_MAKINUM (c
);
1903 else if (SCM_BIGP (n
))
1905 unsigned long count
;
1906 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1907 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1909 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1910 scm_remember_upto_here_1 (n
);
1911 return SCM_I_MAKINUM (count
);
1914 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1919 static const char scm_ilentab
[] = {
1920 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1924 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1926 "Return the number of bits necessary to represent @var{n}.\n"
1929 "(integer-length #b10101010)\n"
1931 "(integer-length 0)\n"
1933 "(integer-length #b1111)\n"
1936 #define FUNC_NAME s_scm_integer_length
1938 if (SCM_I_INUMP (n
))
1940 unsigned long int c
= 0;
1942 long int nn
= SCM_I_INUM (n
);
1948 l
= scm_ilentab
[15 & nn
];
1951 return SCM_I_MAKINUM (c
- 4 + l
);
1953 else if (SCM_BIGP (n
))
1955 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1956 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1957 1 too big, so check for that and adjust. */
1958 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1959 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1960 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1961 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1963 scm_remember_upto_here_1 (n
);
1964 return SCM_I_MAKINUM (size
);
1967 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1971 /*** NUMBERS -> STRINGS ***/
1972 #define SCM_MAX_DBL_PREC 60
1973 #define SCM_MAX_DBL_RADIX 36
1975 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1976 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1977 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1980 void init_dblprec(int *prec
, int radix
) {
1981 /* determine floating point precision by adding successively
1982 smaller increments to 1.0 until it is considered == 1.0 */
1983 double f
= ((double)1.0)/radix
;
1984 double fsum
= 1.0 + f
;
1989 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2001 void init_fx_radix(double *fx_list
, int radix
)
2003 /* initialize a per-radix list of tolerances. When added
2004 to a number < 1.0, we can determine if we should raund
2005 up and quit converting a number to a string. */
2009 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2010 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2013 /* use this array as a way to generate a single digit */
2014 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2017 idbl2str (double f
, char *a
, int radix
)
2019 int efmt
, dpt
, d
, i
, wp
;
2021 #ifdef DBL_MIN_10_EXP
2024 #endif /* DBL_MIN_10_EXP */
2029 radix
> SCM_MAX_DBL_RADIX
)
2031 /* revert to existing behavior */
2035 wp
= scm_dblprec
[radix
-2];
2036 fx
= fx_per_radix
[radix
-2];
2040 #ifdef HAVE_COPYSIGN
2041 double sgn
= copysign (1.0, f
);
2046 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2052 strcpy (a
, "-inf.0");
2054 strcpy (a
, "+inf.0");
2057 else if (xisnan (f
))
2059 strcpy (a
, "+nan.0");
2069 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2070 make-uniform-vector, from causing infinite loops. */
2071 /* just do the checking...if it passes, we do the conversion for our
2072 radix again below */
2079 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2087 while (f_cpy
> 10.0)
2090 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2111 if (f
+ fx
[wp
] >= radix
)
2118 /* adding 9999 makes this equivalent to abs(x) % 3 */
2119 dpt
= (exp
+ 9999) % 3;
2123 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2145 a
[ch
++] = number_chars
[d
];
2148 if (f
+ fx
[wp
] >= 1.0)
2150 a
[ch
- 1] = number_chars
[d
+1];
2162 if ((dpt
> 4) && (exp
> 6))
2164 d
= (a
[0] == '-' ? 2 : 1);
2165 for (i
= ch
++; i
> d
; i
--)
2178 if (a
[ch
- 1] == '.')
2179 a
[ch
++] = '0'; /* trailing zero */
2188 for (i
= radix
; i
<= exp
; i
*= radix
);
2189 for (i
/= radix
; i
; i
/= radix
)
2191 a
[ch
++] = number_chars
[exp
/ i
];
2199 iflo2str (SCM flt
, char *str
, int radix
)
2202 if (SCM_REALP (flt
))
2203 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2206 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
, radix
);
2207 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2209 double imag
= SCM_COMPLEX_IMAG (flt
);
2210 /* Don't output a '+' for negative numbers or for Inf and
2211 NaN. They will provide their own sign. */
2212 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2214 i
+= idbl2str (imag
, &str
[i
], radix
);
2221 /* convert a long to a string (unterminated). returns the number of
2222 characters in the result.
2224 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2226 scm_iint2str (long num
, int rad
, char *p
)
2230 unsigned long n
= (num
< 0) ? -num
: num
;
2232 for (n
/= rad
; n
> 0; n
/= rad
)
2249 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2254 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2256 "Return a string holding the external representation of the\n"
2257 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2258 "inexact, a radix of 10 will be used.")
2259 #define FUNC_NAME s_scm_number_to_string
2263 if (SCM_UNBNDP (radix
))
2266 base
= scm_to_signed_integer (radix
, 2, 36);
2268 if (SCM_I_INUMP (n
))
2270 char num_buf
[SCM_INTBUFLEN
];
2271 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2272 return scm_from_locale_stringn (num_buf
, length
);
2274 else if (SCM_BIGP (n
))
2276 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2277 scm_remember_upto_here_1 (n
);
2278 return scm_take_locale_string (str
);
2280 else if (SCM_FRACTIONP (n
))
2282 scm_i_fraction_reduce (n
);
2283 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2284 scm_from_locale_string ("/"),
2285 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2287 else if (SCM_INEXACTP (n
))
2289 char num_buf
[FLOBUFLEN
];
2290 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2293 SCM_WRONG_TYPE_ARG (1, n
);
2298 /* These print routines used to be stubbed here so that scm_repl.c
2299 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2302 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2304 char num_buf
[FLOBUFLEN
];
2305 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2310 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2313 char num_buf
[FLOBUFLEN
];
2314 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2319 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2322 scm_i_fraction_reduce (sexp
);
2323 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2324 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2325 scm_remember_upto_here_1 (str
);
2330 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2332 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2333 scm_remember_upto_here_1 (exp
);
2334 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2338 /*** END nums->strs ***/
2341 /*** STRINGS -> NUMBERS ***/
2343 /* The following functions implement the conversion from strings to numbers.
2344 * The implementation somehow follows the grammar for numbers as it is given
2345 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2346 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2347 * points should be noted about the implementation:
2348 * * Each function keeps a local index variable 'idx' that points at the
2349 * current position within the parsed string. The global index is only
2350 * updated if the function could parse the corresponding syntactic unit
2352 * * Similarly, the functions keep track of indicators of inexactness ('#',
2353 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2354 * global exactness information is only updated after each part has been
2355 * successfully parsed.
2356 * * Sequences of digits are parsed into temporary variables holding fixnums.
2357 * Only if these fixnums would overflow, the result variables are updated
2358 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2359 * the temporary variables holding the fixnums are cleared, and the process
2360 * starts over again. If for example fixnums were able to store five decimal
2361 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2362 * and the result was computed as 12345 * 100000 + 67890. In other words,
2363 * only every five digits two bignum operations were performed.
2366 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2368 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2370 /* In non ASCII-style encodings the following macro might not work. */
2371 #define XDIGIT2UINT(d) \
2372 (isdigit ((int) (unsigned char) d) \
2374 : tolower ((int) (unsigned char) d) - 'a' + 10)
2377 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2378 unsigned int radix
, enum t_exactness
*p_exactness
)
2380 unsigned int idx
= *p_idx
;
2381 unsigned int hash_seen
= 0;
2382 scm_t_bits shift
= 1;
2384 unsigned int digit_value
;
2392 if (!isxdigit ((int) (unsigned char) c
))
2394 digit_value
= XDIGIT2UINT (c
);
2395 if (digit_value
>= radix
)
2399 result
= SCM_I_MAKINUM (digit_value
);
2403 if (isxdigit ((int) (unsigned char) c
))
2407 digit_value
= XDIGIT2UINT (c
);
2408 if (digit_value
>= radix
)
2420 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2422 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2424 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2431 shift
= shift
* radix
;
2432 add
= add
* radix
+ digit_value
;
2437 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2439 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2443 *p_exactness
= INEXACT
;
2449 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2450 * covers the parts of the rules that start at a potential point. The value
2451 * of the digits up to the point have been parsed by the caller and are given
2452 * in variable result. The content of *p_exactness indicates, whether a hash
2453 * has already been seen in the digits before the point.
2456 /* In non ASCII-style encodings the following macro might not work. */
2457 #define DIGIT2UINT(d) ((d) - '0')
2460 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2461 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2463 unsigned int idx
= *p_idx
;
2464 enum t_exactness x
= *p_exactness
;
2469 if (mem
[idx
] == '.')
2471 scm_t_bits shift
= 1;
2473 unsigned int digit_value
;
2474 SCM big_shift
= SCM_I_MAKINUM (1);
2480 if (isdigit ((int) (unsigned char) c
))
2485 digit_value
= DIGIT2UINT (c
);
2496 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2498 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2499 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2501 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2509 add
= add
* 10 + digit_value
;
2515 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2516 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2517 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2520 result
= scm_divide (result
, big_shift
);
2522 /* We've seen a decimal point, thus the value is implicitly inexact. */
2534 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2561 if (!isdigit ((int) (unsigned char) c
))
2565 exponent
= DIGIT2UINT (c
);
2569 if (isdigit ((int) (unsigned char) c
))
2572 if (exponent
<= SCM_MAXEXP
)
2573 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2579 if (exponent
> SCM_MAXEXP
)
2581 size_t exp_len
= idx
- start
;
2582 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2583 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2584 scm_out_of_range ("string->number", exp_num
);
2587 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2589 result
= scm_product (result
, e
);
2591 result
= scm_divide2real (result
, e
);
2593 /* We've seen an exponent, thus the value is implicitly inexact. */
2611 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2614 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2615 unsigned int radix
, enum t_exactness
*p_exactness
)
2617 unsigned int idx
= *p_idx
;
2623 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2629 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2631 enum t_exactness x
= EXACT
;
2633 /* Cobble up the fractional part. We might want to set the
2634 NaN's mantissa from it. */
2636 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2641 if (mem
[idx
] == '.')
2645 else if (idx
+ 1 == len
)
2647 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2650 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2651 p_idx
, p_exactness
);
2655 enum t_exactness x
= EXACT
;
2658 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2659 if (scm_is_false (uinteger
))
2664 else if (mem
[idx
] == '/')
2670 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2671 if (scm_is_false (divisor
))
2674 /* both are int/big here, I assume */
2675 result
= scm_i_make_ratio (uinteger
, divisor
);
2677 else if (radix
== 10)
2679 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2680 if (scm_is_false (result
))
2691 /* When returning an inexact zero, make sure it is represented as a
2692 floating point value so that we can change its sign.
2694 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2695 result
= scm_from_double (0.0);
2701 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2704 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2705 unsigned int radix
, enum t_exactness
*p_exactness
)
2729 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2730 if (scm_is_false (ureal
))
2732 /* input must be either +i or -i */
2737 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2743 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2750 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2751 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2760 /* either +<ureal>i or -<ureal>i */
2767 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2770 /* polar input: <real>@<real>. */
2795 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2796 if (scm_is_false (angle
))
2801 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2802 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2804 result
= scm_make_polar (ureal
, angle
);
2809 /* expecting input matching <real>[+-]<ureal>?i */
2816 int sign
= (c
== '+') ? 1 : -1;
2817 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2819 if (scm_is_false (imag
))
2820 imag
= SCM_I_MAKINUM (sign
);
2821 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2822 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2826 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2833 return scm_make_rectangular (ureal
, imag
);
2842 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2844 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2847 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2849 unsigned int idx
= 0;
2850 unsigned int radix
= NO_RADIX
;
2851 enum t_exactness forced_x
= NO_EXACTNESS
;
2852 enum t_exactness implicit_x
= EXACT
;
2855 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2856 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2858 switch (mem
[idx
+ 1])
2861 if (radix
!= NO_RADIX
)
2866 if (radix
!= NO_RADIX
)
2871 if (forced_x
!= NO_EXACTNESS
)
2876 if (forced_x
!= NO_EXACTNESS
)
2881 if (radix
!= NO_RADIX
)
2886 if (radix
!= NO_RADIX
)
2896 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2897 if (radix
== NO_RADIX
)
2898 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2900 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2902 if (scm_is_false (result
))
2908 if (SCM_INEXACTP (result
))
2909 return scm_inexact_to_exact (result
);
2913 if (SCM_INEXACTP (result
))
2916 return scm_exact_to_inexact (result
);
2919 if (implicit_x
== INEXACT
)
2921 if (SCM_INEXACTP (result
))
2924 return scm_exact_to_inexact (result
);
2932 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2933 (SCM string
, SCM radix
),
2934 "Return a number of the maximally precise representation\n"
2935 "expressed by the given @var{string}. @var{radix} must be an\n"
2936 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2937 "is a default radix that may be overridden by an explicit radix\n"
2938 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2939 "supplied, then the default radix is 10. If string is not a\n"
2940 "syntactically valid notation for a number, then\n"
2941 "@code{string->number} returns @code{#f}.")
2942 #define FUNC_NAME s_scm_string_to_number
2946 SCM_VALIDATE_STRING (1, string
);
2948 if (SCM_UNBNDP (radix
))
2951 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
2953 answer
= scm_i_mem2number (scm_i_string_chars (string
),
2954 scm_i_string_length (string
),
2956 scm_remember_upto_here_1 (string
);
2962 /*** END strs->nums ***/
2966 scm_bigequal (SCM x
, SCM y
)
2968 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2969 scm_remember_upto_here_2 (x
, y
);
2970 return scm_from_bool (0 == result
);
2974 scm_real_equalp (SCM x
, SCM y
)
2976 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2980 scm_complex_equalp (SCM x
, SCM y
)
2982 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2983 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
2987 scm_i_fraction_equalp (SCM x
, SCM y
)
2989 scm_i_fraction_reduce (x
);
2990 scm_i_fraction_reduce (y
);
2991 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
2992 SCM_FRACTION_NUMERATOR (y
)))
2993 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
2994 SCM_FRACTION_DENOMINATOR (y
))))
3001 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3003 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3005 #define FUNC_NAME s_scm_number_p
3007 return scm_from_bool (SCM_NUMBERP (x
));
3011 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3013 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3014 "otherwise. Note that the sets of real, rational and integer\n"
3015 "values form subsets of the set of complex numbers, i. e. the\n"
3016 "predicate will also be fulfilled if @var{x} is a real,\n"
3017 "rational or integer number.")
3018 #define FUNC_NAME s_scm_complex_p
3020 /* all numbers are complex. */
3021 return scm_number_p (x
);
3025 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3027 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3028 "otherwise. Note that the set of integer values forms a subset of\n"
3029 "the set of real numbers, i. e. the predicate will also be\n"
3030 "fulfilled if @var{x} is an integer number.")
3031 #define FUNC_NAME s_scm_real_p
3033 /* we can't represent irrational numbers. */
3034 return scm_rational_p (x
);
3038 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3040 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3041 "otherwise. Note that the set of integer values forms a subset of\n"
3042 "the set of rational numbers, i. e. the predicate will also be\n"
3043 "fulfilled if @var{x} is an integer number.")
3044 #define FUNC_NAME s_scm_rational_p
3046 if (SCM_I_INUMP (x
))
3048 else if (SCM_IMP (x
))
3050 else if (SCM_BIGP (x
))
3052 else if (SCM_FRACTIONP (x
))
3054 else if (SCM_REALP (x
))
3055 /* due to their limited precision, all floating point numbers are
3056 rational as well. */
3063 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3065 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3067 #define FUNC_NAME s_scm_integer_p
3070 if (SCM_I_INUMP (x
))
3076 if (!SCM_INEXACTP (x
))
3078 if (SCM_COMPLEXP (x
))
3080 r
= SCM_REAL_VALUE (x
);
3081 /* +/-inf passes r==floor(r), making those #t */
3089 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3091 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3093 #define FUNC_NAME s_scm_inexact_p
3095 if (SCM_INEXACTP (x
))
3097 if (SCM_NUMBERP (x
))
3099 SCM_WRONG_TYPE_ARG (1, x
);
3104 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3105 /* "Return @code{#t} if all parameters are numerically equal." */
3107 scm_num_eq_p (SCM x
, SCM y
)
3110 if (SCM_I_INUMP (x
))
3112 long xx
= SCM_I_INUM (x
);
3113 if (SCM_I_INUMP (y
))
3115 long yy
= SCM_I_INUM (y
);
3116 return scm_from_bool (xx
== yy
);
3118 else if (SCM_BIGP (y
))
3120 else if (SCM_REALP (y
))
3121 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3122 else if (SCM_COMPLEXP (y
))
3123 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3124 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3125 else if (SCM_FRACTIONP (y
))
3128 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3130 else if (SCM_BIGP (x
))
3132 if (SCM_I_INUMP (y
))
3134 else if (SCM_BIGP (y
))
3136 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3137 scm_remember_upto_here_2 (x
, y
);
3138 return scm_from_bool (0 == cmp
);
3140 else if (SCM_REALP (y
))
3143 if (xisnan (SCM_REAL_VALUE (y
)))
3145 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3146 scm_remember_upto_here_1 (x
);
3147 return scm_from_bool (0 == cmp
);
3149 else if (SCM_COMPLEXP (y
))
3152 if (0.0 != SCM_COMPLEX_IMAG (y
))
3154 if (xisnan (SCM_COMPLEX_REAL (y
)))
3156 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3157 scm_remember_upto_here_1 (x
);
3158 return scm_from_bool (0 == cmp
);
3160 else if (SCM_FRACTIONP (y
))
3163 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3165 else if (SCM_REALP (x
))
3167 if (SCM_I_INUMP (y
))
3168 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3169 else if (SCM_BIGP (y
))
3172 if (xisnan (SCM_REAL_VALUE (x
)))
3174 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3175 scm_remember_upto_here_1 (y
);
3176 return scm_from_bool (0 == cmp
);
3178 else if (SCM_REALP (y
))
3179 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3180 else if (SCM_COMPLEXP (y
))
3181 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3182 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3183 else if (SCM_FRACTIONP (y
))
3185 double xx
= SCM_REAL_VALUE (x
);
3189 return scm_from_bool (xx
< 0.0);
3190 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3194 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3196 else if (SCM_COMPLEXP (x
))
3198 if (SCM_I_INUMP (y
))
3199 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3200 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3201 else if (SCM_BIGP (y
))
3204 if (0.0 != SCM_COMPLEX_IMAG (x
))
3206 if (xisnan (SCM_COMPLEX_REAL (x
)))
3208 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3209 scm_remember_upto_here_1 (y
);
3210 return scm_from_bool (0 == cmp
);
3212 else if (SCM_REALP (y
))
3213 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3214 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3215 else if (SCM_COMPLEXP (y
))
3216 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3217 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3218 else if (SCM_FRACTIONP (y
))
3221 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3223 xx
= SCM_COMPLEX_REAL (x
);
3227 return scm_from_bool (xx
< 0.0);
3228 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3232 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3234 else if (SCM_FRACTIONP (x
))
3236 if (SCM_I_INUMP (y
))
3238 else if (SCM_BIGP (y
))
3240 else if (SCM_REALP (y
))
3242 double yy
= SCM_REAL_VALUE (y
);
3246 return scm_from_bool (0.0 < yy
);
3247 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3250 else if (SCM_COMPLEXP (y
))
3253 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3255 yy
= SCM_COMPLEX_REAL (y
);
3259 return scm_from_bool (0.0 < yy
);
3260 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3263 else if (SCM_FRACTIONP (y
))
3264 return scm_i_fraction_equalp (x
, y
);
3266 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3269 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3273 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3274 done are good for inums, but for bignums an answer can almost always be
3275 had by just examining a few high bits of the operands, as done by GMP in
3276 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3277 of the float exponent to take into account. */
3279 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3280 /* "Return @code{#t} if the list of parameters is monotonically\n"
3284 scm_less_p (SCM x
, SCM y
)
3287 if (SCM_I_INUMP (x
))
3289 long xx
= SCM_I_INUM (x
);
3290 if (SCM_I_INUMP (y
))
3292 long yy
= SCM_I_INUM (y
);
3293 return scm_from_bool (xx
< yy
);
3295 else if (SCM_BIGP (y
))
3297 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3298 scm_remember_upto_here_1 (y
);
3299 return scm_from_bool (sgn
> 0);
3301 else if (SCM_REALP (y
))
3302 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3303 else if (SCM_FRACTIONP (y
))
3305 /* "x < a/b" becomes "x*b < a" */
3307 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3308 y
= SCM_FRACTION_NUMERATOR (y
);
3312 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3314 else if (SCM_BIGP (x
))
3316 if (SCM_I_INUMP (y
))
3318 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3319 scm_remember_upto_here_1 (x
);
3320 return scm_from_bool (sgn
< 0);
3322 else if (SCM_BIGP (y
))
3324 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3325 scm_remember_upto_here_2 (x
, y
);
3326 return scm_from_bool (cmp
< 0);
3328 else if (SCM_REALP (y
))
3331 if (xisnan (SCM_REAL_VALUE (y
)))
3333 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3334 scm_remember_upto_here_1 (x
);
3335 return scm_from_bool (cmp
< 0);
3337 else if (SCM_FRACTIONP (y
))
3340 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3342 else if (SCM_REALP (x
))
3344 if (SCM_I_INUMP (y
))
3345 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3346 else if (SCM_BIGP (y
))
3349 if (xisnan (SCM_REAL_VALUE (x
)))
3351 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3352 scm_remember_upto_here_1 (y
);
3353 return scm_from_bool (cmp
> 0);
3355 else if (SCM_REALP (y
))
3356 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3357 else if (SCM_FRACTIONP (y
))
3359 double xx
= SCM_REAL_VALUE (x
);
3363 return scm_from_bool (xx
< 0.0);
3364 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3368 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3370 else if (SCM_FRACTIONP (x
))
3372 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3374 /* "a/b < y" becomes "a < y*b" */
3375 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3376 x
= SCM_FRACTION_NUMERATOR (x
);
3379 else if (SCM_REALP (y
))
3381 double yy
= SCM_REAL_VALUE (y
);
3385 return scm_from_bool (0.0 < yy
);
3386 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3389 else if (SCM_FRACTIONP (y
))
3391 /* "a/b < c/d" becomes "a*d < c*b" */
3392 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3393 SCM_FRACTION_DENOMINATOR (y
));
3394 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3395 SCM_FRACTION_DENOMINATOR (x
));
3401 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3404 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3408 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3409 /* "Return @code{#t} if the list of parameters is monotonically\n"
3412 #define FUNC_NAME s_scm_gr_p
3414 scm_gr_p (SCM x
, SCM y
)
3416 if (!SCM_NUMBERP (x
))
3417 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3418 else if (!SCM_NUMBERP (y
))
3419 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3421 return scm_less_p (y
, x
);
3426 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3427 /* "Return @code{#t} if the list of parameters is monotonically\n"
3430 #define FUNC_NAME s_scm_leq_p
3432 scm_leq_p (SCM x
, SCM y
)
3434 if (!SCM_NUMBERP (x
))
3435 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3436 else if (!SCM_NUMBERP (y
))
3437 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3438 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3441 return scm_not (scm_less_p (y
, x
));
3446 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3447 /* "Return @code{#t} if the list of parameters is monotonically\n"
3450 #define FUNC_NAME s_scm_geq_p
3452 scm_geq_p (SCM x
, SCM y
)
3454 if (!SCM_NUMBERP (x
))
3455 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3456 else if (!SCM_NUMBERP (y
))
3457 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3458 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3461 return scm_not (scm_less_p (x
, y
));
3466 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3467 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3473 if (SCM_I_INUMP (z
))
3474 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3475 else if (SCM_BIGP (z
))
3477 else if (SCM_REALP (z
))
3478 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3479 else if (SCM_COMPLEXP (z
))
3480 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3481 && SCM_COMPLEX_IMAG (z
) == 0.0);
3482 else if (SCM_FRACTIONP (z
))
3485 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3489 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3490 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3494 scm_positive_p (SCM x
)
3496 if (SCM_I_INUMP (x
))
3497 return scm_from_bool (SCM_I_INUM (x
) > 0);
3498 else if (SCM_BIGP (x
))
3500 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3501 scm_remember_upto_here_1 (x
);
3502 return scm_from_bool (sgn
> 0);
3504 else if (SCM_REALP (x
))
3505 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3506 else if (SCM_FRACTIONP (x
))
3507 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3509 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3513 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3514 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3518 scm_negative_p (SCM x
)
3520 if (SCM_I_INUMP (x
))
3521 return scm_from_bool (SCM_I_INUM (x
) < 0);
3522 else if (SCM_BIGP (x
))
3524 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3525 scm_remember_upto_here_1 (x
);
3526 return scm_from_bool (sgn
< 0);
3528 else if (SCM_REALP (x
))
3529 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3530 else if (SCM_FRACTIONP (x
))
3531 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3533 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3537 /* scm_min and scm_max return an inexact when either argument is inexact, as
3538 required by r5rs. On that basis, for exact/inexact combinations the
3539 exact is converted to inexact to compare and possibly return. This is
3540 unlike scm_less_p above which takes some trouble to preserve all bits in
3541 its test, such trouble is not required for min and max. */
3543 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3544 /* "Return the maximum of all parameter values."
3547 scm_max (SCM x
, SCM y
)
3552 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3553 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3556 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3559 if (SCM_I_INUMP (x
))
3561 long xx
= SCM_I_INUM (x
);
3562 if (SCM_I_INUMP (y
))
3564 long yy
= SCM_I_INUM (y
);
3565 return (xx
< yy
) ? y
: x
;
3567 else if (SCM_BIGP (y
))
3569 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3570 scm_remember_upto_here_1 (y
);
3571 return (sgn
< 0) ? x
: y
;
3573 else if (SCM_REALP (y
))
3576 /* if y==NaN then ">" is false and we return NaN */
3577 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3579 else if (SCM_FRACTIONP (y
))
3582 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3585 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3587 else if (SCM_BIGP (x
))
3589 if (SCM_I_INUMP (y
))
3591 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3592 scm_remember_upto_here_1 (x
);
3593 return (sgn
< 0) ? y
: x
;
3595 else if (SCM_BIGP (y
))
3597 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3598 scm_remember_upto_here_2 (x
, y
);
3599 return (cmp
> 0) ? x
: y
;
3601 else if (SCM_REALP (y
))
3603 /* if y==NaN then xx>yy is false, so we return the NaN y */
3606 xx
= scm_i_big2dbl (x
);
3607 yy
= SCM_REAL_VALUE (y
);
3608 return (xx
> yy
? scm_from_double (xx
) : y
);
3610 else if (SCM_FRACTIONP (y
))
3615 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3617 else if (SCM_REALP (x
))
3619 if (SCM_I_INUMP (y
))
3621 double z
= SCM_I_INUM (y
);
3622 /* if x==NaN then "<" is false and we return NaN */
3623 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3625 else if (SCM_BIGP (y
))
3630 else if (SCM_REALP (y
))
3632 /* if x==NaN then our explicit check means we return NaN
3633 if y==NaN then ">" is false and we return NaN
3634 calling isnan is unavoidable, since it's the only way to know
3635 which of x or y causes any compares to be false */
3636 double xx
= SCM_REAL_VALUE (x
);
3637 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3639 else if (SCM_FRACTIONP (y
))
3641 double yy
= scm_i_fraction2double (y
);
3642 double xx
= SCM_REAL_VALUE (x
);
3643 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3646 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3648 else if (SCM_FRACTIONP (x
))
3650 if (SCM_I_INUMP (y
))
3654 else if (SCM_BIGP (y
))
3658 else if (SCM_REALP (y
))
3660 double xx
= scm_i_fraction2double (x
);
3661 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3663 else if (SCM_FRACTIONP (y
))
3668 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3671 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3675 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3676 /* "Return the minium of all parameter values."
3679 scm_min (SCM x
, SCM y
)
3684 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3685 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3688 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3691 if (SCM_I_INUMP (x
))
3693 long xx
= SCM_I_INUM (x
);
3694 if (SCM_I_INUMP (y
))
3696 long yy
= SCM_I_INUM (y
);
3697 return (xx
< yy
) ? x
: y
;
3699 else if (SCM_BIGP (y
))
3701 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3702 scm_remember_upto_here_1 (y
);
3703 return (sgn
< 0) ? y
: x
;
3705 else if (SCM_REALP (y
))
3708 /* if y==NaN then "<" is false and we return NaN */
3709 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3711 else if (SCM_FRACTIONP (y
))
3714 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3717 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3719 else if (SCM_BIGP (x
))
3721 if (SCM_I_INUMP (y
))
3723 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3724 scm_remember_upto_here_1 (x
);
3725 return (sgn
< 0) ? x
: y
;
3727 else if (SCM_BIGP (y
))
3729 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3730 scm_remember_upto_here_2 (x
, y
);
3731 return (cmp
> 0) ? y
: x
;
3733 else if (SCM_REALP (y
))
3735 /* if y==NaN then xx<yy is false, so we return the NaN y */
3738 xx
= scm_i_big2dbl (x
);
3739 yy
= SCM_REAL_VALUE (y
);
3740 return (xx
< yy
? scm_from_double (xx
) : y
);
3742 else if (SCM_FRACTIONP (y
))
3747 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3749 else if (SCM_REALP (x
))
3751 if (SCM_I_INUMP (y
))
3753 double z
= SCM_I_INUM (y
);
3754 /* if x==NaN then "<" is false and we return NaN */
3755 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3757 else if (SCM_BIGP (y
))
3762 else if (SCM_REALP (y
))
3764 /* if x==NaN then our explicit check means we return NaN
3765 if y==NaN then "<" is false and we return NaN
3766 calling isnan is unavoidable, since it's the only way to know
3767 which of x or y causes any compares to be false */
3768 double xx
= SCM_REAL_VALUE (x
);
3769 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3771 else if (SCM_FRACTIONP (y
))
3773 double yy
= scm_i_fraction2double (y
);
3774 double xx
= SCM_REAL_VALUE (x
);
3775 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3778 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3780 else if (SCM_FRACTIONP (x
))
3782 if (SCM_I_INUMP (y
))
3786 else if (SCM_BIGP (y
))
3790 else if (SCM_REALP (y
))
3792 double xx
= scm_i_fraction2double (x
);
3793 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3795 else if (SCM_FRACTIONP (y
))
3800 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3803 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3807 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3808 /* "Return the sum of all parameter values. Return 0 if called without\n"
3812 scm_sum (SCM x
, SCM y
)
3816 if (SCM_NUMBERP (x
)) return x
;
3817 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3818 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3821 if (SCM_I_INUMP (x
))
3823 if (SCM_I_INUMP (y
))
3825 long xx
= SCM_I_INUM (x
);
3826 long yy
= SCM_I_INUM (y
);
3827 long int z
= xx
+ yy
;
3828 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3830 else if (SCM_BIGP (y
))
3835 else if (SCM_REALP (y
))
3837 long int xx
= SCM_I_INUM (x
);
3838 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3840 else if (SCM_COMPLEXP (y
))
3842 long int xx
= SCM_I_INUM (x
);
3843 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3844 SCM_COMPLEX_IMAG (y
));
3846 else if (SCM_FRACTIONP (y
))
3847 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3848 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3849 SCM_FRACTION_DENOMINATOR (y
));
3851 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3852 } else if (SCM_BIGP (x
))
3854 if (SCM_I_INUMP (y
))
3859 inum
= SCM_I_INUM (y
);
3862 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3865 SCM result
= scm_i_mkbig ();
3866 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3867 scm_remember_upto_here_1 (x
);
3868 /* we know the result will have to be a bignum */
3871 return scm_i_normbig (result
);
3875 SCM result
= scm_i_mkbig ();
3876 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3877 scm_remember_upto_here_1 (x
);
3878 /* we know the result will have to be a bignum */
3881 return scm_i_normbig (result
);
3884 else if (SCM_BIGP (y
))
3886 SCM result
= scm_i_mkbig ();
3887 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3888 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3889 mpz_add (SCM_I_BIG_MPZ (result
),
3892 scm_remember_upto_here_2 (x
, y
);
3893 /* we know the result will have to be a bignum */
3896 return scm_i_normbig (result
);
3898 else if (SCM_REALP (y
))
3900 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3901 scm_remember_upto_here_1 (x
);
3902 return scm_from_double (result
);
3904 else if (SCM_COMPLEXP (y
))
3906 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3907 + SCM_COMPLEX_REAL (y
));
3908 scm_remember_upto_here_1 (x
);
3909 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
3911 else if (SCM_FRACTIONP (y
))
3912 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3913 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3914 SCM_FRACTION_DENOMINATOR (y
));
3916 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3918 else if (SCM_REALP (x
))
3920 if (SCM_I_INUMP (y
))
3921 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3922 else if (SCM_BIGP (y
))
3924 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3925 scm_remember_upto_here_1 (y
);
3926 return scm_from_double (result
);
3928 else if (SCM_REALP (y
))
3929 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3930 else if (SCM_COMPLEXP (y
))
3931 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3932 SCM_COMPLEX_IMAG (y
));
3933 else if (SCM_FRACTIONP (y
))
3934 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3936 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3938 else if (SCM_COMPLEXP (x
))
3940 if (SCM_I_INUMP (y
))
3941 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3942 SCM_COMPLEX_IMAG (x
));
3943 else if (SCM_BIGP (y
))
3945 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3946 + SCM_COMPLEX_REAL (x
));
3947 scm_remember_upto_here_1 (y
);
3948 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
3950 else if (SCM_REALP (y
))
3951 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3952 SCM_COMPLEX_IMAG (x
));
3953 else if (SCM_COMPLEXP (y
))
3954 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3955 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3956 else if (SCM_FRACTIONP (y
))
3957 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3958 SCM_COMPLEX_IMAG (x
));
3960 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3962 else if (SCM_FRACTIONP (x
))
3964 if (SCM_I_INUMP (y
))
3965 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3966 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3967 SCM_FRACTION_DENOMINATOR (x
));
3968 else if (SCM_BIGP (y
))
3969 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3970 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3971 SCM_FRACTION_DENOMINATOR (x
));
3972 else if (SCM_REALP (y
))
3973 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3974 else if (SCM_COMPLEXP (y
))
3975 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3976 SCM_COMPLEX_IMAG (y
));
3977 else if (SCM_FRACTIONP (y
))
3978 /* a/b + c/d = (ad + bc) / bd */
3979 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3980 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3981 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3983 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3986 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
3990 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
3991 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
3992 * the sum of all but the first argument are subtracted from the first
3994 #define FUNC_NAME s_difference
3996 scm_difference (SCM x
, SCM y
)
4001 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4003 if (SCM_I_INUMP (x
))
4005 long xx
= -SCM_I_INUM (x
);
4006 if (SCM_FIXABLE (xx
))
4007 return SCM_I_MAKINUM (xx
);
4009 return scm_i_long2big (xx
);
4011 else if (SCM_BIGP (x
))
4012 /* FIXME: do we really need to normalize here? */
4013 return scm_i_normbig (scm_i_clonebig (x
, 0));
4014 else if (SCM_REALP (x
))
4015 return scm_from_double (-SCM_REAL_VALUE (x
));
4016 else if (SCM_COMPLEXP (x
))
4017 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4018 -SCM_COMPLEX_IMAG (x
));
4019 else if (SCM_FRACTIONP (x
))
4020 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4021 SCM_FRACTION_DENOMINATOR (x
));
4023 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4026 if (SCM_I_INUMP (x
))
4028 if (SCM_I_INUMP (y
))
4030 long int xx
= SCM_I_INUM (x
);
4031 long int yy
= SCM_I_INUM (y
);
4032 long int z
= xx
- yy
;
4033 if (SCM_FIXABLE (z
))
4034 return SCM_I_MAKINUM (z
);
4036 return scm_i_long2big (z
);
4038 else if (SCM_BIGP (y
))
4040 /* inum-x - big-y */
4041 long xx
= SCM_I_INUM (x
);
4044 return scm_i_clonebig (y
, 0);
4047 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4048 SCM result
= scm_i_mkbig ();
4051 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4054 /* x - y == -(y + -x) */
4055 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4056 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4058 scm_remember_upto_here_1 (y
);
4060 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4061 /* we know the result will have to be a bignum */
4064 return scm_i_normbig (result
);
4067 else if (SCM_REALP (y
))
4069 long int xx
= SCM_I_INUM (x
);
4070 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4072 else if (SCM_COMPLEXP (y
))
4074 long int xx
= SCM_I_INUM (x
);
4075 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4076 - SCM_COMPLEX_IMAG (y
));
4078 else if (SCM_FRACTIONP (y
))
4079 /* a - b/c = (ac - b) / c */
4080 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4081 SCM_FRACTION_NUMERATOR (y
)),
4082 SCM_FRACTION_DENOMINATOR (y
));
4084 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4086 else if (SCM_BIGP (x
))
4088 if (SCM_I_INUMP (y
))
4090 /* big-x - inum-y */
4091 long yy
= SCM_I_INUM (y
);
4092 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4094 scm_remember_upto_here_1 (x
);
4096 return (SCM_FIXABLE (-yy
) ?
4097 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4100 SCM result
= scm_i_mkbig ();
4103 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4105 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4106 scm_remember_upto_here_1 (x
);
4108 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4109 /* we know the result will have to be a bignum */
4112 return scm_i_normbig (result
);
4115 else if (SCM_BIGP (y
))
4117 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4118 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4119 SCM result
= scm_i_mkbig ();
4120 mpz_sub (SCM_I_BIG_MPZ (result
),
4123 scm_remember_upto_here_2 (x
, y
);
4124 /* we know the result will have to be a bignum */
4125 if ((sgn_x
== 1) && (sgn_y
== -1))
4127 if ((sgn_x
== -1) && (sgn_y
== 1))
4129 return scm_i_normbig (result
);
4131 else if (SCM_REALP (y
))
4133 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4134 scm_remember_upto_here_1 (x
);
4135 return scm_from_double (result
);
4137 else if (SCM_COMPLEXP (y
))
4139 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4140 - SCM_COMPLEX_REAL (y
));
4141 scm_remember_upto_here_1 (x
);
4142 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4144 else if (SCM_FRACTIONP (y
))
4145 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4146 SCM_FRACTION_NUMERATOR (y
)),
4147 SCM_FRACTION_DENOMINATOR (y
));
4148 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4150 else if (SCM_REALP (x
))
4152 if (SCM_I_INUMP (y
))
4153 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4154 else if (SCM_BIGP (y
))
4156 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4157 scm_remember_upto_here_1 (x
);
4158 return scm_from_double (result
);
4160 else if (SCM_REALP (y
))
4161 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4162 else if (SCM_COMPLEXP (y
))
4163 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4164 -SCM_COMPLEX_IMAG (y
));
4165 else if (SCM_FRACTIONP (y
))
4166 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4168 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4170 else if (SCM_COMPLEXP (x
))
4172 if (SCM_I_INUMP (y
))
4173 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4174 SCM_COMPLEX_IMAG (x
));
4175 else if (SCM_BIGP (y
))
4177 double real_part
= (SCM_COMPLEX_REAL (x
)
4178 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4179 scm_remember_upto_here_1 (x
);
4180 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4182 else if (SCM_REALP (y
))
4183 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4184 SCM_COMPLEX_IMAG (x
));
4185 else if (SCM_COMPLEXP (y
))
4186 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4187 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4188 else if (SCM_FRACTIONP (y
))
4189 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4190 SCM_COMPLEX_IMAG (x
));
4192 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4194 else if (SCM_FRACTIONP (x
))
4196 if (SCM_I_INUMP (y
))
4197 /* a/b - c = (a - cb) / b */
4198 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4199 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4200 SCM_FRACTION_DENOMINATOR (x
));
4201 else if (SCM_BIGP (y
))
4202 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4203 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4204 SCM_FRACTION_DENOMINATOR (x
));
4205 else if (SCM_REALP (y
))
4206 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4207 else if (SCM_COMPLEXP (y
))
4208 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4209 -SCM_COMPLEX_IMAG (y
));
4210 else if (SCM_FRACTIONP (y
))
4211 /* a/b - c/d = (ad - bc) / bd */
4212 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4213 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4214 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4216 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4219 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4224 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4225 /* "Return the product of all arguments. If called without arguments,\n"
4229 scm_product (SCM x
, SCM y
)
4234 return SCM_I_MAKINUM (1L);
4235 else if (SCM_NUMBERP (x
))
4238 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4241 if (SCM_I_INUMP (x
))
4246 xx
= SCM_I_INUM (x
);
4250 case 0: return x
; break;
4251 case 1: return y
; break;
4254 if (SCM_I_INUMP (y
))
4256 long yy
= SCM_I_INUM (y
);
4258 SCM k
= SCM_I_MAKINUM (kk
);
4259 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4263 SCM result
= scm_i_long2big (xx
);
4264 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4265 return scm_i_normbig (result
);
4268 else if (SCM_BIGP (y
))
4270 SCM result
= scm_i_mkbig ();
4271 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4272 scm_remember_upto_here_1 (y
);
4275 else if (SCM_REALP (y
))
4276 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4277 else if (SCM_COMPLEXP (y
))
4278 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4279 xx
* SCM_COMPLEX_IMAG (y
));
4280 else if (SCM_FRACTIONP (y
))
4281 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4282 SCM_FRACTION_DENOMINATOR (y
));
4284 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4286 else if (SCM_BIGP (x
))
4288 if (SCM_I_INUMP (y
))
4293 else if (SCM_BIGP (y
))
4295 SCM result
= scm_i_mkbig ();
4296 mpz_mul (SCM_I_BIG_MPZ (result
),
4299 scm_remember_upto_here_2 (x
, y
);
4302 else if (SCM_REALP (y
))
4304 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4305 scm_remember_upto_here_1 (x
);
4306 return scm_from_double (result
);
4308 else if (SCM_COMPLEXP (y
))
4310 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4311 scm_remember_upto_here_1 (x
);
4312 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4313 z
* SCM_COMPLEX_IMAG (y
));
4315 else if (SCM_FRACTIONP (y
))
4316 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4317 SCM_FRACTION_DENOMINATOR (y
));
4319 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4321 else if (SCM_REALP (x
))
4323 if (SCM_I_INUMP (y
))
4324 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4325 else if (SCM_BIGP (y
))
4327 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4328 scm_remember_upto_here_1 (y
);
4329 return scm_from_double (result
);
4331 else if (SCM_REALP (y
))
4332 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4333 else if (SCM_COMPLEXP (y
))
4334 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4335 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4336 else if (SCM_FRACTIONP (y
))
4337 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4339 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4341 else if (SCM_COMPLEXP (x
))
4343 if (SCM_I_INUMP (y
))
4344 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4345 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4346 else if (SCM_BIGP (y
))
4348 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4349 scm_remember_upto_here_1 (y
);
4350 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4351 z
* SCM_COMPLEX_IMAG (x
));
4353 else if (SCM_REALP (y
))
4354 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4355 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4356 else if (SCM_COMPLEXP (y
))
4358 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4359 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4360 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4361 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4363 else if (SCM_FRACTIONP (y
))
4365 double yy
= scm_i_fraction2double (y
);
4366 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4367 yy
* SCM_COMPLEX_IMAG (x
));
4370 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4372 else if (SCM_FRACTIONP (x
))
4374 if (SCM_I_INUMP (y
))
4375 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4376 SCM_FRACTION_DENOMINATOR (x
));
4377 else if (SCM_BIGP (y
))
4378 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4379 SCM_FRACTION_DENOMINATOR (x
));
4380 else if (SCM_REALP (y
))
4381 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4382 else if (SCM_COMPLEXP (y
))
4384 double xx
= scm_i_fraction2double (x
);
4385 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4386 xx
* SCM_COMPLEX_IMAG (y
));
4388 else if (SCM_FRACTIONP (y
))
4389 /* a/b * c/d = ac / bd */
4390 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4391 SCM_FRACTION_NUMERATOR (y
)),
4392 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4393 SCM_FRACTION_DENOMINATOR (y
)));
4395 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4398 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4401 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4402 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4403 #define ALLOW_DIVIDE_BY_ZERO
4404 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4407 /* The code below for complex division is adapted from the GNU
4408 libstdc++, which adapted it from f2c's libF77, and is subject to
4411 /****************************************************************
4412 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4414 Permission to use, copy, modify, and distribute this software
4415 and its documentation for any purpose and without fee is hereby
4416 granted, provided that the above copyright notice appear in all
4417 copies and that both that the copyright notice and this
4418 permission notice and warranty disclaimer appear in supporting
4419 documentation, and that the names of AT&T Bell Laboratories or
4420 Bellcore or any of their entities not be used in advertising or
4421 publicity pertaining to distribution of the software without
4422 specific, written prior permission.
4424 AT&T and Bellcore disclaim all warranties with regard to this
4425 software, including all implied warranties of merchantability
4426 and fitness. In no event shall AT&T or Bellcore be liable for
4427 any special, indirect or consequential damages or any damages
4428 whatsoever resulting from loss of use, data or profits, whether
4429 in an action of contract, negligence or other tortious action,
4430 arising out of or in connection with the use or performance of
4432 ****************************************************************/
4434 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4435 /* Divide the first argument by the product of the remaining
4436 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4438 #define FUNC_NAME s_divide
4440 scm_i_divide (SCM x
, SCM y
, int inexact
)
4447 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4448 else if (SCM_I_INUMP (x
))
4450 long xx
= SCM_I_INUM (x
);
4451 if (xx
== 1 || xx
== -1)
4453 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4455 scm_num_overflow (s_divide
);
4460 return scm_from_double (1.0 / (double) xx
);
4461 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4464 else if (SCM_BIGP (x
))
4467 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4468 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4470 else if (SCM_REALP (x
))
4472 double xx
= SCM_REAL_VALUE (x
);
4473 #ifndef ALLOW_DIVIDE_BY_ZERO
4475 scm_num_overflow (s_divide
);
4478 return scm_from_double (1.0 / xx
);
4480 else if (SCM_COMPLEXP (x
))
4482 double r
= SCM_COMPLEX_REAL (x
);
4483 double i
= SCM_COMPLEX_IMAG (x
);
4487 double d
= i
* (1.0 + t
* t
);
4488 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4493 double d
= r
* (1.0 + t
* t
);
4494 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4497 else if (SCM_FRACTIONP (x
))
4498 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4499 SCM_FRACTION_NUMERATOR (x
));
4501 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4504 if (SCM_I_INUMP (x
))
4506 long xx
= SCM_I_INUM (x
);
4507 if (SCM_I_INUMP (y
))
4509 long yy
= SCM_I_INUM (y
);
4512 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4513 scm_num_overflow (s_divide
);
4515 return scm_from_double ((double) xx
/ (double) yy
);
4518 else if (xx
% yy
!= 0)
4521 return scm_from_double ((double) xx
/ (double) yy
);
4522 else return scm_i_make_ratio (x
, y
);
4527 if (SCM_FIXABLE (z
))
4528 return SCM_I_MAKINUM (z
);
4530 return scm_i_long2big (z
);
4533 else if (SCM_BIGP (y
))
4536 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4537 else return scm_i_make_ratio (x
, y
);
4539 else if (SCM_REALP (y
))
4541 double yy
= SCM_REAL_VALUE (y
);
4542 #ifndef ALLOW_DIVIDE_BY_ZERO
4544 scm_num_overflow (s_divide
);
4547 return scm_from_double ((double) xx
/ yy
);
4549 else if (SCM_COMPLEXP (y
))
4552 complex_div
: /* y _must_ be a complex number */
4554 double r
= SCM_COMPLEX_REAL (y
);
4555 double i
= SCM_COMPLEX_IMAG (y
);
4559 double d
= i
* (1.0 + t
* t
);
4560 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4565 double d
= r
* (1.0 + t
* t
);
4566 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4570 else if (SCM_FRACTIONP (y
))
4571 /* a / b/c = ac / b */
4572 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4573 SCM_FRACTION_NUMERATOR (y
));
4575 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4577 else if (SCM_BIGP (x
))
4579 if (SCM_I_INUMP (y
))
4581 long int yy
= SCM_I_INUM (y
);
4584 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4585 scm_num_overflow (s_divide
);
4587 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4588 scm_remember_upto_here_1 (x
);
4589 return (sgn
== 0) ? scm_nan () : scm_inf ();
4596 /* FIXME: HMM, what are the relative performance issues here?
4597 We need to test. Is it faster on average to test
4598 divisible_p, then perform whichever operation, or is it
4599 faster to perform the integer div opportunistically and
4600 switch to real if there's a remainder? For now we take the
4601 middle ground: test, then if divisible, use the faster div
4604 long abs_yy
= yy
< 0 ? -yy
: yy
;
4605 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4609 SCM result
= scm_i_mkbig ();
4610 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4611 scm_remember_upto_here_1 (x
);
4613 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4614 return scm_i_normbig (result
);
4619 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4620 else return scm_i_make_ratio (x
, y
);
4624 else if (SCM_BIGP (y
))
4626 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4629 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4630 scm_num_overflow (s_divide
);
4632 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4633 scm_remember_upto_here_1 (x
);
4634 return (sgn
== 0) ? scm_nan () : scm_inf ();
4640 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4644 SCM result
= scm_i_mkbig ();
4645 mpz_divexact (SCM_I_BIG_MPZ (result
),
4648 scm_remember_upto_here_2 (x
, y
);
4649 return scm_i_normbig (result
);
4655 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4656 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4657 scm_remember_upto_here_2 (x
, y
);
4658 return scm_from_double (dbx
/ dby
);
4660 else return scm_i_make_ratio (x
, y
);
4664 else if (SCM_REALP (y
))
4666 double yy
= SCM_REAL_VALUE (y
);
4667 #ifndef ALLOW_DIVIDE_BY_ZERO
4669 scm_num_overflow (s_divide
);
4672 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4674 else if (SCM_COMPLEXP (y
))
4676 a
= scm_i_big2dbl (x
);
4679 else if (SCM_FRACTIONP (y
))
4680 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4681 SCM_FRACTION_NUMERATOR (y
));
4683 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4685 else if (SCM_REALP (x
))
4687 double rx
= SCM_REAL_VALUE (x
);
4688 if (SCM_I_INUMP (y
))
4690 long int yy
= SCM_I_INUM (y
);
4691 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4693 scm_num_overflow (s_divide
);
4696 return scm_from_double (rx
/ (double) yy
);
4698 else if (SCM_BIGP (y
))
4700 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4701 scm_remember_upto_here_1 (y
);
4702 return scm_from_double (rx
/ dby
);
4704 else if (SCM_REALP (y
))
4706 double yy
= SCM_REAL_VALUE (y
);
4707 #ifndef ALLOW_DIVIDE_BY_ZERO
4709 scm_num_overflow (s_divide
);
4712 return scm_from_double (rx
/ yy
);
4714 else if (SCM_COMPLEXP (y
))
4719 else if (SCM_FRACTIONP (y
))
4720 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4722 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4724 else if (SCM_COMPLEXP (x
))
4726 double rx
= SCM_COMPLEX_REAL (x
);
4727 double ix
= SCM_COMPLEX_IMAG (x
);
4728 if (SCM_I_INUMP (y
))
4730 long int yy
= SCM_I_INUM (y
);
4731 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4733 scm_num_overflow (s_divide
);
4738 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4741 else if (SCM_BIGP (y
))
4743 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4744 scm_remember_upto_here_1 (y
);
4745 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4747 else if (SCM_REALP (y
))
4749 double yy
= SCM_REAL_VALUE (y
);
4750 #ifndef ALLOW_DIVIDE_BY_ZERO
4752 scm_num_overflow (s_divide
);
4755 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4757 else if (SCM_COMPLEXP (y
))
4759 double ry
= SCM_COMPLEX_REAL (y
);
4760 double iy
= SCM_COMPLEX_IMAG (y
);
4764 double d
= iy
* (1.0 + t
* t
);
4765 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4770 double d
= ry
* (1.0 + t
* t
);
4771 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4774 else if (SCM_FRACTIONP (y
))
4776 double yy
= scm_i_fraction2double (y
);
4777 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4780 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4782 else if (SCM_FRACTIONP (x
))
4784 if (SCM_I_INUMP (y
))
4786 long int yy
= SCM_I_INUM (y
);
4787 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4789 scm_num_overflow (s_divide
);
4792 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4793 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4795 else if (SCM_BIGP (y
))
4797 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4798 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4800 else if (SCM_REALP (y
))
4802 double yy
= SCM_REAL_VALUE (y
);
4803 #ifndef ALLOW_DIVIDE_BY_ZERO
4805 scm_num_overflow (s_divide
);
4808 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4810 else if (SCM_COMPLEXP (y
))
4812 a
= scm_i_fraction2double (x
);
4815 else if (SCM_FRACTIONP (y
))
4816 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4817 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4819 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4822 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4826 scm_divide (SCM x
, SCM y
)
4828 return scm_i_divide (x
, y
, 0);
4831 static SCM
scm_divide2real (SCM x
, SCM y
)
4833 return scm_i_divide (x
, y
, 1);
4839 scm_asinh (double x
)
4844 #define asinh scm_asinh
4845 return log (x
+ sqrt (x
* x
+ 1));
4848 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4849 /* "Return the inverse hyperbolic sine of @var{x}."
4854 scm_acosh (double x
)
4859 #define acosh scm_acosh
4860 return log (x
+ sqrt (x
* x
- 1));
4863 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4864 /* "Return the inverse hyperbolic cosine of @var{x}."
4869 scm_atanh (double x
)
4874 #define atanh scm_atanh
4875 return 0.5 * log ((1 + x
) / (1 - x
));
4878 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4879 /* "Return the inverse hyperbolic tangent of @var{x}."
4884 scm_c_truncate (double x
)
4895 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4896 half-way case (ie. when x is an integer plus 0.5) going upwards.
4897 Then half-way cases are identified and adjusted down if the
4898 round-upwards didn't give the desired even integer.
4900 "plus_half == result" identifies a half-way case. If plus_half, which is
4901 x + 0.5, is an integer then x must be an integer plus 0.5.
4903 An odd "result" value is identified with result/2 != floor(result/2).
4904 This is done with plus_half, since that value is ready for use sooner in
4905 a pipelined cpu, and we're already requiring plus_half == result.
4907 Note however that we need to be careful when x is big and already an
4908 integer. In that case "x+0.5" may round to an adjacent integer, causing
4909 us to return such a value, incorrectly. For instance if the hardware is
4910 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4911 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4912 returned. Or if the hardware is in round-upwards mode, then other bigger
4913 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4914 representable value, 2^128+2^76 (or whatever), again incorrect.
4916 These bad roundings of x+0.5 are avoided by testing at the start whether
4917 x is already an integer. If it is then clearly that's the desired result
4918 already. And if it's not then the exponent must be small enough to allow
4919 an 0.5 to be represented, and hence added without a bad rounding. */
4922 scm_c_round (double x
)
4924 double plus_half
, result
;
4929 plus_half
= x
+ 0.5;
4930 result
= floor (plus_half
);
4931 /* Adjust so that the rounding is towards even. */
4932 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4937 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4939 "Round the number @var{x} towards zero.")
4940 #define FUNC_NAME s_scm_truncate_number
4942 if (scm_is_false (scm_negative_p (x
)))
4943 return scm_floor (x
);
4945 return scm_ceiling (x
);
4949 static SCM exactly_one_half
;
4951 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4953 "Round the number @var{x} towards the nearest integer. "
4954 "When it is exactly halfway between two integers, "
4955 "round towards the even one.")
4956 #define FUNC_NAME s_scm_round_number
4958 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4960 else if (SCM_REALP (x
))
4961 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
4964 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
4965 single quotient+remainder division then examining to see which way
4966 the rounding should go. */
4967 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4968 SCM result
= scm_floor (plus_half
);
4969 /* Adjust so that the rounding is towards even. */
4970 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
4971 && scm_is_true (scm_odd_p (result
)))
4972 return scm_difference (result
, SCM_I_MAKINUM (1));
4979 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4981 "Round the number @var{x} towards minus infinity.")
4982 #define FUNC_NAME s_scm_floor
4984 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4986 else if (SCM_REALP (x
))
4987 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
4988 else if (SCM_FRACTIONP (x
))
4990 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4991 SCM_FRACTION_DENOMINATOR (x
));
4992 if (scm_is_false (scm_negative_p (x
)))
4994 /* For positive x, rounding towards zero is correct. */
4999 /* For negative x, we need to return q-1 unless x is an
5000 integer. But fractions are never integer, per our
5002 return scm_difference (q
, SCM_I_MAKINUM (1));
5006 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5010 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5012 "Round the number @var{x} towards infinity.")
5013 #define FUNC_NAME s_scm_ceiling
5015 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5017 else if (SCM_REALP (x
))
5018 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5019 else if (SCM_FRACTIONP (x
))
5021 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5022 SCM_FRACTION_DENOMINATOR (x
));
5023 if (scm_is_false (scm_positive_p (x
)))
5025 /* For negative x, rounding towards zero is correct. */
5030 /* For positive x, we need to return q+1 unless x is an
5031 integer. But fractions are never integer, per our
5033 return scm_sum (q
, SCM_I_MAKINUM (1));
5037 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5041 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5042 /* "Return the square root of the real number @var{x}."
5044 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5045 /* "Return the absolute value of the real number @var{x}."
5047 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5048 /* "Return the @var{x}th power of e."
5050 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5051 /* "Return the natural logarithm of the real number @var{x}."
5053 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5054 /* "Return the sine of the real number @var{x}."
5056 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5057 /* "Return the cosine of the real number @var{x}."
5059 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5060 /* "Return the tangent of the real number @var{x}."
5062 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5063 /* "Return the arc sine of the real number @var{x}."
5065 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5066 /* "Return the arc cosine of the real number @var{x}."
5068 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5069 /* "Return the arc tangent of the real number @var{x}."
5071 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5072 /* "Return the hyperbolic sine of the real number @var{x}."
5074 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5075 /* "Return the hyperbolic cosine of the real number @var{x}."
5077 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5078 /* "Return the hyperbolic tangent of the real number @var{x}."
5086 static void scm_two_doubles (SCM x
,
5088 const char *sstring
,
5092 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5094 if (SCM_I_INUMP (x
))
5095 xy
->x
= SCM_I_INUM (x
);
5096 else if (SCM_BIGP (x
))
5097 xy
->x
= scm_i_big2dbl (x
);
5098 else if (SCM_REALP (x
))
5099 xy
->x
= SCM_REAL_VALUE (x
);
5100 else if (SCM_FRACTIONP (x
))
5101 xy
->x
= scm_i_fraction2double (x
);
5103 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5105 if (SCM_I_INUMP (y
))
5106 xy
->y
= SCM_I_INUM (y
);
5107 else if (SCM_BIGP (y
))
5108 xy
->y
= scm_i_big2dbl (y
);
5109 else if (SCM_REALP (y
))
5110 xy
->y
= SCM_REAL_VALUE (y
);
5111 else if (SCM_FRACTIONP (y
))
5112 xy
->y
= scm_i_fraction2double (y
);
5114 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5118 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5120 "Return @var{x} raised to the power of @var{y}. This\n"
5121 "procedure does not accept complex arguments.")
5122 #define FUNC_NAME s_scm_sys_expt
5125 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5126 return scm_from_double (pow (xy
.x
, xy
.y
));
5131 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5133 "Return the arc tangent of the two arguments @var{x} and\n"
5134 "@var{y}. This is similar to calculating the arc tangent of\n"
5135 "@var{x} / @var{y}, except that the signs of both arguments\n"
5136 "are used to determine the quadrant of the result. This\n"
5137 "procedure does not accept complex arguments.")
5138 #define FUNC_NAME s_scm_sys_atan2
5141 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5142 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5147 scm_c_make_rectangular (double re
, double im
)
5150 return scm_from_double (re
);
5154 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5156 SCM_COMPLEX_REAL (z
) = re
;
5157 SCM_COMPLEX_IMAG (z
) = im
;
5162 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5163 (SCM real
, SCM imaginary
),
5164 "Return a complex number constructed of the given @var{real} and\n"
5165 "@var{imaginary} parts.")
5166 #define FUNC_NAME s_scm_make_rectangular
5169 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5170 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5175 scm_c_make_polar (double mag
, double ang
)
5179 sincos (ang
, &s
, &c
);
5184 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5187 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5189 "Return the complex number @var{x} * e^(i * @var{y}).")
5190 #define FUNC_NAME s_scm_make_polar
5193 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5194 return scm_c_make_polar (xy
.x
, xy
.y
);
5199 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5200 /* "Return the real part of the number @var{z}."
5203 scm_real_part (SCM z
)
5205 if (SCM_I_INUMP (z
))
5207 else if (SCM_BIGP (z
))
5209 else if (SCM_REALP (z
))
5211 else if (SCM_COMPLEXP (z
))
5212 return scm_from_double (SCM_COMPLEX_REAL (z
));
5213 else if (SCM_FRACTIONP (z
))
5216 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5220 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5221 /* "Return the imaginary part of the number @var{z}."
5224 scm_imag_part (SCM z
)
5226 if (SCM_I_INUMP (z
))
5228 else if (SCM_BIGP (z
))
5230 else if (SCM_REALP (z
))
5232 else if (SCM_COMPLEXP (z
))
5233 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5234 else if (SCM_FRACTIONP (z
))
5237 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5240 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5241 /* "Return the numerator of the number @var{z}."
5244 scm_numerator (SCM z
)
5246 if (SCM_I_INUMP (z
))
5248 else if (SCM_BIGP (z
))
5250 else if (SCM_FRACTIONP (z
))
5252 scm_i_fraction_reduce (z
);
5253 return SCM_FRACTION_NUMERATOR (z
);
5255 else if (SCM_REALP (z
))
5256 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5258 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5262 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5263 /* "Return the denominator of the number @var{z}."
5266 scm_denominator (SCM z
)
5268 if (SCM_I_INUMP (z
))
5269 return SCM_I_MAKINUM (1);
5270 else if (SCM_BIGP (z
))
5271 return SCM_I_MAKINUM (1);
5272 else if (SCM_FRACTIONP (z
))
5274 scm_i_fraction_reduce (z
);
5275 return SCM_FRACTION_DENOMINATOR (z
);
5277 else if (SCM_REALP (z
))
5278 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5280 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5283 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5284 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5285 * "@code{abs} for real arguments, but also allows complex numbers."
5288 scm_magnitude (SCM z
)
5290 if (SCM_I_INUMP (z
))
5292 long int zz
= SCM_I_INUM (z
);
5295 else if (SCM_POSFIXABLE (-zz
))
5296 return SCM_I_MAKINUM (-zz
);
5298 return scm_i_long2big (-zz
);
5300 else if (SCM_BIGP (z
))
5302 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5303 scm_remember_upto_here_1 (z
);
5305 return scm_i_clonebig (z
, 0);
5309 else if (SCM_REALP (z
))
5310 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5311 else if (SCM_COMPLEXP (z
))
5312 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5313 else if (SCM_FRACTIONP (z
))
5315 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5317 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5318 SCM_FRACTION_DENOMINATOR (z
));
5321 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5325 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5326 /* "Return the angle of the complex number @var{z}."
5331 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5332 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5333 But if atan2 follows the floating point rounding mode, then the value
5334 is not a constant. Maybe it'd be close enough though. */
5335 if (SCM_I_INUMP (z
))
5337 if (SCM_I_INUM (z
) >= 0)
5340 return scm_from_double (atan2 (0.0, -1.0));
5342 else if (SCM_BIGP (z
))
5344 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5345 scm_remember_upto_here_1 (z
);
5347 return scm_from_double (atan2 (0.0, -1.0));
5351 else if (SCM_REALP (z
))
5353 if (SCM_REAL_VALUE (z
) >= 0)
5356 return scm_from_double (atan2 (0.0, -1.0));
5358 else if (SCM_COMPLEXP (z
))
5359 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5360 else if (SCM_FRACTIONP (z
))
5362 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5364 else return scm_from_double (atan2 (0.0, -1.0));
5367 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5371 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5372 /* Convert the number @var{x} to its inexact representation.\n"
5375 scm_exact_to_inexact (SCM z
)
5377 if (SCM_I_INUMP (z
))
5378 return scm_from_double ((double) SCM_I_INUM (z
));
5379 else if (SCM_BIGP (z
))
5380 return scm_from_double (scm_i_big2dbl (z
));
5381 else if (SCM_FRACTIONP (z
))
5382 return scm_from_double (scm_i_fraction2double (z
));
5383 else if (SCM_INEXACTP (z
))
5386 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5390 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5392 "Return an exact number that is numerically closest to @var{z}.")
5393 #define FUNC_NAME s_scm_inexact_to_exact
5395 if (SCM_I_INUMP (z
))
5397 else if (SCM_BIGP (z
))
5399 else if (SCM_REALP (z
))
5401 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5402 SCM_OUT_OF_RANGE (1, z
);
5409 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5410 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5411 scm_i_mpz2num (mpq_denref (frac
)));
5413 /* When scm_i_make_ratio throws, we leak the memory allocated
5420 else if (SCM_FRACTIONP (z
))
5423 SCM_WRONG_TYPE_ARG (1, z
);
5427 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5429 "Return an exact number that is within @var{err} of @var{x}.")
5430 #define FUNC_NAME s_scm_rationalize
5432 if (SCM_I_INUMP (x
))
5434 else if (SCM_BIGP (x
))
5436 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5438 /* Use continued fractions to find closest ratio. All
5439 arithmetic is done with exact numbers.
5442 SCM ex
= scm_inexact_to_exact (x
);
5443 SCM int_part
= scm_floor (ex
);
5444 SCM tt
= SCM_I_MAKINUM (1);
5445 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5446 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5450 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5453 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5454 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5456 /* We stop after a million iterations just to be absolutely sure
5457 that we don't go into an infinite loop. The process normally
5458 converges after less than a dozen iterations.
5461 err
= scm_abs (err
);
5462 while (++i
< 1000000)
5464 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5465 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5466 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5468 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5469 err
))) /* abs(x-a/b) <= err */
5471 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5472 if (scm_is_false (scm_exact_p (x
))
5473 || scm_is_false (scm_exact_p (err
)))
5474 return scm_exact_to_inexact (res
);
5478 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5480 tt
= scm_floor (rx
); /* tt = floor (rx) */
5486 scm_num_overflow (s_scm_rationalize
);
5489 SCM_WRONG_TYPE_ARG (1, x
);
5493 /* conversion functions */
5496 scm_is_integer (SCM val
)
5498 return scm_is_true (scm_integer_p (val
));
5502 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5504 if (SCM_I_INUMP (val
))
5506 scm_t_signed_bits n
= SCM_I_INUM (val
);
5507 return n
>= min
&& n
<= max
;
5509 else if (SCM_BIGP (val
))
5511 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5513 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5515 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5517 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5518 return n
>= min
&& n
<= max
;
5528 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5529 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5532 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5533 SCM_I_BIG_MPZ (val
));
5535 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5547 return n
>= min
&& n
<= max
;
5555 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5557 if (SCM_I_INUMP (val
))
5559 scm_t_signed_bits n
= SCM_I_INUM (val
);
5560 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5562 else if (SCM_BIGP (val
))
5564 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5566 else if (max
<= ULONG_MAX
)
5568 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5570 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5571 return n
>= min
&& n
<= max
;
5581 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5584 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5585 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5588 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5589 SCM_I_BIG_MPZ (val
));
5591 return n
>= min
&& n
<= max
;
5598 #define TYPE scm_t_intmax
5599 #define TYPE_MIN min
5600 #define TYPE_MAX max
5601 #define SIZEOF_TYPE 0
5602 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5603 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5604 #include "libguile/conv-integer.i.c"
5606 #define TYPE scm_t_uintmax
5607 #define TYPE_MIN min
5608 #define TYPE_MAX max
5609 #define SIZEOF_TYPE 0
5610 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5611 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5612 #include "libguile/conv-uinteger.i.c"
5614 #define TYPE scm_t_int8
5615 #define TYPE_MIN SCM_T_INT8_MIN
5616 #define TYPE_MAX SCM_T_INT8_MAX
5617 #define SIZEOF_TYPE 1
5618 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5619 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5620 #include "libguile/conv-integer.i.c"
5622 #define TYPE scm_t_uint8
5624 #define TYPE_MAX SCM_T_UINT8_MAX
5625 #define SIZEOF_TYPE 1
5626 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5627 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5628 #include "libguile/conv-uinteger.i.c"
5630 #define TYPE scm_t_int16
5631 #define TYPE_MIN SCM_T_INT16_MIN
5632 #define TYPE_MAX SCM_T_INT16_MAX
5633 #define SIZEOF_TYPE 2
5634 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5635 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5636 #include "libguile/conv-integer.i.c"
5638 #define TYPE scm_t_uint16
5640 #define TYPE_MAX SCM_T_UINT16_MAX
5641 #define SIZEOF_TYPE 2
5642 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5643 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5644 #include "libguile/conv-uinteger.i.c"
5646 #define TYPE scm_t_int32
5647 #define TYPE_MIN SCM_T_INT32_MIN
5648 #define TYPE_MAX SCM_T_INT32_MAX
5649 #define SIZEOF_TYPE 4
5650 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5651 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5652 #include "libguile/conv-integer.i.c"
5654 #define TYPE scm_t_uint32
5656 #define TYPE_MAX SCM_T_UINT32_MAX
5657 #define SIZEOF_TYPE 4
5658 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5659 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5660 #include "libguile/conv-uinteger.i.c"
5662 #if SCM_HAVE_T_INT64
5664 #define TYPE scm_t_int64
5665 #define TYPE_MIN SCM_T_INT64_MIN
5666 #define TYPE_MAX SCM_T_INT64_MAX
5667 #define SIZEOF_TYPE 8
5668 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5669 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5670 #include "libguile/conv-integer.i.c"
5672 #define TYPE scm_t_uint64
5674 #define TYPE_MAX SCM_T_UINT64_MAX
5675 #define SIZEOF_TYPE 8
5676 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5677 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5678 #include "libguile/conv-uinteger.i.c"
5683 scm_to_mpz (SCM val
, mpz_t rop
)
5685 if (SCM_I_INUMP (val
))
5686 mpz_set_si (rop
, SCM_I_INUM (val
));
5687 else if (SCM_BIGP (val
))
5688 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5690 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5694 scm_from_mpz (mpz_t val
)
5696 return scm_i_mpz2num (val
);
5700 scm_is_real (SCM val
)
5702 return scm_is_true (scm_real_p (val
));
5706 scm_is_rational (SCM val
)
5708 return scm_is_true (scm_rational_p (val
));
5712 scm_to_double (SCM val
)
5714 if (SCM_I_INUMP (val
))
5715 return SCM_I_INUM (val
);
5716 else if (SCM_BIGP (val
))
5717 return scm_i_big2dbl (val
);
5718 else if (SCM_FRACTIONP (val
))
5719 return scm_i_fraction2double (val
);
5720 else if (SCM_REALP (val
))
5721 return SCM_REAL_VALUE (val
);
5723 scm_wrong_type_arg (NULL
, 0, val
);
5727 scm_from_double (double val
)
5729 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5730 SCM_REAL_VALUE (z
) = val
;
5734 #if SCM_ENABLE_DISCOURAGED == 1
5737 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5741 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5745 scm_out_of_range (NULL
, num
);
5748 return scm_to_double (num
);
5752 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5756 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5760 scm_out_of_range (NULL
, num
);
5763 return scm_to_double (num
);
5769 scm_is_complex (SCM val
)
5771 return scm_is_true (scm_complex_p (val
));
5775 scm_c_real_part (SCM z
)
5777 if (SCM_COMPLEXP (z
))
5778 return SCM_COMPLEX_REAL (z
);
5781 /* Use the scm_real_part to get proper error checking and
5784 return scm_to_double (scm_real_part (z
));
5789 scm_c_imag_part (SCM z
)
5791 if (SCM_COMPLEXP (z
))
5792 return SCM_COMPLEX_IMAG (z
);
5795 /* Use the scm_imag_part to get proper error checking and
5796 dispatching. The result will almost always be 0.0, but not
5799 return scm_to_double (scm_imag_part (z
));
5804 scm_c_magnitude (SCM z
)
5806 return scm_to_double (scm_magnitude (z
));
5812 return scm_to_double (scm_angle (z
));
5816 scm_is_number (SCM z
)
5818 return scm_is_true (scm_number_p (z
));
5826 mpz_init_set_si (z_negative_one
, -1);
5828 /* It may be possible to tune the performance of some algorithms by using
5829 * the following constants to avoid the creation of bignums. Please, before
5830 * using these values, remember the two rules of program optimization:
5831 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5832 scm_c_define ("most-positive-fixnum",
5833 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5834 scm_c_define ("most-negative-fixnum",
5835 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5837 scm_add_feature ("complex");
5838 scm_add_feature ("inexact");
5839 scm_flo0
= scm_from_double (0.0);
5841 /* determine floating point precision */
5842 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5844 init_dblprec(&scm_dblprec
[i
-2],i
);
5845 init_fx_radix(fx_per_radix
[i
-2],i
);
5848 /* hard code precision for base 10 if the preprocessor tells us to... */
5849 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5852 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5853 SCM_I_MAKINUM (2)));
5854 #include "libguile/numbers.x"