1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005 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 SCM
187 scm_i_clonebig (SCM src_big
, int same_sign_p
)
189 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
190 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
191 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
193 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
197 SCM_C_INLINE_KEYWORD
int
198 scm_i_bigcmp (SCM x
, SCM y
)
200 /* Return neg if x < y, pos if x > y, and 0 if x == y */
201 /* presume we already know x and y are bignums */
202 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
203 scm_remember_upto_here_2 (x
, y
);
207 SCM_C_INLINE_KEYWORD SCM
208 scm_i_dbl2big (double d
)
210 /* results are only defined if d is an integer */
211 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
212 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
216 /* Convert a integer in double representation to a SCM number. */
218 SCM_C_INLINE_KEYWORD SCM
219 scm_i_dbl2num (double u
)
221 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
222 powers of 2, so there's no rounding when making "double" values
223 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
224 get rounded on a 64-bit machine, hence the "+1".
226 The use of floor() to force to an integer value ensures we get a
227 "numerically closest" value without depending on how a
228 double->long cast or how mpz_set_d will round. For reference,
229 double->long probably follows the hardware rounding mode,
230 mpz_set_d truncates towards zero. */
232 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
233 representable as a double? */
235 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
236 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
237 return SCM_I_MAKINUM ((long) u
);
239 return scm_i_dbl2big (u
);
242 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
243 with R5RS exact->inexact.
245 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
246 (ie. truncate towards zero), then adjust to get the closest double by
247 examining the next lower bit and adding 1 (to the absolute value) if
250 Bignums exactly half way between representable doubles are rounded to the
251 next higher absolute value (ie. away from zero). This seems like an
252 adequate interpretation of R5RS "numerically closest", and it's easier
253 and faster than a full "nearest-even" style.
255 The bit test must be done on the absolute value of the mpz_t, which means
256 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
257 negatives as twos complement.
259 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
260 following the hardware rounding mode, but applied to the absolute value
261 of the mpz_t operand. This is not what we want so we put the high
262 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
263 mpz_get_d is supposed to always truncate towards zero.
265 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
266 is a slowdown. It'd be faster to pick out the relevant high bits with
267 mpz_getlimbn if we could be bothered coding that, and if the new
268 truncating gmp doesn't come out. */
271 scm_i_big2dbl (SCM b
)
276 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
280 /* Current GMP, eg. 4.1.3, force truncation towards zero */
282 if (bits
> DBL_MANT_DIG
)
284 size_t shift
= bits
- DBL_MANT_DIG
;
285 mpz_init2 (tmp
, DBL_MANT_DIG
);
286 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
287 result
= ldexp (mpz_get_d (tmp
), shift
);
292 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
297 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
300 if (bits
> DBL_MANT_DIG
)
302 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
303 /* test bit number "pos" in absolute value */
304 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
305 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
307 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
311 scm_remember_upto_here_1 (b
);
315 SCM_C_INLINE_KEYWORD SCM
316 scm_i_normbig (SCM b
)
318 /* convert a big back to a fixnum if it'll fit */
319 /* presume b is a bignum */
320 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
322 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
323 if (SCM_FIXABLE (val
))
324 b
= SCM_I_MAKINUM (val
);
329 static SCM_C_INLINE_KEYWORD SCM
330 scm_i_mpz2num (mpz_t b
)
332 /* convert a mpz number to a SCM number. */
333 if (mpz_fits_slong_p (b
))
335 long val
= mpz_get_si (b
);
336 if (SCM_FIXABLE (val
))
337 return SCM_I_MAKINUM (val
);
341 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
342 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
347 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
348 static SCM
scm_divide2real (SCM x
, SCM y
);
351 scm_i_make_ratio (SCM numerator
, SCM denominator
)
352 #define FUNC_NAME "make-ratio"
354 /* First make sure the arguments are proper.
356 if (SCM_I_INUMP (denominator
))
358 if (scm_is_eq (denominator
, SCM_INUM0
))
359 scm_num_overflow ("make-ratio");
360 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
365 if (!(SCM_BIGP(denominator
)))
366 SCM_WRONG_TYPE_ARG (2, denominator
);
368 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
369 SCM_WRONG_TYPE_ARG (1, numerator
);
371 /* Then flip signs so that the denominator is positive.
373 if (scm_is_true (scm_negative_p (denominator
)))
375 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
376 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
379 /* Now consider for each of the four fixnum/bignum combinations
380 whether the rational number is really an integer.
382 if (SCM_I_INUMP (numerator
))
384 long x
= SCM_I_INUM (numerator
);
385 if (scm_is_eq (numerator
, SCM_INUM0
))
387 if (SCM_I_INUMP (denominator
))
390 y
= SCM_I_INUM (denominator
);
392 return SCM_I_MAKINUM(1);
394 return SCM_I_MAKINUM (x
/ y
);
398 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
399 of that value for the denominator, as a bignum. Apart from
400 that case, abs(bignum) > abs(inum) so inum/bignum is not an
402 if (x
== SCM_MOST_NEGATIVE_FIXNUM
403 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
404 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
405 return SCM_I_MAKINUM(-1);
408 else if (SCM_BIGP (numerator
))
410 if (SCM_I_INUMP (denominator
))
412 long yy
= SCM_I_INUM (denominator
);
413 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
414 return scm_divide (numerator
, denominator
);
418 if (scm_is_eq (numerator
, denominator
))
419 return SCM_I_MAKINUM(1);
420 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
421 SCM_I_BIG_MPZ (denominator
)))
422 return scm_divide(numerator
, denominator
);
426 /* No, it's a proper fraction.
428 return scm_double_cell (scm_tc16_fraction
,
429 SCM_UNPACK (numerator
),
430 SCM_UNPACK (denominator
), 0);
434 static void scm_i_fraction_reduce (SCM z
)
436 if (!(SCM_FRACTION_REDUCED (z
)))
439 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
440 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
443 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
444 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
446 SCM_FRACTION_REDUCED_SET (z
);
451 scm_i_fraction2double (SCM z
)
453 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
454 SCM_FRACTION_DENOMINATOR (z
)));
457 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
459 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
461 #define FUNC_NAME s_scm_exact_p
467 if (SCM_FRACTIONP (x
))
471 SCM_WRONG_TYPE_ARG (1, x
);
476 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
478 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
480 #define FUNC_NAME s_scm_odd_p
484 long val
= SCM_I_INUM (n
);
485 return scm_from_bool ((val
& 1L) != 0);
487 else if (SCM_BIGP (n
))
489 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
490 scm_remember_upto_here_1 (n
);
491 return scm_from_bool (odd_p
);
493 else if (scm_is_true (scm_inf_p (n
)))
495 else if (SCM_REALP (n
))
497 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
503 SCM_WRONG_TYPE_ARG (1, n
);
506 SCM_WRONG_TYPE_ARG (1, n
);
511 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
513 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
515 #define FUNC_NAME s_scm_even_p
519 long val
= SCM_I_INUM (n
);
520 return scm_from_bool ((val
& 1L) == 0);
522 else if (SCM_BIGP (n
))
524 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
525 scm_remember_upto_here_1 (n
);
526 return scm_from_bool (even_p
);
528 else if (scm_is_true (scm_inf_p (n
)))
530 else if (SCM_REALP (n
))
532 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
538 SCM_WRONG_TYPE_ARG (1, n
);
541 SCM_WRONG_TYPE_ARG (1, n
);
545 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
547 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
548 "or @samp{-inf.0}, @code{#f} otherwise.")
549 #define FUNC_NAME s_scm_inf_p
552 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
553 else if (SCM_COMPLEXP (x
))
554 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
555 || xisinf (SCM_COMPLEX_IMAG (x
)));
561 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
563 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
565 #define FUNC_NAME s_scm_nan_p
568 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
569 else if (SCM_COMPLEXP (n
))
570 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
571 || xisnan (SCM_COMPLEX_IMAG (n
)));
577 /* Guile's idea of infinity. */
578 static double guile_Inf
;
580 /* Guile's idea of not a number. */
581 static double guile_NaN
;
584 guile_ieee_init (void)
586 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
588 /* Some version of gcc on some old version of Linux used to crash when
589 trying to make Inf and NaN. */
592 /* C99 INFINITY, when available.
593 FIXME: The standard allows for INFINITY to be something that overflows
594 at compile time. We ought to have a configure test to check for that
595 before trying to use it. (But in practice we believe this is not a
596 problem on any system guile is likely to target.) */
597 guile_Inf
= INFINITY
;
600 extern unsigned int DINFINITY
[2];
601 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
608 if (guile_Inf
== tmp
)
616 #if defined (HAVE_ISNAN)
619 /* C99 NAN, when available */
623 extern unsigned int DQNAN
[2];
624 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
626 guile_NaN
= guile_Inf
/ guile_Inf
;
632 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
635 #define FUNC_NAME s_scm_inf
637 static int initialized
= 0;
643 return scm_from_double (guile_Inf
);
647 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
650 #define FUNC_NAME s_scm_nan
652 static int initialized
= 0;
658 return scm_from_double (guile_NaN
);
663 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
665 "Return the absolute value of @var{x}.")
670 long int xx
= SCM_I_INUM (x
);
673 else if (SCM_POSFIXABLE (-xx
))
674 return SCM_I_MAKINUM (-xx
);
676 return scm_i_long2big (-xx
);
678 else if (SCM_BIGP (x
))
680 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
682 return scm_i_clonebig (x
, 0);
686 else if (SCM_REALP (x
))
688 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
689 double xx
= SCM_REAL_VALUE (x
);
691 return scm_from_double (-xx
);
695 else if (SCM_FRACTIONP (x
))
697 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
699 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
700 SCM_FRACTION_DENOMINATOR (x
));
703 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
708 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
709 /* "Return the quotient of the numbers @var{x} and @var{y}."
712 scm_quotient (SCM x
, SCM y
)
716 long xx
= SCM_I_INUM (x
);
719 long yy
= SCM_I_INUM (y
);
721 scm_num_overflow (s_quotient
);
726 return SCM_I_MAKINUM (z
);
728 return scm_i_long2big (z
);
731 else if (SCM_BIGP (y
))
733 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
734 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
735 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
737 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
738 scm_remember_upto_here_1 (y
);
739 return SCM_I_MAKINUM (-1);
742 return SCM_I_MAKINUM (0);
745 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
747 else if (SCM_BIGP (x
))
751 long yy
= SCM_I_INUM (y
);
753 scm_num_overflow (s_quotient
);
758 SCM result
= scm_i_mkbig ();
761 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
764 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
767 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
768 scm_remember_upto_here_1 (x
);
769 return scm_i_normbig (result
);
772 else if (SCM_BIGP (y
))
774 SCM result
= scm_i_mkbig ();
775 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
778 scm_remember_upto_here_2 (x
, y
);
779 return scm_i_normbig (result
);
782 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
785 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
788 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
789 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
791 * "(remainder 13 4) @result{} 1\n"
792 * "(remainder -13 4) @result{} -1\n"
796 scm_remainder (SCM x
, SCM y
)
802 long yy
= SCM_I_INUM (y
);
804 scm_num_overflow (s_remainder
);
807 long z
= SCM_I_INUM (x
) % yy
;
808 return SCM_I_MAKINUM (z
);
811 else if (SCM_BIGP (y
))
813 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
814 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
815 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
817 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
818 scm_remember_upto_here_1 (y
);
819 return SCM_I_MAKINUM (0);
825 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
827 else if (SCM_BIGP (x
))
831 long yy
= SCM_I_INUM (y
);
833 scm_num_overflow (s_remainder
);
836 SCM result
= scm_i_mkbig ();
839 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
840 scm_remember_upto_here_1 (x
);
841 return scm_i_normbig (result
);
844 else if (SCM_BIGP (y
))
846 SCM result
= scm_i_mkbig ();
847 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
850 scm_remember_upto_here_2 (x
, y
);
851 return scm_i_normbig (result
);
854 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
857 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
861 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
862 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
864 * "(modulo 13 4) @result{} 1\n"
865 * "(modulo -13 4) @result{} 3\n"
869 scm_modulo (SCM x
, SCM y
)
873 long xx
= SCM_I_INUM (x
);
876 long yy
= SCM_I_INUM (y
);
878 scm_num_overflow (s_modulo
);
881 /* C99 specifies that "%" is the remainder corresponding to a
882 quotient rounded towards zero, and that's also traditional
883 for machine division, so z here should be well defined. */
901 return SCM_I_MAKINUM (result
);
904 else if (SCM_BIGP (y
))
906 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
913 SCM pos_y
= scm_i_clonebig (y
, 0);
914 /* do this after the last scm_op */
915 mpz_init_set_si (z_x
, xx
);
916 result
= pos_y
; /* re-use this bignum */
917 mpz_mod (SCM_I_BIG_MPZ (result
),
919 SCM_I_BIG_MPZ (pos_y
));
920 scm_remember_upto_here_1 (pos_y
);
924 result
= scm_i_mkbig ();
925 /* do this after the last scm_op */
926 mpz_init_set_si (z_x
, xx
);
927 mpz_mod (SCM_I_BIG_MPZ (result
),
930 scm_remember_upto_here_1 (y
);
933 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
934 mpz_add (SCM_I_BIG_MPZ (result
),
936 SCM_I_BIG_MPZ (result
));
937 scm_remember_upto_here_1 (y
);
938 /* and do this before the next one */
940 return scm_i_normbig (result
);
944 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
946 else if (SCM_BIGP (x
))
950 long yy
= SCM_I_INUM (y
);
952 scm_num_overflow (s_modulo
);
955 SCM result
= scm_i_mkbig ();
956 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
958 (yy
< 0) ? - yy
: yy
);
959 scm_remember_upto_here_1 (x
);
960 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
961 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
962 SCM_I_BIG_MPZ (result
),
964 return scm_i_normbig (result
);
967 else if (SCM_BIGP (y
))
970 SCM result
= scm_i_mkbig ();
971 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
972 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
973 mpz_mod (SCM_I_BIG_MPZ (result
),
975 SCM_I_BIG_MPZ (pos_y
));
977 scm_remember_upto_here_1 (x
);
978 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
979 mpz_add (SCM_I_BIG_MPZ (result
),
981 SCM_I_BIG_MPZ (result
));
982 scm_remember_upto_here_2 (y
, pos_y
);
983 return scm_i_normbig (result
);
987 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
990 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
993 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
994 /* "Return the greatest common divisor of all arguments.\n"
995 * "If called without arguments, 0 is returned."
998 scm_gcd (SCM x
, SCM y
)
1001 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1003 if (SCM_I_INUMP (x
))
1005 if (SCM_I_INUMP (y
))
1007 long xx
= SCM_I_INUM (x
);
1008 long yy
= SCM_I_INUM (y
);
1009 long u
= xx
< 0 ? -xx
: xx
;
1010 long v
= yy
< 0 ? -yy
: yy
;
1020 /* Determine a common factor 2^k */
1021 while (!(1 & (u
| v
)))
1027 /* Now, any factor 2^n can be eliminated */
1047 return (SCM_POSFIXABLE (result
)
1048 ? SCM_I_MAKINUM (result
)
1049 : scm_i_long2big (result
));
1051 else if (SCM_BIGP (y
))
1057 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1059 else if (SCM_BIGP (x
))
1061 if (SCM_I_INUMP (y
))
1063 unsigned long result
;
1066 yy
= SCM_I_INUM (y
);
1071 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1072 scm_remember_upto_here_1 (x
);
1073 return (SCM_POSFIXABLE (result
)
1074 ? SCM_I_MAKINUM (result
)
1075 : scm_from_ulong (result
));
1077 else if (SCM_BIGP (y
))
1079 SCM result
= scm_i_mkbig ();
1080 mpz_gcd (SCM_I_BIG_MPZ (result
),
1083 scm_remember_upto_here_2 (x
, y
);
1084 return scm_i_normbig (result
);
1087 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1090 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1093 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1094 /* "Return the least common multiple of the arguments.\n"
1095 * "If called without arguments, 1 is returned."
1098 scm_lcm (SCM n1
, SCM n2
)
1100 if (SCM_UNBNDP (n2
))
1102 if (SCM_UNBNDP (n1
))
1103 return SCM_I_MAKINUM (1L);
1104 n2
= SCM_I_MAKINUM (1L);
1107 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1108 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1109 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1110 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1112 if (SCM_I_INUMP (n1
))
1114 if (SCM_I_INUMP (n2
))
1116 SCM d
= scm_gcd (n1
, n2
);
1117 if (scm_is_eq (d
, SCM_INUM0
))
1120 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1124 /* inum n1, big n2 */
1127 SCM result
= scm_i_mkbig ();
1128 long nn1
= SCM_I_INUM (n1
);
1129 if (nn1
== 0) return SCM_INUM0
;
1130 if (nn1
< 0) nn1
= - nn1
;
1131 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1132 scm_remember_upto_here_1 (n2
);
1140 if (SCM_I_INUMP (n2
))
1147 SCM result
= scm_i_mkbig ();
1148 mpz_lcm(SCM_I_BIG_MPZ (result
),
1150 SCM_I_BIG_MPZ (n2
));
1151 scm_remember_upto_here_2(n1
, n2
);
1152 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1158 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1163 + + + x (map digit:logand X Y)
1164 + - + x (map digit:logand X (lognot (+ -1 Y)))
1165 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1166 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1171 + + + (map digit:logior X Y)
1172 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1173 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1174 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1179 + + + (map digit:logxor X Y)
1180 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1181 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1182 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1187 + + (any digit:logand X Y)
1188 + - (any digit:logand X (lognot (+ -1 Y)))
1189 - + (any digit:logand (lognot (+ -1 X)) Y)
1194 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1196 "Return the bitwise AND of the integer arguments.\n\n"
1198 "(logand) @result{} -1\n"
1199 "(logand 7) @result{} 7\n"
1200 "(logand #b111 #b011 #b001) @result{} 1\n"
1202 #define FUNC_NAME s_scm_logand
1206 if (SCM_UNBNDP (n2
))
1208 if (SCM_UNBNDP (n1
))
1209 return SCM_I_MAKINUM (-1);
1210 else if (!SCM_NUMBERP (n1
))
1211 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1212 else if (SCM_NUMBERP (n1
))
1215 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1218 if (SCM_I_INUMP (n1
))
1220 nn1
= SCM_I_INUM (n1
);
1221 if (SCM_I_INUMP (n2
))
1223 long nn2
= SCM_I_INUM (n2
);
1224 return SCM_I_MAKINUM (nn1
& nn2
);
1226 else if SCM_BIGP (n2
)
1232 SCM result_z
= scm_i_mkbig ();
1234 mpz_init_set_si (nn1_z
, nn1
);
1235 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1236 scm_remember_upto_here_1 (n2
);
1238 return scm_i_normbig (result_z
);
1242 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1244 else if (SCM_BIGP (n1
))
1246 if (SCM_I_INUMP (n2
))
1249 nn1
= SCM_I_INUM (n1
);
1252 else if (SCM_BIGP (n2
))
1254 SCM result_z
= scm_i_mkbig ();
1255 mpz_and (SCM_I_BIG_MPZ (result_z
),
1257 SCM_I_BIG_MPZ (n2
));
1258 scm_remember_upto_here_2 (n1
, n2
);
1259 return scm_i_normbig (result_z
);
1262 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1265 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1270 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1272 "Return the bitwise OR of the integer arguments.\n\n"
1274 "(logior) @result{} 0\n"
1275 "(logior 7) @result{} 7\n"
1276 "(logior #b000 #b001 #b011) @result{} 3\n"
1278 #define FUNC_NAME s_scm_logior
1282 if (SCM_UNBNDP (n2
))
1284 if (SCM_UNBNDP (n1
))
1286 else if (SCM_NUMBERP (n1
))
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1292 if (SCM_I_INUMP (n1
))
1294 nn1
= SCM_I_INUM (n1
);
1295 if (SCM_I_INUMP (n2
))
1297 long nn2
= SCM_I_INUM (n2
);
1298 return SCM_I_MAKINUM (nn1
| nn2
);
1300 else if (SCM_BIGP (n2
))
1306 SCM result_z
= scm_i_mkbig ();
1308 mpz_init_set_si (nn1_z
, nn1
);
1309 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1310 scm_remember_upto_here_1 (n2
);
1312 return scm_i_normbig (result_z
);
1316 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1318 else if (SCM_BIGP (n1
))
1320 if (SCM_I_INUMP (n2
))
1323 nn1
= SCM_I_INUM (n1
);
1326 else if (SCM_BIGP (n2
))
1328 SCM result_z
= scm_i_mkbig ();
1329 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1331 SCM_I_BIG_MPZ (n2
));
1332 scm_remember_upto_here_2 (n1
, n2
);
1333 return scm_i_normbig (result_z
);
1336 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1339 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1344 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1346 "Return the bitwise XOR of the integer arguments. A bit is\n"
1347 "set in the result if it is set in an odd number of arguments.\n"
1349 "(logxor) @result{} 0\n"
1350 "(logxor 7) @result{} 7\n"
1351 "(logxor #b000 #b001 #b011) @result{} 2\n"
1352 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1354 #define FUNC_NAME s_scm_logxor
1358 if (SCM_UNBNDP (n2
))
1360 if (SCM_UNBNDP (n1
))
1362 else if (SCM_NUMBERP (n1
))
1365 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1368 if (SCM_I_INUMP (n1
))
1370 nn1
= SCM_I_INUM (n1
);
1371 if (SCM_I_INUMP (n2
))
1373 long nn2
= SCM_I_INUM (n2
);
1374 return SCM_I_MAKINUM (nn1
^ nn2
);
1376 else if (SCM_BIGP (n2
))
1380 SCM result_z
= scm_i_mkbig ();
1382 mpz_init_set_si (nn1_z
, nn1
);
1383 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1384 scm_remember_upto_here_1 (n2
);
1386 return scm_i_normbig (result_z
);
1390 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1392 else if (SCM_BIGP (n1
))
1394 if (SCM_I_INUMP (n2
))
1397 nn1
= SCM_I_INUM (n1
);
1400 else if (SCM_BIGP (n2
))
1402 SCM result_z
= scm_i_mkbig ();
1403 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1405 SCM_I_BIG_MPZ (n2
));
1406 scm_remember_upto_here_2 (n1
, n2
);
1407 return scm_i_normbig (result_z
);
1410 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1413 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1418 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1420 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1421 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1422 "without actually calculating the @code{logand}, just testing\n"
1426 "(logtest #b0100 #b1011) @result{} #f\n"
1427 "(logtest #b0100 #b0111) @result{} #t\n"
1429 #define FUNC_NAME s_scm_logtest
1433 if (SCM_I_INUMP (j
))
1435 nj
= SCM_I_INUM (j
);
1436 if (SCM_I_INUMP (k
))
1438 long nk
= SCM_I_INUM (k
);
1439 return scm_from_bool (nj
& nk
);
1441 else if (SCM_BIGP (k
))
1449 mpz_init_set_si (nj_z
, nj
);
1450 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1451 scm_remember_upto_here_1 (k
);
1452 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1458 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1460 else if (SCM_BIGP (j
))
1462 if (SCM_I_INUMP (k
))
1465 nj
= SCM_I_INUM (j
);
1468 else if (SCM_BIGP (k
))
1472 mpz_init (result_z
);
1476 scm_remember_upto_here_2 (j
, k
);
1477 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1478 mpz_clear (result_z
);
1482 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1485 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1490 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1492 "Test whether bit number @var{index} in @var{j} is set.\n"
1493 "@var{index} starts from 0 for the least significant bit.\n"
1496 "(logbit? 0 #b1101) @result{} #t\n"
1497 "(logbit? 1 #b1101) @result{} #f\n"
1498 "(logbit? 2 #b1101) @result{} #t\n"
1499 "(logbit? 3 #b1101) @result{} #t\n"
1500 "(logbit? 4 #b1101) @result{} #f\n"
1502 #define FUNC_NAME s_scm_logbit_p
1504 unsigned long int iindex
;
1505 iindex
= scm_to_ulong (index
);
1507 if (SCM_I_INUMP (j
))
1509 /* bits above what's in an inum follow the sign bit */
1510 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1511 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1513 else if (SCM_BIGP (j
))
1515 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1516 scm_remember_upto_here_1 (j
);
1517 return scm_from_bool (val
);
1520 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1525 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1527 "Return the integer which is the ones-complement of the integer\n"
1531 "(number->string (lognot #b10000000) 2)\n"
1532 " @result{} \"-10000001\"\n"
1533 "(number->string (lognot #b0) 2)\n"
1534 " @result{} \"-1\"\n"
1536 #define FUNC_NAME s_scm_lognot
1538 if (SCM_I_INUMP (n
)) {
1539 /* No overflow here, just need to toggle all the bits making up the inum.
1540 Enhancement: No need to strip the tag and add it back, could just xor
1541 a block of 1 bits, if that worked with the various debug versions of
1543 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1545 } else if (SCM_BIGP (n
)) {
1546 SCM result
= scm_i_mkbig ();
1547 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1548 scm_remember_upto_here_1 (n
);
1552 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1557 /* returns 0 if IN is not an integer. OUT must already be
1560 coerce_to_big (SCM in
, mpz_t out
)
1563 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1564 else if (SCM_I_INUMP (in
))
1565 mpz_set_si (out
, SCM_I_INUM (in
));
1572 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1573 (SCM n
, SCM k
, SCM m
),
1574 "Return @var{n} raised to the integer exponent\n"
1575 "@var{k}, modulo @var{m}.\n"
1578 "(modulo-expt 2 3 5)\n"
1581 #define FUNC_NAME s_scm_modulo_expt
1587 /* There are two classes of error we might encounter --
1588 1) Math errors, which we'll report by calling scm_num_overflow,
1590 2) wrong-type errors, which of course we'll report by calling
1592 We don't report those errors immediately, however; instead we do
1593 some cleanup first. These variables tell us which error (if
1594 any) we should report after cleaning up.
1596 int report_overflow
= 0;
1598 int position_of_wrong_type
= 0;
1599 SCM value_of_wrong_type
= SCM_INUM0
;
1601 SCM result
= SCM_UNDEFINED
;
1607 if (scm_is_eq (m
, SCM_INUM0
))
1609 report_overflow
= 1;
1613 if (!coerce_to_big (n
, n_tmp
))
1615 value_of_wrong_type
= n
;
1616 position_of_wrong_type
= 1;
1620 if (!coerce_to_big (k
, k_tmp
))
1622 value_of_wrong_type
= k
;
1623 position_of_wrong_type
= 2;
1627 if (!coerce_to_big (m
, m_tmp
))
1629 value_of_wrong_type
= m
;
1630 position_of_wrong_type
= 3;
1634 /* if the exponent K is negative, and we simply call mpz_powm, we
1635 will get a divide-by-zero exception when an inverse 1/n mod m
1636 doesn't exist (or is not unique). Since exceptions are hard to
1637 handle, we'll attempt the inversion "by hand" -- that way, we get
1638 a simple failure code, which is easy to handle. */
1640 if (-1 == mpz_sgn (k_tmp
))
1642 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1644 report_overflow
= 1;
1647 mpz_neg (k_tmp
, k_tmp
);
1650 result
= scm_i_mkbig ();
1651 mpz_powm (SCM_I_BIG_MPZ (result
),
1656 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1657 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1664 if (report_overflow
)
1665 scm_num_overflow (FUNC_NAME
);
1667 if (position_of_wrong_type
)
1668 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1669 value_of_wrong_type
);
1671 return scm_i_normbig (result
);
1675 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1677 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1678 "exact integer, @var{n} can be any number.\n"
1680 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1681 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1682 "includes @math{0^0} is 1.\n"
1685 "(integer-expt 2 5) @result{} 32\n"
1686 "(integer-expt -3 3) @result{} -27\n"
1687 "(integer-expt 5 -3) @result{} 1/125\n"
1688 "(integer-expt 0 0) @result{} 1\n"
1690 #define FUNC_NAME s_scm_integer_expt
1693 SCM z_i2
= SCM_BOOL_F
;
1695 SCM acc
= SCM_I_MAKINUM (1L);
1697 /* 0^0 == 1 according to R5RS */
1698 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1699 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1700 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1701 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1703 if (SCM_I_INUMP (k
))
1704 i2
= SCM_I_INUM (k
);
1705 else if (SCM_BIGP (k
))
1707 z_i2
= scm_i_clonebig (k
, 1);
1708 scm_remember_upto_here_1 (k
);
1712 SCM_WRONG_TYPE_ARG (2, k
);
1716 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1718 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1719 n
= scm_divide (n
, SCM_UNDEFINED
);
1723 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1727 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1729 return scm_product (acc
, n
);
1731 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1732 acc
= scm_product (acc
, n
);
1733 n
= scm_product (n
, n
);
1734 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1742 n
= scm_divide (n
, SCM_UNDEFINED
);
1749 return scm_product (acc
, n
);
1751 acc
= scm_product (acc
, n
);
1752 n
= scm_product (n
, n
);
1759 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1761 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1762 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1764 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1765 "@var{cnt} is negative it's a division, rounded towards negative\n"
1766 "infinity. (Note that this is not the same rounding as\n"
1767 "@code{quotient} does.)\n"
1769 "With @var{n} viewed as an infinite precision twos complement,\n"
1770 "@code{ash} means a left shift introducing zero bits, or a right\n"
1771 "shift dropping bits.\n"
1774 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1775 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1777 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1778 "(ash -23 -2) @result{} -6\n"
1780 #define FUNC_NAME s_scm_ash
1783 bits_to_shift
= scm_to_long (cnt
);
1785 if (SCM_I_INUMP (n
))
1787 long nn
= SCM_I_INUM (n
);
1789 if (bits_to_shift
> 0)
1791 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1792 overflow a non-zero fixnum. For smaller shifts we check the
1793 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1794 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1795 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1801 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1803 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1806 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1810 SCM result
= scm_i_long2big (nn
);
1811 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1818 bits_to_shift
= -bits_to_shift
;
1819 if (bits_to_shift
>= SCM_LONG_BIT
)
1820 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1822 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1826 else if (SCM_BIGP (n
))
1830 if (bits_to_shift
== 0)
1833 result
= scm_i_mkbig ();
1834 if (bits_to_shift
>= 0)
1836 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1842 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1843 we have to allocate a bignum even if the result is going to be a
1845 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1847 return scm_i_normbig (result
);
1853 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1859 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1860 (SCM n
, SCM start
, SCM end
),
1861 "Return the integer composed of the @var{start} (inclusive)\n"
1862 "through @var{end} (exclusive) bits of @var{n}. The\n"
1863 "@var{start}th bit becomes the 0-th bit in the result.\n"
1866 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1867 " @result{} \"1010\"\n"
1868 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1869 " @result{} \"10110\"\n"
1871 #define FUNC_NAME s_scm_bit_extract
1873 unsigned long int istart
, iend
, bits
;
1874 istart
= scm_to_ulong (start
);
1875 iend
= scm_to_ulong (end
);
1876 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1878 /* how many bits to keep */
1879 bits
= iend
- istart
;
1881 if (SCM_I_INUMP (n
))
1883 long int in
= SCM_I_INUM (n
);
1885 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1886 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1887 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1889 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1891 /* Since we emulate two's complement encoded numbers, this
1892 * special case requires us to produce a result that has
1893 * more bits than can be stored in a fixnum.
1895 SCM result
= scm_i_long2big (in
);
1896 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1901 /* mask down to requisite bits */
1902 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1903 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1905 else if (SCM_BIGP (n
))
1910 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1914 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1915 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1916 such bits into a ulong. */
1917 result
= scm_i_mkbig ();
1918 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1919 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1920 result
= scm_i_normbig (result
);
1922 scm_remember_upto_here_1 (n
);
1926 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1931 static const char scm_logtab
[] = {
1932 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1935 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1937 "Return the number of bits in integer @var{n}. If integer is\n"
1938 "positive, the 1-bits in its binary representation are counted.\n"
1939 "If negative, the 0-bits in its two's-complement binary\n"
1940 "representation are counted. If 0, 0 is returned.\n"
1943 "(logcount #b10101010)\n"
1950 #define FUNC_NAME s_scm_logcount
1952 if (SCM_I_INUMP (n
))
1954 unsigned long int c
= 0;
1955 long int nn
= SCM_I_INUM (n
);
1960 c
+= scm_logtab
[15 & nn
];
1963 return SCM_I_MAKINUM (c
);
1965 else if (SCM_BIGP (n
))
1967 unsigned long count
;
1968 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1969 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1971 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1972 scm_remember_upto_here_1 (n
);
1973 return SCM_I_MAKINUM (count
);
1976 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1981 static const char scm_ilentab
[] = {
1982 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1986 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1988 "Return the number of bits necessary to represent @var{n}.\n"
1991 "(integer-length #b10101010)\n"
1993 "(integer-length 0)\n"
1995 "(integer-length #b1111)\n"
1998 #define FUNC_NAME s_scm_integer_length
2000 if (SCM_I_INUMP (n
))
2002 unsigned long int c
= 0;
2004 long int nn
= SCM_I_INUM (n
);
2010 l
= scm_ilentab
[15 & nn
];
2013 return SCM_I_MAKINUM (c
- 4 + l
);
2015 else if (SCM_BIGP (n
))
2017 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2018 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2019 1 too big, so check for that and adjust. */
2020 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2021 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2022 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2023 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2025 scm_remember_upto_here_1 (n
);
2026 return SCM_I_MAKINUM (size
);
2029 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2033 /*** NUMBERS -> STRINGS ***/
2034 #define SCM_MAX_DBL_PREC 60
2035 #define SCM_MAX_DBL_RADIX 36
2037 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2038 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2039 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2042 void init_dblprec(int *prec
, int radix
) {
2043 /* determine floating point precision by adding successively
2044 smaller increments to 1.0 until it is considered == 1.0 */
2045 double f
= ((double)1.0)/radix
;
2046 double fsum
= 1.0 + f
;
2051 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2063 void init_fx_radix(double *fx_list
, int radix
)
2065 /* initialize a per-radix list of tolerances. When added
2066 to a number < 1.0, we can determine if we should raund
2067 up and quit converting a number to a string. */
2071 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2072 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2075 /* use this array as a way to generate a single digit */
2076 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2079 idbl2str (double f
, char *a
, int radix
)
2081 int efmt
, dpt
, d
, i
, wp
;
2083 #ifdef DBL_MIN_10_EXP
2086 #endif /* DBL_MIN_10_EXP */
2091 radix
> SCM_MAX_DBL_RADIX
)
2093 /* revert to existing behavior */
2097 wp
= scm_dblprec
[radix
-2];
2098 fx
= fx_per_radix
[radix
-2];
2102 #ifdef HAVE_COPYSIGN
2103 double sgn
= copysign (1.0, f
);
2108 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2114 strcpy (a
, "-inf.0");
2116 strcpy (a
, "+inf.0");
2119 else if (xisnan (f
))
2121 strcpy (a
, "+nan.0");
2131 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2132 make-uniform-vector, from causing infinite loops. */
2133 /* just do the checking...if it passes, we do the conversion for our
2134 radix again below */
2141 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2149 while (f_cpy
> 10.0)
2152 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2173 if (f
+ fx
[wp
] >= radix
)
2180 /* adding 9999 makes this equivalent to abs(x) % 3 */
2181 dpt
= (exp
+ 9999) % 3;
2185 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2207 a
[ch
++] = number_chars
[d
];
2210 if (f
+ fx
[wp
] >= 1.0)
2212 a
[ch
- 1] = number_chars
[d
+1];
2224 if ((dpt
> 4) && (exp
> 6))
2226 d
= (a
[0] == '-' ? 2 : 1);
2227 for (i
= ch
++; i
> d
; i
--)
2240 if (a
[ch
- 1] == '.')
2241 a
[ch
++] = '0'; /* trailing zero */
2250 for (i
= radix
; i
<= exp
; i
*= radix
);
2251 for (i
/= radix
; i
; i
/= radix
)
2253 a
[ch
++] = number_chars
[exp
/ i
];
2262 icmplx2str (double real
, double imag
, char *str
, int radix
)
2266 i
= idbl2str (real
, str
, radix
);
2269 /* Don't output a '+' for negative numbers or for Inf and
2270 NaN. They will provide their own sign. */
2271 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2273 i
+= idbl2str (imag
, &str
[i
], radix
);
2280 iflo2str (SCM flt
, char *str
, int radix
)
2283 if (SCM_REALP (flt
))
2284 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2286 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2291 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2292 characters in the result.
2294 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2296 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2301 return scm_iuint2str (-num
, rad
, p
) + 1;
2304 return scm_iuint2str (num
, rad
, p
);
2307 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2308 characters in the result.
2310 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2312 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2316 scm_t_uintmax n
= num
;
2318 for (n
/= rad
; n
> 0; n
/= rad
)
2328 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2333 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2335 "Return a string holding the external representation of the\n"
2336 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2337 "inexact, a radix of 10 will be used.")
2338 #define FUNC_NAME s_scm_number_to_string
2342 if (SCM_UNBNDP (radix
))
2345 base
= scm_to_signed_integer (radix
, 2, 36);
2347 if (SCM_I_INUMP (n
))
2349 char num_buf
[SCM_INTBUFLEN
];
2350 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2351 return scm_from_locale_stringn (num_buf
, length
);
2353 else if (SCM_BIGP (n
))
2355 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2356 scm_remember_upto_here_1 (n
);
2357 return scm_take_locale_string (str
);
2359 else if (SCM_FRACTIONP (n
))
2361 scm_i_fraction_reduce (n
);
2362 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2363 scm_from_locale_string ("/"),
2364 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2366 else if (SCM_INEXACTP (n
))
2368 char num_buf
[FLOBUFLEN
];
2369 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2372 SCM_WRONG_TYPE_ARG (1, n
);
2377 /* These print routines used to be stubbed here so that scm_repl.c
2378 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2381 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2383 char num_buf
[FLOBUFLEN
];
2384 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2389 scm_i_print_double (double val
, SCM port
)
2391 char num_buf
[FLOBUFLEN
];
2392 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2396 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2399 char num_buf
[FLOBUFLEN
];
2400 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2405 scm_i_print_complex (double real
, double imag
, SCM port
)
2407 char num_buf
[FLOBUFLEN
];
2408 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2412 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2415 scm_i_fraction_reduce (sexp
);
2416 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2417 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2418 scm_remember_upto_here_1 (str
);
2423 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2425 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2426 scm_remember_upto_here_1 (exp
);
2427 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2431 /*** END nums->strs ***/
2434 /*** STRINGS -> NUMBERS ***/
2436 /* The following functions implement the conversion from strings to numbers.
2437 * The implementation somehow follows the grammar for numbers as it is given
2438 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2439 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2440 * points should be noted about the implementation:
2441 * * Each function keeps a local index variable 'idx' that points at the
2442 * current position within the parsed string. The global index is only
2443 * updated if the function could parse the corresponding syntactic unit
2445 * * Similarly, the functions keep track of indicators of inexactness ('#',
2446 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2447 * global exactness information is only updated after each part has been
2448 * successfully parsed.
2449 * * Sequences of digits are parsed into temporary variables holding fixnums.
2450 * Only if these fixnums would overflow, the result variables are updated
2451 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2452 * the temporary variables holding the fixnums are cleared, and the process
2453 * starts over again. If for example fixnums were able to store five decimal
2454 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2455 * and the result was computed as 12345 * 100000 + 67890. In other words,
2456 * only every five digits two bignum operations were performed.
2459 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2461 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2463 /* In non ASCII-style encodings the following macro might not work. */
2464 #define XDIGIT2UINT(d) \
2465 (isdigit ((int) (unsigned char) d) \
2467 : tolower ((int) (unsigned char) d) - 'a' + 10)
2470 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2471 unsigned int radix
, enum t_exactness
*p_exactness
)
2473 unsigned int idx
= *p_idx
;
2474 unsigned int hash_seen
= 0;
2475 scm_t_bits shift
= 1;
2477 unsigned int digit_value
;
2485 if (!isxdigit ((int) (unsigned char) c
))
2487 digit_value
= XDIGIT2UINT (c
);
2488 if (digit_value
>= radix
)
2492 result
= SCM_I_MAKINUM (digit_value
);
2496 if (isxdigit ((int) (unsigned char) c
))
2500 digit_value
= XDIGIT2UINT (c
);
2501 if (digit_value
>= radix
)
2513 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2515 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2517 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2524 shift
= shift
* radix
;
2525 add
= add
* radix
+ digit_value
;
2530 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2532 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2536 *p_exactness
= INEXACT
;
2542 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2543 * covers the parts of the rules that start at a potential point. The value
2544 * of the digits up to the point have been parsed by the caller and are given
2545 * in variable result. The content of *p_exactness indicates, whether a hash
2546 * has already been seen in the digits before the point.
2549 /* In non ASCII-style encodings the following macro might not work. */
2550 #define DIGIT2UINT(d) ((d) - '0')
2553 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2554 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2556 unsigned int idx
= *p_idx
;
2557 enum t_exactness x
= *p_exactness
;
2562 if (mem
[idx
] == '.')
2564 scm_t_bits shift
= 1;
2566 unsigned int digit_value
;
2567 SCM big_shift
= SCM_I_MAKINUM (1);
2573 if (isdigit ((int) (unsigned char) c
))
2578 digit_value
= DIGIT2UINT (c
);
2589 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2591 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2592 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2594 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2602 add
= add
* 10 + digit_value
;
2608 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2609 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2610 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2613 result
= scm_divide (result
, big_shift
);
2615 /* We've seen a decimal point, thus the value is implicitly inexact. */
2627 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2654 if (!isdigit ((int) (unsigned char) c
))
2658 exponent
= DIGIT2UINT (c
);
2662 if (isdigit ((int) (unsigned char) c
))
2665 if (exponent
<= SCM_MAXEXP
)
2666 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2672 if (exponent
> SCM_MAXEXP
)
2674 size_t exp_len
= idx
- start
;
2675 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2676 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2677 scm_out_of_range ("string->number", exp_num
);
2680 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2682 result
= scm_product (result
, e
);
2684 result
= scm_divide2real (result
, e
);
2686 /* We've seen an exponent, thus the value is implicitly inexact. */
2704 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2707 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2708 unsigned int radix
, enum t_exactness
*p_exactness
)
2710 unsigned int idx
= *p_idx
;
2716 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2722 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2724 enum t_exactness x
= EXACT
;
2726 /* Cobble up the fractional part. We might want to set the
2727 NaN's mantissa from it. */
2729 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2734 if (mem
[idx
] == '.')
2738 else if (idx
+ 1 == len
)
2740 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2743 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2744 p_idx
, p_exactness
);
2748 enum t_exactness x
= EXACT
;
2751 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2752 if (scm_is_false (uinteger
))
2757 else if (mem
[idx
] == '/')
2763 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2764 if (scm_is_false (divisor
))
2767 /* both are int/big here, I assume */
2768 result
= scm_i_make_ratio (uinteger
, divisor
);
2770 else if (radix
== 10)
2772 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2773 if (scm_is_false (result
))
2784 /* When returning an inexact zero, make sure it is represented as a
2785 floating point value so that we can change its sign.
2787 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2788 result
= scm_from_double (0.0);
2794 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2797 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2798 unsigned int radix
, enum t_exactness
*p_exactness
)
2822 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2823 if (scm_is_false (ureal
))
2825 /* input must be either +i or -i */
2830 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2836 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2843 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2844 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2853 /* either +<ureal>i or -<ureal>i */
2860 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2863 /* polar input: <real>@<real>. */
2888 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2889 if (scm_is_false (angle
))
2894 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2895 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2897 result
= scm_make_polar (ureal
, angle
);
2902 /* expecting input matching <real>[+-]<ureal>?i */
2909 int sign
= (c
== '+') ? 1 : -1;
2910 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2912 if (scm_is_false (imag
))
2913 imag
= SCM_I_MAKINUM (sign
);
2914 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2915 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2919 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2926 return scm_make_rectangular (ureal
, imag
);
2935 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2937 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2940 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2942 unsigned int idx
= 0;
2943 unsigned int radix
= NO_RADIX
;
2944 enum t_exactness forced_x
= NO_EXACTNESS
;
2945 enum t_exactness implicit_x
= EXACT
;
2948 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2949 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2951 switch (mem
[idx
+ 1])
2954 if (radix
!= NO_RADIX
)
2959 if (radix
!= NO_RADIX
)
2964 if (forced_x
!= NO_EXACTNESS
)
2969 if (forced_x
!= NO_EXACTNESS
)
2974 if (radix
!= NO_RADIX
)
2979 if (radix
!= NO_RADIX
)
2989 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2990 if (radix
== NO_RADIX
)
2991 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2993 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2995 if (scm_is_false (result
))
3001 if (SCM_INEXACTP (result
))
3002 return scm_inexact_to_exact (result
);
3006 if (SCM_INEXACTP (result
))
3009 return scm_exact_to_inexact (result
);
3012 if (implicit_x
== INEXACT
)
3014 if (SCM_INEXACTP (result
))
3017 return scm_exact_to_inexact (result
);
3025 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3026 (SCM string
, SCM radix
),
3027 "Return a number of the maximally precise representation\n"
3028 "expressed by the given @var{string}. @var{radix} must be an\n"
3029 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3030 "is a default radix that may be overridden by an explicit radix\n"
3031 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3032 "supplied, then the default radix is 10. If string is not a\n"
3033 "syntactically valid notation for a number, then\n"
3034 "@code{string->number} returns @code{#f}.")
3035 #define FUNC_NAME s_scm_string_to_number
3039 SCM_VALIDATE_STRING (1, string
);
3041 if (SCM_UNBNDP (radix
))
3044 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3046 answer
= scm_i_mem2number (scm_i_string_chars (string
),
3047 scm_i_string_length (string
),
3049 scm_remember_upto_here_1 (string
);
3055 /*** END strs->nums ***/
3059 scm_bigequal (SCM x
, SCM y
)
3061 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3062 scm_remember_upto_here_2 (x
, y
);
3063 return scm_from_bool (0 == result
);
3067 scm_real_equalp (SCM x
, SCM y
)
3069 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3073 scm_complex_equalp (SCM x
, SCM y
)
3075 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3076 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3080 scm_i_fraction_equalp (SCM x
, SCM y
)
3082 scm_i_fraction_reduce (x
);
3083 scm_i_fraction_reduce (y
);
3084 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3085 SCM_FRACTION_NUMERATOR (y
)))
3086 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3087 SCM_FRACTION_DENOMINATOR (y
))))
3094 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3096 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3098 #define FUNC_NAME s_scm_number_p
3100 return scm_from_bool (SCM_NUMBERP (x
));
3104 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3106 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3107 "otherwise. Note that the sets of real, rational and integer\n"
3108 "values form subsets of the set of complex numbers, i. e. the\n"
3109 "predicate will also be fulfilled if @var{x} is a real,\n"
3110 "rational or integer number.")
3111 #define FUNC_NAME s_scm_complex_p
3113 /* all numbers are complex. */
3114 return scm_number_p (x
);
3118 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3120 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3121 "otherwise. Note that the set of integer values forms a subset of\n"
3122 "the set of real numbers, i. e. the predicate will also be\n"
3123 "fulfilled if @var{x} is an integer number.")
3124 #define FUNC_NAME s_scm_real_p
3126 /* we can't represent irrational numbers. */
3127 return scm_rational_p (x
);
3131 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3133 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3134 "otherwise. Note that the set of integer values forms a subset of\n"
3135 "the set of rational numbers, i. e. the predicate will also be\n"
3136 "fulfilled if @var{x} is an integer number.")
3137 #define FUNC_NAME s_scm_rational_p
3139 if (SCM_I_INUMP (x
))
3141 else if (SCM_IMP (x
))
3143 else if (SCM_BIGP (x
))
3145 else if (SCM_FRACTIONP (x
))
3147 else if (SCM_REALP (x
))
3148 /* due to their limited precision, all floating point numbers are
3149 rational as well. */
3156 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3158 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3160 #define FUNC_NAME s_scm_integer_p
3163 if (SCM_I_INUMP (x
))
3169 if (!SCM_INEXACTP (x
))
3171 if (SCM_COMPLEXP (x
))
3173 r
= SCM_REAL_VALUE (x
);
3174 /* +/-inf passes r==floor(r), making those #t */
3182 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3184 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3186 #define FUNC_NAME s_scm_inexact_p
3188 if (SCM_INEXACTP (x
))
3190 if (SCM_NUMBERP (x
))
3192 SCM_WRONG_TYPE_ARG (1, x
);
3197 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3198 /* "Return @code{#t} if all parameters are numerically equal." */
3200 scm_num_eq_p (SCM x
, SCM y
)
3203 if (SCM_I_INUMP (x
))
3205 long xx
= SCM_I_INUM (x
);
3206 if (SCM_I_INUMP (y
))
3208 long yy
= SCM_I_INUM (y
);
3209 return scm_from_bool (xx
== yy
);
3211 else if (SCM_BIGP (y
))
3213 else if (SCM_REALP (y
))
3215 /* On a 32-bit system an inum fits a double, we can cast the inum
3216 to a double and compare.
3218 But on a 64-bit system an inum is bigger than a double and
3219 casting it to a double (call that dxx) will round. dxx is at
3220 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3221 an integer and fits a long. So we cast yy to a long and
3222 compare with plain xx.
3224 An alternative (for any size system actually) would be to check
3225 yy is an integer (with floor) and is in range of an inum
3226 (compare against appropriate powers of 2) then test
3227 xx==(long)yy. It's just a matter of which casts/comparisons
3228 might be fastest or easiest for the cpu. */
3230 double yy
= SCM_REAL_VALUE (y
);
3231 return SCM_BOOL ((double) xx
== yy
3232 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3233 || xx
== (long) yy
));
3235 else if (SCM_COMPLEXP (y
))
3236 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3237 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3238 else if (SCM_FRACTIONP (y
))
3241 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3243 else if (SCM_BIGP (x
))
3245 if (SCM_I_INUMP (y
))
3247 else if (SCM_BIGP (y
))
3249 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3250 scm_remember_upto_here_2 (x
, y
);
3251 return scm_from_bool (0 == cmp
);
3253 else if (SCM_REALP (y
))
3256 if (xisnan (SCM_REAL_VALUE (y
)))
3258 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3259 scm_remember_upto_here_1 (x
);
3260 return scm_from_bool (0 == cmp
);
3262 else if (SCM_COMPLEXP (y
))
3265 if (0.0 != SCM_COMPLEX_IMAG (y
))
3267 if (xisnan (SCM_COMPLEX_REAL (y
)))
3269 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3270 scm_remember_upto_here_1 (x
);
3271 return scm_from_bool (0 == cmp
);
3273 else if (SCM_FRACTIONP (y
))
3276 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3278 else if (SCM_REALP (x
))
3280 double xx
= SCM_REAL_VALUE (x
);
3281 if (SCM_I_INUMP (y
))
3283 /* see comments with inum/real above */
3284 long yy
= SCM_I_INUM (y
);
3285 return SCM_BOOL (xx
== (double) yy
3286 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3287 || (long) xx
== yy
));
3289 else if (SCM_BIGP (y
))
3292 if (xisnan (SCM_REAL_VALUE (x
)))
3294 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3295 scm_remember_upto_here_1 (y
);
3296 return scm_from_bool (0 == cmp
);
3298 else if (SCM_REALP (y
))
3299 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3300 else if (SCM_COMPLEXP (y
))
3301 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3302 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3303 else if (SCM_FRACTIONP (y
))
3305 double xx
= SCM_REAL_VALUE (x
);
3309 return scm_from_bool (xx
< 0.0);
3310 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3314 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3316 else if (SCM_COMPLEXP (x
))
3318 if (SCM_I_INUMP (y
))
3319 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3320 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3321 else if (SCM_BIGP (y
))
3324 if (0.0 != SCM_COMPLEX_IMAG (x
))
3326 if (xisnan (SCM_COMPLEX_REAL (x
)))
3328 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3329 scm_remember_upto_here_1 (y
);
3330 return scm_from_bool (0 == cmp
);
3332 else if (SCM_REALP (y
))
3333 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3334 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3335 else if (SCM_COMPLEXP (y
))
3336 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3337 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3338 else if (SCM_FRACTIONP (y
))
3341 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3343 xx
= SCM_COMPLEX_REAL (x
);
3347 return scm_from_bool (xx
< 0.0);
3348 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3352 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3354 else if (SCM_FRACTIONP (x
))
3356 if (SCM_I_INUMP (y
))
3358 else if (SCM_BIGP (y
))
3360 else if (SCM_REALP (y
))
3362 double yy
= SCM_REAL_VALUE (y
);
3366 return scm_from_bool (0.0 < yy
);
3367 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3370 else if (SCM_COMPLEXP (y
))
3373 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3375 yy
= SCM_COMPLEX_REAL (y
);
3379 return scm_from_bool (0.0 < yy
);
3380 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3383 else if (SCM_FRACTIONP (y
))
3384 return scm_i_fraction_equalp (x
, y
);
3386 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3389 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3393 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3394 done are good for inums, but for bignums an answer can almost always be
3395 had by just examining a few high bits of the operands, as done by GMP in
3396 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3397 of the float exponent to take into account. */
3399 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3400 /* "Return @code{#t} if the list of parameters is monotonically\n"
3404 scm_less_p (SCM x
, SCM y
)
3407 if (SCM_I_INUMP (x
))
3409 long xx
= SCM_I_INUM (x
);
3410 if (SCM_I_INUMP (y
))
3412 long yy
= SCM_I_INUM (y
);
3413 return scm_from_bool (xx
< yy
);
3415 else if (SCM_BIGP (y
))
3417 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3418 scm_remember_upto_here_1 (y
);
3419 return scm_from_bool (sgn
> 0);
3421 else if (SCM_REALP (y
))
3422 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3423 else if (SCM_FRACTIONP (y
))
3425 /* "x < a/b" becomes "x*b < a" */
3427 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3428 y
= SCM_FRACTION_NUMERATOR (y
);
3432 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3434 else if (SCM_BIGP (x
))
3436 if (SCM_I_INUMP (y
))
3438 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3439 scm_remember_upto_here_1 (x
);
3440 return scm_from_bool (sgn
< 0);
3442 else if (SCM_BIGP (y
))
3444 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3445 scm_remember_upto_here_2 (x
, y
);
3446 return scm_from_bool (cmp
< 0);
3448 else if (SCM_REALP (y
))
3451 if (xisnan (SCM_REAL_VALUE (y
)))
3453 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3454 scm_remember_upto_here_1 (x
);
3455 return scm_from_bool (cmp
< 0);
3457 else if (SCM_FRACTIONP (y
))
3460 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3462 else if (SCM_REALP (x
))
3464 if (SCM_I_INUMP (y
))
3465 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3466 else if (SCM_BIGP (y
))
3469 if (xisnan (SCM_REAL_VALUE (x
)))
3471 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3472 scm_remember_upto_here_1 (y
);
3473 return scm_from_bool (cmp
> 0);
3475 else if (SCM_REALP (y
))
3476 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3477 else if (SCM_FRACTIONP (y
))
3479 double xx
= SCM_REAL_VALUE (x
);
3483 return scm_from_bool (xx
< 0.0);
3484 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3488 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3490 else if (SCM_FRACTIONP (x
))
3492 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3494 /* "a/b < y" becomes "a < y*b" */
3495 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3496 x
= SCM_FRACTION_NUMERATOR (x
);
3499 else if (SCM_REALP (y
))
3501 double yy
= SCM_REAL_VALUE (y
);
3505 return scm_from_bool (0.0 < yy
);
3506 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3509 else if (SCM_FRACTIONP (y
))
3511 /* "a/b < c/d" becomes "a*d < c*b" */
3512 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3513 SCM_FRACTION_DENOMINATOR (y
));
3514 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3515 SCM_FRACTION_DENOMINATOR (x
));
3521 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3524 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3528 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3529 /* "Return @code{#t} if the list of parameters is monotonically\n"
3532 #define FUNC_NAME s_scm_gr_p
3534 scm_gr_p (SCM x
, SCM y
)
3536 if (!SCM_NUMBERP (x
))
3537 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3538 else if (!SCM_NUMBERP (y
))
3539 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3541 return scm_less_p (y
, x
);
3546 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3547 /* "Return @code{#t} if the list of parameters is monotonically\n"
3550 #define FUNC_NAME s_scm_leq_p
3552 scm_leq_p (SCM x
, SCM y
)
3554 if (!SCM_NUMBERP (x
))
3555 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3556 else if (!SCM_NUMBERP (y
))
3557 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3558 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3561 return scm_not (scm_less_p (y
, x
));
3566 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3567 /* "Return @code{#t} if the list of parameters is monotonically\n"
3570 #define FUNC_NAME s_scm_geq_p
3572 scm_geq_p (SCM x
, SCM y
)
3574 if (!SCM_NUMBERP (x
))
3575 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3576 else if (!SCM_NUMBERP (y
))
3577 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3578 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3581 return scm_not (scm_less_p (x
, y
));
3586 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3587 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3593 if (SCM_I_INUMP (z
))
3594 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3595 else if (SCM_BIGP (z
))
3597 else if (SCM_REALP (z
))
3598 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3599 else if (SCM_COMPLEXP (z
))
3600 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3601 && SCM_COMPLEX_IMAG (z
) == 0.0);
3602 else if (SCM_FRACTIONP (z
))
3605 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3609 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3610 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3614 scm_positive_p (SCM x
)
3616 if (SCM_I_INUMP (x
))
3617 return scm_from_bool (SCM_I_INUM (x
) > 0);
3618 else if (SCM_BIGP (x
))
3620 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3621 scm_remember_upto_here_1 (x
);
3622 return scm_from_bool (sgn
> 0);
3624 else if (SCM_REALP (x
))
3625 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3626 else if (SCM_FRACTIONP (x
))
3627 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3629 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3633 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3634 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3638 scm_negative_p (SCM x
)
3640 if (SCM_I_INUMP (x
))
3641 return scm_from_bool (SCM_I_INUM (x
) < 0);
3642 else if (SCM_BIGP (x
))
3644 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3645 scm_remember_upto_here_1 (x
);
3646 return scm_from_bool (sgn
< 0);
3648 else if (SCM_REALP (x
))
3649 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3650 else if (SCM_FRACTIONP (x
))
3651 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3653 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3657 /* scm_min and scm_max return an inexact when either argument is inexact, as
3658 required by r5rs. On that basis, for exact/inexact combinations the
3659 exact is converted to inexact to compare and possibly return. This is
3660 unlike scm_less_p above which takes some trouble to preserve all bits in
3661 its test, such trouble is not required for min and max. */
3663 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3664 /* "Return the maximum of all parameter values."
3667 scm_max (SCM x
, SCM y
)
3672 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3673 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3676 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3679 if (SCM_I_INUMP (x
))
3681 long xx
= SCM_I_INUM (x
);
3682 if (SCM_I_INUMP (y
))
3684 long yy
= SCM_I_INUM (y
);
3685 return (xx
< yy
) ? y
: x
;
3687 else if (SCM_BIGP (y
))
3689 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3690 scm_remember_upto_here_1 (y
);
3691 return (sgn
< 0) ? x
: y
;
3693 else if (SCM_REALP (y
))
3696 /* if y==NaN then ">" is false and we return NaN */
3697 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3699 else if (SCM_FRACTIONP (y
))
3702 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3705 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3707 else if (SCM_BIGP (x
))
3709 if (SCM_I_INUMP (y
))
3711 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3712 scm_remember_upto_here_1 (x
);
3713 return (sgn
< 0) ? y
: x
;
3715 else if (SCM_BIGP (y
))
3717 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3718 scm_remember_upto_here_2 (x
, y
);
3719 return (cmp
> 0) ? x
: y
;
3721 else if (SCM_REALP (y
))
3723 /* if y==NaN then xx>yy is false, so we return the NaN y */
3726 xx
= scm_i_big2dbl (x
);
3727 yy
= SCM_REAL_VALUE (y
);
3728 return (xx
> yy
? scm_from_double (xx
) : y
);
3730 else if (SCM_FRACTIONP (y
))
3735 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3737 else if (SCM_REALP (x
))
3739 if (SCM_I_INUMP (y
))
3741 double z
= SCM_I_INUM (y
);
3742 /* if x==NaN then "<" is false and we return NaN */
3743 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3745 else if (SCM_BIGP (y
))
3750 else if (SCM_REALP (y
))
3752 /* if x==NaN then our explicit check means we return NaN
3753 if y==NaN then ">" is false and we return NaN
3754 calling isnan is unavoidable, since it's the only way to know
3755 which of x or y causes any compares to be false */
3756 double xx
= SCM_REAL_VALUE (x
);
3757 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3759 else if (SCM_FRACTIONP (y
))
3761 double yy
= scm_i_fraction2double (y
);
3762 double xx
= SCM_REAL_VALUE (x
);
3763 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3766 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3768 else if (SCM_FRACTIONP (x
))
3770 if (SCM_I_INUMP (y
))
3774 else if (SCM_BIGP (y
))
3778 else if (SCM_REALP (y
))
3780 double xx
= scm_i_fraction2double (x
);
3781 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3783 else if (SCM_FRACTIONP (y
))
3788 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3791 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3795 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3796 /* "Return the minium of all parameter values."
3799 scm_min (SCM x
, SCM y
)
3804 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3805 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3808 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3811 if (SCM_I_INUMP (x
))
3813 long xx
= SCM_I_INUM (x
);
3814 if (SCM_I_INUMP (y
))
3816 long yy
= SCM_I_INUM (y
);
3817 return (xx
< yy
) ? x
: y
;
3819 else if (SCM_BIGP (y
))
3821 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3822 scm_remember_upto_here_1 (y
);
3823 return (sgn
< 0) ? y
: x
;
3825 else if (SCM_REALP (y
))
3828 /* if y==NaN then "<" is false and we return NaN */
3829 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3831 else if (SCM_FRACTIONP (y
))
3834 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3837 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3839 else if (SCM_BIGP (x
))
3841 if (SCM_I_INUMP (y
))
3843 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3844 scm_remember_upto_here_1 (x
);
3845 return (sgn
< 0) ? x
: y
;
3847 else if (SCM_BIGP (y
))
3849 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3850 scm_remember_upto_here_2 (x
, y
);
3851 return (cmp
> 0) ? y
: x
;
3853 else if (SCM_REALP (y
))
3855 /* if y==NaN then xx<yy is false, so we return the NaN y */
3858 xx
= scm_i_big2dbl (x
);
3859 yy
= SCM_REAL_VALUE (y
);
3860 return (xx
< yy
? scm_from_double (xx
) : y
);
3862 else if (SCM_FRACTIONP (y
))
3867 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3869 else if (SCM_REALP (x
))
3871 if (SCM_I_INUMP (y
))
3873 double z
= SCM_I_INUM (y
);
3874 /* if x==NaN then "<" is false and we return NaN */
3875 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3877 else if (SCM_BIGP (y
))
3882 else if (SCM_REALP (y
))
3884 /* if x==NaN then our explicit check means we return NaN
3885 if y==NaN then "<" is false and we return NaN
3886 calling isnan is unavoidable, since it's the only way to know
3887 which of x or y causes any compares to be false */
3888 double xx
= SCM_REAL_VALUE (x
);
3889 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3891 else if (SCM_FRACTIONP (y
))
3893 double yy
= scm_i_fraction2double (y
);
3894 double xx
= SCM_REAL_VALUE (x
);
3895 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3898 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3900 else if (SCM_FRACTIONP (x
))
3902 if (SCM_I_INUMP (y
))
3906 else if (SCM_BIGP (y
))
3910 else if (SCM_REALP (y
))
3912 double xx
= scm_i_fraction2double (x
);
3913 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3915 else if (SCM_FRACTIONP (y
))
3920 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3923 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3927 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3928 /* "Return the sum of all parameter values. Return 0 if called without\n"
3932 scm_sum (SCM x
, SCM y
)
3936 if (SCM_NUMBERP (x
)) return x
;
3937 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3938 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3941 if (SCM_I_INUMP (x
))
3943 if (SCM_I_INUMP (y
))
3945 long xx
= SCM_I_INUM (x
);
3946 long yy
= SCM_I_INUM (y
);
3947 long int z
= xx
+ yy
;
3948 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3950 else if (SCM_BIGP (y
))
3955 else if (SCM_REALP (y
))
3957 long int xx
= SCM_I_INUM (x
);
3958 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3960 else if (SCM_COMPLEXP (y
))
3962 long int xx
= SCM_I_INUM (x
);
3963 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3964 SCM_COMPLEX_IMAG (y
));
3966 else if (SCM_FRACTIONP (y
))
3967 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3968 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3969 SCM_FRACTION_DENOMINATOR (y
));
3971 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3972 } else if (SCM_BIGP (x
))
3974 if (SCM_I_INUMP (y
))
3979 inum
= SCM_I_INUM (y
);
3982 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3985 SCM result
= scm_i_mkbig ();
3986 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3987 scm_remember_upto_here_1 (x
);
3988 /* we know the result will have to be a bignum */
3991 return scm_i_normbig (result
);
3995 SCM result
= scm_i_mkbig ();
3996 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3997 scm_remember_upto_here_1 (x
);
3998 /* we know the result will have to be a bignum */
4001 return scm_i_normbig (result
);
4004 else if (SCM_BIGP (y
))
4006 SCM result
= scm_i_mkbig ();
4007 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4008 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4009 mpz_add (SCM_I_BIG_MPZ (result
),
4012 scm_remember_upto_here_2 (x
, y
);
4013 /* we know the result will have to be a bignum */
4016 return scm_i_normbig (result
);
4018 else if (SCM_REALP (y
))
4020 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4021 scm_remember_upto_here_1 (x
);
4022 return scm_from_double (result
);
4024 else if (SCM_COMPLEXP (y
))
4026 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4027 + SCM_COMPLEX_REAL (y
));
4028 scm_remember_upto_here_1 (x
);
4029 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4031 else if (SCM_FRACTIONP (y
))
4032 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4033 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4034 SCM_FRACTION_DENOMINATOR (y
));
4036 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4038 else if (SCM_REALP (x
))
4040 if (SCM_I_INUMP (y
))
4041 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4042 else if (SCM_BIGP (y
))
4044 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4045 scm_remember_upto_here_1 (y
);
4046 return scm_from_double (result
);
4048 else if (SCM_REALP (y
))
4049 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4050 else if (SCM_COMPLEXP (y
))
4051 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4052 SCM_COMPLEX_IMAG (y
));
4053 else if (SCM_FRACTIONP (y
))
4054 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4056 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4058 else if (SCM_COMPLEXP (x
))
4060 if (SCM_I_INUMP (y
))
4061 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4062 SCM_COMPLEX_IMAG (x
));
4063 else if (SCM_BIGP (y
))
4065 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4066 + SCM_COMPLEX_REAL (x
));
4067 scm_remember_upto_here_1 (y
);
4068 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4070 else if (SCM_REALP (y
))
4071 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4072 SCM_COMPLEX_IMAG (x
));
4073 else if (SCM_COMPLEXP (y
))
4074 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4075 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4076 else if (SCM_FRACTIONP (y
))
4077 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4078 SCM_COMPLEX_IMAG (x
));
4080 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4082 else if (SCM_FRACTIONP (x
))
4084 if (SCM_I_INUMP (y
))
4085 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4086 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4087 SCM_FRACTION_DENOMINATOR (x
));
4088 else if (SCM_BIGP (y
))
4089 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4090 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4091 SCM_FRACTION_DENOMINATOR (x
));
4092 else if (SCM_REALP (y
))
4093 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4094 else if (SCM_COMPLEXP (y
))
4095 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4096 SCM_COMPLEX_IMAG (y
));
4097 else if (SCM_FRACTIONP (y
))
4098 /* a/b + c/d = (ad + bc) / bd */
4099 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4100 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4101 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4103 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4106 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4110 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4112 "Return @math{@var{x}+1}.")
4113 #define FUNC_NAME s_scm_oneplus
4115 return scm_sum (x
, SCM_I_MAKINUM (1));
4120 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4121 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4122 * the sum of all but the first argument are subtracted from the first
4124 #define FUNC_NAME s_difference
4126 scm_difference (SCM x
, SCM y
)
4131 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4133 if (SCM_I_INUMP (x
))
4135 long xx
= -SCM_I_INUM (x
);
4136 if (SCM_FIXABLE (xx
))
4137 return SCM_I_MAKINUM (xx
);
4139 return scm_i_long2big (xx
);
4141 else if (SCM_BIGP (x
))
4142 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4143 bignum, but negating that gives a fixnum. */
4144 return scm_i_normbig (scm_i_clonebig (x
, 0));
4145 else if (SCM_REALP (x
))
4146 return scm_from_double (-SCM_REAL_VALUE (x
));
4147 else if (SCM_COMPLEXP (x
))
4148 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4149 -SCM_COMPLEX_IMAG (x
));
4150 else if (SCM_FRACTIONP (x
))
4151 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4152 SCM_FRACTION_DENOMINATOR (x
));
4154 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4157 if (SCM_I_INUMP (x
))
4159 if (SCM_I_INUMP (y
))
4161 long int xx
= SCM_I_INUM (x
);
4162 long int yy
= SCM_I_INUM (y
);
4163 long int z
= xx
- yy
;
4164 if (SCM_FIXABLE (z
))
4165 return SCM_I_MAKINUM (z
);
4167 return scm_i_long2big (z
);
4169 else if (SCM_BIGP (y
))
4171 /* inum-x - big-y */
4172 long xx
= SCM_I_INUM (x
);
4175 return scm_i_clonebig (y
, 0);
4178 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4179 SCM result
= scm_i_mkbig ();
4182 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4185 /* x - y == -(y + -x) */
4186 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4187 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4189 scm_remember_upto_here_1 (y
);
4191 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4192 /* we know the result will have to be a bignum */
4195 return scm_i_normbig (result
);
4198 else if (SCM_REALP (y
))
4200 long int xx
= SCM_I_INUM (x
);
4201 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4203 else if (SCM_COMPLEXP (y
))
4205 long int xx
= SCM_I_INUM (x
);
4206 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4207 - SCM_COMPLEX_IMAG (y
));
4209 else if (SCM_FRACTIONP (y
))
4210 /* a - b/c = (ac - b) / c */
4211 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4212 SCM_FRACTION_NUMERATOR (y
)),
4213 SCM_FRACTION_DENOMINATOR (y
));
4215 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4217 else if (SCM_BIGP (x
))
4219 if (SCM_I_INUMP (y
))
4221 /* big-x - inum-y */
4222 long yy
= SCM_I_INUM (y
);
4223 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4225 scm_remember_upto_here_1 (x
);
4227 return (SCM_FIXABLE (-yy
) ?
4228 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4231 SCM result
= scm_i_mkbig ();
4234 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4236 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4237 scm_remember_upto_here_1 (x
);
4239 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4240 /* we know the result will have to be a bignum */
4243 return scm_i_normbig (result
);
4246 else if (SCM_BIGP (y
))
4248 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4249 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4250 SCM result
= scm_i_mkbig ();
4251 mpz_sub (SCM_I_BIG_MPZ (result
),
4254 scm_remember_upto_here_2 (x
, y
);
4255 /* we know the result will have to be a bignum */
4256 if ((sgn_x
== 1) && (sgn_y
== -1))
4258 if ((sgn_x
== -1) && (sgn_y
== 1))
4260 return scm_i_normbig (result
);
4262 else if (SCM_REALP (y
))
4264 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4265 scm_remember_upto_here_1 (x
);
4266 return scm_from_double (result
);
4268 else if (SCM_COMPLEXP (y
))
4270 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4271 - SCM_COMPLEX_REAL (y
));
4272 scm_remember_upto_here_1 (x
);
4273 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4275 else if (SCM_FRACTIONP (y
))
4276 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4277 SCM_FRACTION_NUMERATOR (y
)),
4278 SCM_FRACTION_DENOMINATOR (y
));
4279 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4281 else if (SCM_REALP (x
))
4283 if (SCM_I_INUMP (y
))
4284 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4285 else if (SCM_BIGP (y
))
4287 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4288 scm_remember_upto_here_1 (x
);
4289 return scm_from_double (result
);
4291 else if (SCM_REALP (y
))
4292 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4293 else if (SCM_COMPLEXP (y
))
4294 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4295 -SCM_COMPLEX_IMAG (y
));
4296 else if (SCM_FRACTIONP (y
))
4297 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4299 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4301 else if (SCM_COMPLEXP (x
))
4303 if (SCM_I_INUMP (y
))
4304 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4305 SCM_COMPLEX_IMAG (x
));
4306 else if (SCM_BIGP (y
))
4308 double real_part
= (SCM_COMPLEX_REAL (x
)
4309 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4310 scm_remember_upto_here_1 (x
);
4311 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4313 else if (SCM_REALP (y
))
4314 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4315 SCM_COMPLEX_IMAG (x
));
4316 else if (SCM_COMPLEXP (y
))
4317 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4318 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4319 else if (SCM_FRACTIONP (y
))
4320 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4321 SCM_COMPLEX_IMAG (x
));
4323 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4325 else if (SCM_FRACTIONP (x
))
4327 if (SCM_I_INUMP (y
))
4328 /* a/b - c = (a - cb) / b */
4329 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4330 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4331 SCM_FRACTION_DENOMINATOR (x
));
4332 else if (SCM_BIGP (y
))
4333 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4334 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4335 SCM_FRACTION_DENOMINATOR (x
));
4336 else if (SCM_REALP (y
))
4337 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4338 else if (SCM_COMPLEXP (y
))
4339 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4340 -SCM_COMPLEX_IMAG (y
));
4341 else if (SCM_FRACTIONP (y
))
4342 /* a/b - c/d = (ad - bc) / bd */
4343 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4344 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4345 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4347 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4350 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4355 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4357 "Return @math{@var{x}-1}.")
4358 #define FUNC_NAME s_scm_oneminus
4360 return scm_difference (x
, SCM_I_MAKINUM (1));
4365 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4366 /* "Return the product of all arguments. If called without arguments,\n"
4370 scm_product (SCM x
, SCM y
)
4375 return SCM_I_MAKINUM (1L);
4376 else if (SCM_NUMBERP (x
))
4379 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4382 if (SCM_I_INUMP (x
))
4387 xx
= SCM_I_INUM (x
);
4391 case 0: return x
; break;
4392 case 1: return y
; break;
4395 if (SCM_I_INUMP (y
))
4397 long yy
= SCM_I_INUM (y
);
4399 SCM k
= SCM_I_MAKINUM (kk
);
4400 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4404 SCM result
= scm_i_long2big (xx
);
4405 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4406 return scm_i_normbig (result
);
4409 else if (SCM_BIGP (y
))
4411 SCM result
= scm_i_mkbig ();
4412 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4413 scm_remember_upto_here_1 (y
);
4416 else if (SCM_REALP (y
))
4417 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4418 else if (SCM_COMPLEXP (y
))
4419 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4420 xx
* SCM_COMPLEX_IMAG (y
));
4421 else if (SCM_FRACTIONP (y
))
4422 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4423 SCM_FRACTION_DENOMINATOR (y
));
4425 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4427 else if (SCM_BIGP (x
))
4429 if (SCM_I_INUMP (y
))
4434 else if (SCM_BIGP (y
))
4436 SCM result
= scm_i_mkbig ();
4437 mpz_mul (SCM_I_BIG_MPZ (result
),
4440 scm_remember_upto_here_2 (x
, y
);
4443 else if (SCM_REALP (y
))
4445 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4446 scm_remember_upto_here_1 (x
);
4447 return scm_from_double (result
);
4449 else if (SCM_COMPLEXP (y
))
4451 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4452 scm_remember_upto_here_1 (x
);
4453 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4454 z
* SCM_COMPLEX_IMAG (y
));
4456 else if (SCM_FRACTIONP (y
))
4457 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4458 SCM_FRACTION_DENOMINATOR (y
));
4460 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4462 else if (SCM_REALP (x
))
4464 if (SCM_I_INUMP (y
))
4465 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4466 else if (SCM_BIGP (y
))
4468 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4469 scm_remember_upto_here_1 (y
);
4470 return scm_from_double (result
);
4472 else if (SCM_REALP (y
))
4473 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4474 else if (SCM_COMPLEXP (y
))
4475 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4476 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4477 else if (SCM_FRACTIONP (y
))
4478 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4480 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4482 else if (SCM_COMPLEXP (x
))
4484 if (SCM_I_INUMP (y
))
4485 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4486 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4487 else if (SCM_BIGP (y
))
4489 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4490 scm_remember_upto_here_1 (y
);
4491 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4492 z
* SCM_COMPLEX_IMAG (x
));
4494 else if (SCM_REALP (y
))
4495 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4496 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4497 else if (SCM_COMPLEXP (y
))
4499 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4500 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4501 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4502 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4504 else if (SCM_FRACTIONP (y
))
4506 double yy
= scm_i_fraction2double (y
);
4507 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4508 yy
* SCM_COMPLEX_IMAG (x
));
4511 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4513 else if (SCM_FRACTIONP (x
))
4515 if (SCM_I_INUMP (y
))
4516 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4517 SCM_FRACTION_DENOMINATOR (x
));
4518 else if (SCM_BIGP (y
))
4519 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4520 SCM_FRACTION_DENOMINATOR (x
));
4521 else if (SCM_REALP (y
))
4522 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4523 else if (SCM_COMPLEXP (y
))
4525 double xx
= scm_i_fraction2double (x
);
4526 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4527 xx
* SCM_COMPLEX_IMAG (y
));
4529 else if (SCM_FRACTIONP (y
))
4530 /* a/b * c/d = ac / bd */
4531 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4532 SCM_FRACTION_NUMERATOR (y
)),
4533 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4534 SCM_FRACTION_DENOMINATOR (y
)));
4536 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4539 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4542 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4543 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4544 #define ALLOW_DIVIDE_BY_ZERO
4545 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4548 /* The code below for complex division is adapted from the GNU
4549 libstdc++, which adapted it from f2c's libF77, and is subject to
4552 /****************************************************************
4553 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4555 Permission to use, copy, modify, and distribute this software
4556 and its documentation for any purpose and without fee is hereby
4557 granted, provided that the above copyright notice appear in all
4558 copies and that both that the copyright notice and this
4559 permission notice and warranty disclaimer appear in supporting
4560 documentation, and that the names of AT&T Bell Laboratories or
4561 Bellcore or any of their entities not be used in advertising or
4562 publicity pertaining to distribution of the software without
4563 specific, written prior permission.
4565 AT&T and Bellcore disclaim all warranties with regard to this
4566 software, including all implied warranties of merchantability
4567 and fitness. In no event shall AT&T or Bellcore be liable for
4568 any special, indirect or consequential damages or any damages
4569 whatsoever resulting from loss of use, data or profits, whether
4570 in an action of contract, negligence or other tortious action,
4571 arising out of or in connection with the use or performance of
4573 ****************************************************************/
4575 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4576 /* Divide the first argument by the product of the remaining
4577 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4579 #define FUNC_NAME s_divide
4581 scm_i_divide (SCM x
, SCM y
, int inexact
)
4588 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4589 else if (SCM_I_INUMP (x
))
4591 long xx
= SCM_I_INUM (x
);
4592 if (xx
== 1 || xx
== -1)
4594 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4596 scm_num_overflow (s_divide
);
4601 return scm_from_double (1.0 / (double) xx
);
4602 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4605 else if (SCM_BIGP (x
))
4608 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4609 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4611 else if (SCM_REALP (x
))
4613 double xx
= SCM_REAL_VALUE (x
);
4614 #ifndef ALLOW_DIVIDE_BY_ZERO
4616 scm_num_overflow (s_divide
);
4619 return scm_from_double (1.0 / xx
);
4621 else if (SCM_COMPLEXP (x
))
4623 double r
= SCM_COMPLEX_REAL (x
);
4624 double i
= SCM_COMPLEX_IMAG (x
);
4628 double d
= i
* (1.0 + t
* t
);
4629 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4634 double d
= r
* (1.0 + t
* t
);
4635 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4638 else if (SCM_FRACTIONP (x
))
4639 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4640 SCM_FRACTION_NUMERATOR (x
));
4642 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4645 if (SCM_I_INUMP (x
))
4647 long xx
= SCM_I_INUM (x
);
4648 if (SCM_I_INUMP (y
))
4650 long yy
= SCM_I_INUM (y
);
4653 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4654 scm_num_overflow (s_divide
);
4656 return scm_from_double ((double) xx
/ (double) yy
);
4659 else if (xx
% yy
!= 0)
4662 return scm_from_double ((double) xx
/ (double) yy
);
4663 else return scm_i_make_ratio (x
, y
);
4668 if (SCM_FIXABLE (z
))
4669 return SCM_I_MAKINUM (z
);
4671 return scm_i_long2big (z
);
4674 else if (SCM_BIGP (y
))
4677 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4678 else return scm_i_make_ratio (x
, y
);
4680 else if (SCM_REALP (y
))
4682 double yy
= SCM_REAL_VALUE (y
);
4683 #ifndef ALLOW_DIVIDE_BY_ZERO
4685 scm_num_overflow (s_divide
);
4688 return scm_from_double ((double) xx
/ yy
);
4690 else if (SCM_COMPLEXP (y
))
4693 complex_div
: /* y _must_ be a complex number */
4695 double r
= SCM_COMPLEX_REAL (y
);
4696 double i
= SCM_COMPLEX_IMAG (y
);
4700 double d
= i
* (1.0 + t
* t
);
4701 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4706 double d
= r
* (1.0 + t
* t
);
4707 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4711 else if (SCM_FRACTIONP (y
))
4712 /* a / b/c = ac / b */
4713 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4714 SCM_FRACTION_NUMERATOR (y
));
4716 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4718 else if (SCM_BIGP (x
))
4720 if (SCM_I_INUMP (y
))
4722 long int yy
= SCM_I_INUM (y
);
4725 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4726 scm_num_overflow (s_divide
);
4728 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4729 scm_remember_upto_here_1 (x
);
4730 return (sgn
== 0) ? scm_nan () : scm_inf ();
4737 /* FIXME: HMM, what are the relative performance issues here?
4738 We need to test. Is it faster on average to test
4739 divisible_p, then perform whichever operation, or is it
4740 faster to perform the integer div opportunistically and
4741 switch to real if there's a remainder? For now we take the
4742 middle ground: test, then if divisible, use the faster div
4745 long abs_yy
= yy
< 0 ? -yy
: yy
;
4746 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4750 SCM result
= scm_i_mkbig ();
4751 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4752 scm_remember_upto_here_1 (x
);
4754 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4755 return scm_i_normbig (result
);
4760 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4761 else return scm_i_make_ratio (x
, y
);
4765 else if (SCM_BIGP (y
))
4767 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4770 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4771 scm_num_overflow (s_divide
);
4773 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4774 scm_remember_upto_here_1 (x
);
4775 return (sgn
== 0) ? scm_nan () : scm_inf ();
4781 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4785 SCM result
= scm_i_mkbig ();
4786 mpz_divexact (SCM_I_BIG_MPZ (result
),
4789 scm_remember_upto_here_2 (x
, y
);
4790 return scm_i_normbig (result
);
4796 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4797 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4798 scm_remember_upto_here_2 (x
, y
);
4799 return scm_from_double (dbx
/ dby
);
4801 else return scm_i_make_ratio (x
, y
);
4805 else if (SCM_REALP (y
))
4807 double yy
= SCM_REAL_VALUE (y
);
4808 #ifndef ALLOW_DIVIDE_BY_ZERO
4810 scm_num_overflow (s_divide
);
4813 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4815 else if (SCM_COMPLEXP (y
))
4817 a
= scm_i_big2dbl (x
);
4820 else if (SCM_FRACTIONP (y
))
4821 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4822 SCM_FRACTION_NUMERATOR (y
));
4824 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4826 else if (SCM_REALP (x
))
4828 double rx
= SCM_REAL_VALUE (x
);
4829 if (SCM_I_INUMP (y
))
4831 long int yy
= SCM_I_INUM (y
);
4832 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4834 scm_num_overflow (s_divide
);
4837 return scm_from_double (rx
/ (double) yy
);
4839 else if (SCM_BIGP (y
))
4841 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4842 scm_remember_upto_here_1 (y
);
4843 return scm_from_double (rx
/ dby
);
4845 else if (SCM_REALP (y
))
4847 double yy
= SCM_REAL_VALUE (y
);
4848 #ifndef ALLOW_DIVIDE_BY_ZERO
4850 scm_num_overflow (s_divide
);
4853 return scm_from_double (rx
/ yy
);
4855 else if (SCM_COMPLEXP (y
))
4860 else if (SCM_FRACTIONP (y
))
4861 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4863 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4865 else if (SCM_COMPLEXP (x
))
4867 double rx
= SCM_COMPLEX_REAL (x
);
4868 double ix
= SCM_COMPLEX_IMAG (x
);
4869 if (SCM_I_INUMP (y
))
4871 long int yy
= SCM_I_INUM (y
);
4872 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4874 scm_num_overflow (s_divide
);
4879 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4882 else if (SCM_BIGP (y
))
4884 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4885 scm_remember_upto_here_1 (y
);
4886 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4888 else if (SCM_REALP (y
))
4890 double yy
= SCM_REAL_VALUE (y
);
4891 #ifndef ALLOW_DIVIDE_BY_ZERO
4893 scm_num_overflow (s_divide
);
4896 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4898 else if (SCM_COMPLEXP (y
))
4900 double ry
= SCM_COMPLEX_REAL (y
);
4901 double iy
= SCM_COMPLEX_IMAG (y
);
4905 double d
= iy
* (1.0 + t
* t
);
4906 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4911 double d
= ry
* (1.0 + t
* t
);
4912 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4915 else if (SCM_FRACTIONP (y
))
4917 double yy
= scm_i_fraction2double (y
);
4918 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4921 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4923 else if (SCM_FRACTIONP (x
))
4925 if (SCM_I_INUMP (y
))
4927 long int yy
= SCM_I_INUM (y
);
4928 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4930 scm_num_overflow (s_divide
);
4933 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4934 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4936 else if (SCM_BIGP (y
))
4938 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4939 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4941 else if (SCM_REALP (y
))
4943 double yy
= SCM_REAL_VALUE (y
);
4944 #ifndef ALLOW_DIVIDE_BY_ZERO
4946 scm_num_overflow (s_divide
);
4949 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4951 else if (SCM_COMPLEXP (y
))
4953 a
= scm_i_fraction2double (x
);
4956 else if (SCM_FRACTIONP (y
))
4957 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4958 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4960 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4963 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4967 scm_divide (SCM x
, SCM y
)
4969 return scm_i_divide (x
, y
, 0);
4972 static SCM
scm_divide2real (SCM x
, SCM y
)
4974 return scm_i_divide (x
, y
, 1);
4980 scm_asinh (double x
)
4985 #define asinh scm_asinh
4986 return log (x
+ sqrt (x
* x
+ 1));
4989 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4990 /* "Return the inverse hyperbolic sine of @var{x}."
4995 scm_acosh (double x
)
5000 #define acosh scm_acosh
5001 return log (x
+ sqrt (x
* x
- 1));
5004 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5005 /* "Return the inverse hyperbolic cosine of @var{x}."
5010 scm_atanh (double x
)
5015 #define atanh scm_atanh
5016 return 0.5 * log ((1 + x
) / (1 - x
));
5019 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5020 /* "Return the inverse hyperbolic tangent of @var{x}."
5025 scm_c_truncate (double x
)
5036 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5037 half-way case (ie. when x is an integer plus 0.5) going upwards.
5038 Then half-way cases are identified and adjusted down if the
5039 round-upwards didn't give the desired even integer.
5041 "plus_half == result" identifies a half-way case. If plus_half, which is
5042 x + 0.5, is an integer then x must be an integer plus 0.5.
5044 An odd "result" value is identified with result/2 != floor(result/2).
5045 This is done with plus_half, since that value is ready for use sooner in
5046 a pipelined cpu, and we're already requiring plus_half == result.
5048 Note however that we need to be careful when x is big and already an
5049 integer. In that case "x+0.5" may round to an adjacent integer, causing
5050 us to return such a value, incorrectly. For instance if the hardware is
5051 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5052 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5053 returned. Or if the hardware is in round-upwards mode, then other bigger
5054 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5055 representable value, 2^128+2^76 (or whatever), again incorrect.
5057 These bad roundings of x+0.5 are avoided by testing at the start whether
5058 x is already an integer. If it is then clearly that's the desired result
5059 already. And if it's not then the exponent must be small enough to allow
5060 an 0.5 to be represented, and hence added without a bad rounding. */
5063 scm_c_round (double x
)
5065 double plus_half
, result
;
5070 plus_half
= x
+ 0.5;
5071 result
= floor (plus_half
);
5072 /* Adjust so that the rounding is towards even. */
5073 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5078 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5080 "Round the number @var{x} towards zero.")
5081 #define FUNC_NAME s_scm_truncate_number
5083 if (scm_is_false (scm_negative_p (x
)))
5084 return scm_floor (x
);
5086 return scm_ceiling (x
);
5090 static SCM exactly_one_half
;
5092 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5094 "Round the number @var{x} towards the nearest integer. "
5095 "When it is exactly halfway between two integers, "
5096 "round towards the even one.")
5097 #define FUNC_NAME s_scm_round_number
5099 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5101 else if (SCM_REALP (x
))
5102 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5105 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5106 single quotient+remainder division then examining to see which way
5107 the rounding should go. */
5108 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5109 SCM result
= scm_floor (plus_half
);
5110 /* Adjust so that the rounding is towards even. */
5111 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5112 && scm_is_true (scm_odd_p (result
)))
5113 return scm_difference (result
, SCM_I_MAKINUM (1));
5120 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5122 "Round the number @var{x} towards minus infinity.")
5123 #define FUNC_NAME s_scm_floor
5125 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5127 else if (SCM_REALP (x
))
5128 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5129 else if (SCM_FRACTIONP (x
))
5131 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5132 SCM_FRACTION_DENOMINATOR (x
));
5133 if (scm_is_false (scm_negative_p (x
)))
5135 /* For positive x, rounding towards zero is correct. */
5140 /* For negative x, we need to return q-1 unless x is an
5141 integer. But fractions are never integer, per our
5143 return scm_difference (q
, SCM_I_MAKINUM (1));
5147 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5151 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5153 "Round the number @var{x} towards infinity.")
5154 #define FUNC_NAME s_scm_ceiling
5156 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5158 else if (SCM_REALP (x
))
5159 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5160 else if (SCM_FRACTIONP (x
))
5162 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5163 SCM_FRACTION_DENOMINATOR (x
));
5164 if (scm_is_false (scm_positive_p (x
)))
5166 /* For negative x, rounding towards zero is correct. */
5171 /* For positive x, we need to return q+1 unless x is an
5172 integer. But fractions are never integer, per our
5174 return scm_sum (q
, SCM_I_MAKINUM (1));
5178 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5182 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5183 /* "Return the square root of the real number @var{x}."
5185 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5186 /* "Return the absolute value of the real number @var{x}."
5188 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5189 /* "Return the @var{x}th power of e."
5191 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5192 /* "Return the natural logarithm of the real number @var{x}."
5194 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5195 /* "Return the sine of the real number @var{x}."
5197 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5198 /* "Return the cosine of the real number @var{x}."
5200 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5201 /* "Return the tangent of the real number @var{x}."
5203 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5204 /* "Return the arc sine of the real number @var{x}."
5206 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5207 /* "Return the arc cosine of the real number @var{x}."
5209 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5210 /* "Return the arc tangent of the real number @var{x}."
5212 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5213 /* "Return the hyperbolic sine of the real number @var{x}."
5215 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5216 /* "Return the hyperbolic cosine of the real number @var{x}."
5218 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5219 /* "Return the hyperbolic tangent of the real number @var{x}."
5227 static void scm_two_doubles (SCM x
,
5229 const char *sstring
,
5233 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5235 if (SCM_I_INUMP (x
))
5236 xy
->x
= SCM_I_INUM (x
);
5237 else if (SCM_BIGP (x
))
5238 xy
->x
= scm_i_big2dbl (x
);
5239 else if (SCM_REALP (x
))
5240 xy
->x
= SCM_REAL_VALUE (x
);
5241 else if (SCM_FRACTIONP (x
))
5242 xy
->x
= scm_i_fraction2double (x
);
5244 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5246 if (SCM_I_INUMP (y
))
5247 xy
->y
= SCM_I_INUM (y
);
5248 else if (SCM_BIGP (y
))
5249 xy
->y
= scm_i_big2dbl (y
);
5250 else if (SCM_REALP (y
))
5251 xy
->y
= SCM_REAL_VALUE (y
);
5252 else if (SCM_FRACTIONP (y
))
5253 xy
->y
= scm_i_fraction2double (y
);
5255 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5259 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5261 "Return @var{x} raised to the power of @var{y}. This\n"
5262 "procedure does not accept complex arguments.")
5263 #define FUNC_NAME s_scm_sys_expt
5266 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5267 return scm_from_double (pow (xy
.x
, xy
.y
));
5272 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5274 "Return the arc tangent of the two arguments @var{x} and\n"
5275 "@var{y}. This is similar to calculating the arc tangent of\n"
5276 "@var{x} / @var{y}, except that the signs of both arguments\n"
5277 "are used to determine the quadrant of the result. This\n"
5278 "procedure does not accept complex arguments.")
5279 #define FUNC_NAME s_scm_sys_atan2
5282 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5283 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5288 scm_c_make_rectangular (double re
, double im
)
5291 return scm_from_double (re
);
5295 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5297 SCM_COMPLEX_REAL (z
) = re
;
5298 SCM_COMPLEX_IMAG (z
) = im
;
5303 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5304 (SCM real
, SCM imaginary
),
5305 "Return a complex number constructed of the given @var{real} and\n"
5306 "@var{imaginary} parts.")
5307 #define FUNC_NAME s_scm_make_rectangular
5310 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5311 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5316 scm_c_make_polar (double mag
, double ang
)
5320 sincos (ang
, &s
, &c
);
5325 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5328 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5330 "Return the complex number @var{x} * e^(i * @var{y}).")
5331 #define FUNC_NAME s_scm_make_polar
5334 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5335 return scm_c_make_polar (xy
.x
, xy
.y
);
5340 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5341 /* "Return the real part of the number @var{z}."
5344 scm_real_part (SCM z
)
5346 if (SCM_I_INUMP (z
))
5348 else if (SCM_BIGP (z
))
5350 else if (SCM_REALP (z
))
5352 else if (SCM_COMPLEXP (z
))
5353 return scm_from_double (SCM_COMPLEX_REAL (z
));
5354 else if (SCM_FRACTIONP (z
))
5357 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5361 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5362 /* "Return the imaginary part of the number @var{z}."
5365 scm_imag_part (SCM z
)
5367 if (SCM_I_INUMP (z
))
5369 else if (SCM_BIGP (z
))
5371 else if (SCM_REALP (z
))
5373 else if (SCM_COMPLEXP (z
))
5374 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5375 else if (SCM_FRACTIONP (z
))
5378 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5381 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5382 /* "Return the numerator of the number @var{z}."
5385 scm_numerator (SCM z
)
5387 if (SCM_I_INUMP (z
))
5389 else if (SCM_BIGP (z
))
5391 else if (SCM_FRACTIONP (z
))
5393 scm_i_fraction_reduce (z
);
5394 return SCM_FRACTION_NUMERATOR (z
);
5396 else if (SCM_REALP (z
))
5397 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5399 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5403 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5404 /* "Return the denominator of the number @var{z}."
5407 scm_denominator (SCM z
)
5409 if (SCM_I_INUMP (z
))
5410 return SCM_I_MAKINUM (1);
5411 else if (SCM_BIGP (z
))
5412 return SCM_I_MAKINUM (1);
5413 else if (SCM_FRACTIONP (z
))
5415 scm_i_fraction_reduce (z
);
5416 return SCM_FRACTION_DENOMINATOR (z
);
5418 else if (SCM_REALP (z
))
5419 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5421 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5424 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5425 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5426 * "@code{abs} for real arguments, but also allows complex numbers."
5429 scm_magnitude (SCM z
)
5431 if (SCM_I_INUMP (z
))
5433 long int zz
= SCM_I_INUM (z
);
5436 else if (SCM_POSFIXABLE (-zz
))
5437 return SCM_I_MAKINUM (-zz
);
5439 return scm_i_long2big (-zz
);
5441 else if (SCM_BIGP (z
))
5443 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5444 scm_remember_upto_here_1 (z
);
5446 return scm_i_clonebig (z
, 0);
5450 else if (SCM_REALP (z
))
5451 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5452 else if (SCM_COMPLEXP (z
))
5453 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5454 else if (SCM_FRACTIONP (z
))
5456 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5458 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5459 SCM_FRACTION_DENOMINATOR (z
));
5462 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5466 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5467 /* "Return the angle of the complex number @var{z}."
5472 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5473 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5474 But if atan2 follows the floating point rounding mode, then the value
5475 is not a constant. Maybe it'd be close enough though. */
5476 if (SCM_I_INUMP (z
))
5478 if (SCM_I_INUM (z
) >= 0)
5481 return scm_from_double (atan2 (0.0, -1.0));
5483 else if (SCM_BIGP (z
))
5485 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5486 scm_remember_upto_here_1 (z
);
5488 return scm_from_double (atan2 (0.0, -1.0));
5492 else if (SCM_REALP (z
))
5494 if (SCM_REAL_VALUE (z
) >= 0)
5497 return scm_from_double (atan2 (0.0, -1.0));
5499 else if (SCM_COMPLEXP (z
))
5500 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5501 else if (SCM_FRACTIONP (z
))
5503 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5505 else return scm_from_double (atan2 (0.0, -1.0));
5508 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5512 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5513 /* Convert the number @var{x} to its inexact representation.\n"
5516 scm_exact_to_inexact (SCM z
)
5518 if (SCM_I_INUMP (z
))
5519 return scm_from_double ((double) SCM_I_INUM (z
));
5520 else if (SCM_BIGP (z
))
5521 return scm_from_double (scm_i_big2dbl (z
));
5522 else if (SCM_FRACTIONP (z
))
5523 return scm_from_double (scm_i_fraction2double (z
));
5524 else if (SCM_INEXACTP (z
))
5527 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5531 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5533 "Return an exact number that is numerically closest to @var{z}.")
5534 #define FUNC_NAME s_scm_inexact_to_exact
5536 if (SCM_I_INUMP (z
))
5538 else if (SCM_BIGP (z
))
5540 else if (SCM_REALP (z
))
5542 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5543 SCM_OUT_OF_RANGE (1, z
);
5550 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5551 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5552 scm_i_mpz2num (mpq_denref (frac
)));
5554 /* When scm_i_make_ratio throws, we leak the memory allocated
5561 else if (SCM_FRACTIONP (z
))
5564 SCM_WRONG_TYPE_ARG (1, z
);
5568 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5570 "Return an exact number that is within @var{err} of @var{x}.")
5571 #define FUNC_NAME s_scm_rationalize
5573 if (SCM_I_INUMP (x
))
5575 else if (SCM_BIGP (x
))
5577 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5579 /* Use continued fractions to find closest ratio. All
5580 arithmetic is done with exact numbers.
5583 SCM ex
= scm_inexact_to_exact (x
);
5584 SCM int_part
= scm_floor (ex
);
5585 SCM tt
= SCM_I_MAKINUM (1);
5586 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5587 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5591 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5594 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5595 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5597 /* We stop after a million iterations just to be absolutely sure
5598 that we don't go into an infinite loop. The process normally
5599 converges after less than a dozen iterations.
5602 err
= scm_abs (err
);
5603 while (++i
< 1000000)
5605 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5606 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5607 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5609 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5610 err
))) /* abs(x-a/b) <= err */
5612 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5613 if (scm_is_false (scm_exact_p (x
))
5614 || scm_is_false (scm_exact_p (err
)))
5615 return scm_exact_to_inexact (res
);
5619 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5621 tt
= scm_floor (rx
); /* tt = floor (rx) */
5627 scm_num_overflow (s_scm_rationalize
);
5630 SCM_WRONG_TYPE_ARG (1, x
);
5634 /* conversion functions */
5637 scm_is_integer (SCM val
)
5639 return scm_is_true (scm_integer_p (val
));
5643 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5645 if (SCM_I_INUMP (val
))
5647 scm_t_signed_bits n
= SCM_I_INUM (val
);
5648 return n
>= min
&& n
<= max
;
5650 else if (SCM_BIGP (val
))
5652 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5654 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5656 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5658 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5659 return n
>= min
&& n
<= max
;
5669 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5670 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5673 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5674 SCM_I_BIG_MPZ (val
));
5676 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5688 return n
>= min
&& n
<= max
;
5696 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5698 if (SCM_I_INUMP (val
))
5700 scm_t_signed_bits n
= SCM_I_INUM (val
);
5701 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5703 else if (SCM_BIGP (val
))
5705 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5707 else if (max
<= ULONG_MAX
)
5709 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5711 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5712 return n
>= min
&& n
<= max
;
5722 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5725 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5726 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5729 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5730 SCM_I_BIG_MPZ (val
));
5732 return n
>= min
&& n
<= max
;
5740 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5742 scm_error (scm_out_of_range_key
,
5744 "Value out of range ~S to ~S: ~S",
5745 scm_list_3 (min
, max
, bad_val
),
5746 scm_list_1 (bad_val
));
5749 #define TYPE scm_t_intmax
5750 #define TYPE_MIN min
5751 #define TYPE_MAX max
5752 #define SIZEOF_TYPE 0
5753 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5754 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5755 #include "libguile/conv-integer.i.c"
5757 #define TYPE scm_t_uintmax
5758 #define TYPE_MIN min
5759 #define TYPE_MAX max
5760 #define SIZEOF_TYPE 0
5761 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5762 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5763 #include "libguile/conv-uinteger.i.c"
5765 #define TYPE scm_t_int8
5766 #define TYPE_MIN SCM_T_INT8_MIN
5767 #define TYPE_MAX SCM_T_INT8_MAX
5768 #define SIZEOF_TYPE 1
5769 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5770 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5771 #include "libguile/conv-integer.i.c"
5773 #define TYPE scm_t_uint8
5775 #define TYPE_MAX SCM_T_UINT8_MAX
5776 #define SIZEOF_TYPE 1
5777 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5778 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5779 #include "libguile/conv-uinteger.i.c"
5781 #define TYPE scm_t_int16
5782 #define TYPE_MIN SCM_T_INT16_MIN
5783 #define TYPE_MAX SCM_T_INT16_MAX
5784 #define SIZEOF_TYPE 2
5785 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5786 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5787 #include "libguile/conv-integer.i.c"
5789 #define TYPE scm_t_uint16
5791 #define TYPE_MAX SCM_T_UINT16_MAX
5792 #define SIZEOF_TYPE 2
5793 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5794 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5795 #include "libguile/conv-uinteger.i.c"
5797 #define TYPE scm_t_int32
5798 #define TYPE_MIN SCM_T_INT32_MIN
5799 #define TYPE_MAX SCM_T_INT32_MAX
5800 #define SIZEOF_TYPE 4
5801 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5802 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5803 #include "libguile/conv-integer.i.c"
5805 #define TYPE scm_t_uint32
5807 #define TYPE_MAX SCM_T_UINT32_MAX
5808 #define SIZEOF_TYPE 4
5809 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5810 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5811 #include "libguile/conv-uinteger.i.c"
5813 #if SCM_HAVE_T_INT64
5815 #define TYPE scm_t_int64
5816 #define TYPE_MIN SCM_T_INT64_MIN
5817 #define TYPE_MAX SCM_T_INT64_MAX
5818 #define SIZEOF_TYPE 8
5819 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5820 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5821 #include "libguile/conv-integer.i.c"
5823 #define TYPE scm_t_uint64
5825 #define TYPE_MAX SCM_T_UINT64_MAX
5826 #define SIZEOF_TYPE 8
5827 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5828 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5829 #include "libguile/conv-uinteger.i.c"
5834 scm_to_mpz (SCM val
, mpz_t rop
)
5836 if (SCM_I_INUMP (val
))
5837 mpz_set_si (rop
, SCM_I_INUM (val
));
5838 else if (SCM_BIGP (val
))
5839 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5841 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5845 scm_from_mpz (mpz_t val
)
5847 return scm_i_mpz2num (val
);
5851 scm_is_real (SCM val
)
5853 return scm_is_true (scm_real_p (val
));
5857 scm_is_rational (SCM val
)
5859 return scm_is_true (scm_rational_p (val
));
5863 scm_to_double (SCM val
)
5865 if (SCM_I_INUMP (val
))
5866 return SCM_I_INUM (val
);
5867 else if (SCM_BIGP (val
))
5868 return scm_i_big2dbl (val
);
5869 else if (SCM_FRACTIONP (val
))
5870 return scm_i_fraction2double (val
);
5871 else if (SCM_REALP (val
))
5872 return SCM_REAL_VALUE (val
);
5874 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5878 scm_from_double (double val
)
5880 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5881 SCM_REAL_VALUE (z
) = val
;
5885 #if SCM_ENABLE_DISCOURAGED == 1
5888 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5892 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5896 scm_out_of_range (NULL
, num
);
5899 return scm_to_double (num
);
5903 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5907 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5911 scm_out_of_range (NULL
, num
);
5914 return scm_to_double (num
);
5920 scm_is_complex (SCM val
)
5922 return scm_is_true (scm_complex_p (val
));
5926 scm_c_real_part (SCM z
)
5928 if (SCM_COMPLEXP (z
))
5929 return SCM_COMPLEX_REAL (z
);
5932 /* Use the scm_real_part to get proper error checking and
5935 return scm_to_double (scm_real_part (z
));
5940 scm_c_imag_part (SCM z
)
5942 if (SCM_COMPLEXP (z
))
5943 return SCM_COMPLEX_IMAG (z
);
5946 /* Use the scm_imag_part to get proper error checking and
5947 dispatching. The result will almost always be 0.0, but not
5950 return scm_to_double (scm_imag_part (z
));
5955 scm_c_magnitude (SCM z
)
5957 return scm_to_double (scm_magnitude (z
));
5963 return scm_to_double (scm_angle (z
));
5967 scm_is_number (SCM z
)
5969 return scm_is_true (scm_number_p (z
));
5977 mpz_init_set_si (z_negative_one
, -1);
5979 /* It may be possible to tune the performance of some algorithms by using
5980 * the following constants to avoid the creation of bignums. Please, before
5981 * using these values, remember the two rules of program optimization:
5982 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5983 scm_c_define ("most-positive-fixnum",
5984 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5985 scm_c_define ("most-negative-fixnum",
5986 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5988 scm_add_feature ("complex");
5989 scm_add_feature ("inexact");
5990 scm_flo0
= scm_from_double (0.0);
5992 /* determine floating point precision */
5993 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5995 init_dblprec(&scm_dblprec
[i
-2],i
);
5996 init_fx_radix(fx_per_radix
[i
-2],i
);
5999 /* hard code precision for base 10 if the preprocessor tells us to... */
6000 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6003 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6004 SCM_I_MAKINUM (2)));
6005 #include "libguile/numbers.x"