20fda02daae4041683d8a4e80e30e946fdcca24b
1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public License
9 * as published by the Free Software Foundation; either version 3 of
10 * the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
24 /* General assumptions:
25 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
26 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
27 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
28 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
29 * All objects satisfying SCM_FRACTIONP are never an integer.
34 - see if special casing bignums and reals in integer-exponent when
35 possible (to use mpz_pow and mpf_pow_ui) is faster.
37 - look in to better short-circuiting of common cases in
38 integer-expt and elsewhere.
40 - see if direct mpz operations can help in ash and elsewhere.
57 #include "libguile/_scm.h"
58 #include "libguile/feature.h"
59 #include "libguile/ports.h"
60 #include "libguile/root.h"
61 #include "libguile/smob.h"
62 #include "libguile/strings.h"
64 #include "libguile/validate.h"
65 #include "libguile/numbers.h"
66 #include "libguile/deprecation.h"
68 #include "libguile/eq.h"
70 #include "libguile/discouraged.h"
72 /* values per glibc, if not already defined */
74 #define M_LOG10E 0.43429448190325182765
77 #define M_PI 3.14159265358979323846
83 Wonder if this might be faster for some of our code? A switch on
84 the numtag would jump directly to the right case, and the
85 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
87 #define SCM_I_NUMTAG_NOTNUM 0
88 #define SCM_I_NUMTAG_INUM 1
89 #define SCM_I_NUMTAG_BIG scm_tc16_big
90 #define SCM_I_NUMTAG_REAL scm_tc16_real
91 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
92 #define SCM_I_NUMTAG(x) \
93 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
94 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
95 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
96 : SCM_I_NUMTAG_NOTNUM)))
98 /* the macro above will not work as is with fractions */
101 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
103 /* FLOBUFLEN is the maximum number of characters neccessary for the
104 * printed or scm_string representation of an inexact number.
106 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
109 #if ! defined (HAVE_ISNAN)
114 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
117 #if ! defined (HAVE_ISINF)
122 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
129 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
130 an explicit check. In some future gmp (don't know what version number),
131 mpz_cmp_d is supposed to do this itself. */
133 #define xmpz_cmp_d(z, d) \
134 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
136 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
139 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
140 isinf. It does have finite and isnan though, hence the use of those.
141 fpclass would be a possibility on that system too. */
145 #if defined (HAVE_ISINF)
147 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
148 return (! (finite (x
) || isnan (x
)));
157 #if defined (HAVE_ISNAN)
164 #if defined (GUILE_I)
165 #if HAVE_COMPLEX_DOUBLE
167 /* For an SCM object Z which is a complex number (ie. satisfies
168 SCM_COMPLEXP), return its value as a C level "complex double". */
169 #define SCM_COMPLEX_VALUE(z) \
170 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
172 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
174 /* Convert a C "complex double" to an SCM value. */
176 scm_from_complex_double (complex double z
)
178 return scm_c_make_rectangular (creal (z
), cimag (z
));
181 #endif /* HAVE_COMPLEX_DOUBLE */
186 static mpz_t z_negative_one
;
193 /* Return a newly created bignum. */
194 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
195 mpz_init (SCM_I_BIG_MPZ (z
));
200 scm_i_long2big (long x
)
202 /* Return a newly created bignum initialized to X. */
203 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
204 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
209 scm_i_ulong2big (unsigned long x
)
211 /* Return a newly created bignum initialized to X. */
212 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
213 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
218 scm_i_clonebig (SCM src_big
, int same_sign_p
)
220 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
221 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
222 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
224 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
229 scm_i_bigcmp (SCM x
, SCM y
)
231 /* Return neg if x < y, pos if x > y, and 0 if x == y */
232 /* presume we already know x and y are bignums */
233 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
234 scm_remember_upto_here_2 (x
, y
);
239 scm_i_dbl2big (double d
)
241 /* results are only defined if d is an integer */
242 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
243 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
247 /* Convert a integer in double representation to a SCM number. */
250 scm_i_dbl2num (double u
)
252 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
253 powers of 2, so there's no rounding when making "double" values
254 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
255 get rounded on a 64-bit machine, hence the "+1".
257 The use of floor() to force to an integer value ensures we get a
258 "numerically closest" value without depending on how a
259 double->long cast or how mpz_set_d will round. For reference,
260 double->long probably follows the hardware rounding mode,
261 mpz_set_d truncates towards zero. */
263 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
264 representable as a double? */
266 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
267 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
268 return SCM_I_MAKINUM ((long) u
);
270 return scm_i_dbl2big (u
);
273 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
274 with R5RS exact->inexact.
276 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
277 (ie. truncate towards zero), then adjust to get the closest double by
278 examining the next lower bit and adding 1 (to the absolute value) if
281 Bignums exactly half way between representable doubles are rounded to the
282 next higher absolute value (ie. away from zero). This seems like an
283 adequate interpretation of R5RS "numerically closest", and it's easier
284 and faster than a full "nearest-even" style.
286 The bit test must be done on the absolute value of the mpz_t, which means
287 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
288 negatives as twos complement.
290 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
291 following the hardware rounding mode, but applied to the absolute value
292 of the mpz_t operand. This is not what we want so we put the high
293 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
294 mpz_get_d is supposed to always truncate towards zero.
296 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
297 is a slowdown. It'd be faster to pick out the relevant high bits with
298 mpz_getlimbn if we could be bothered coding that, and if the new
299 truncating gmp doesn't come out. */
302 scm_i_big2dbl (SCM b
)
307 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
311 /* Current GMP, eg. 4.1.3, force truncation towards zero */
313 if (bits
> DBL_MANT_DIG
)
315 size_t shift
= bits
- DBL_MANT_DIG
;
316 mpz_init2 (tmp
, DBL_MANT_DIG
);
317 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
318 result
= ldexp (mpz_get_d (tmp
), shift
);
323 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
328 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
331 if (bits
> DBL_MANT_DIG
)
333 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
334 /* test bit number "pos" in absolute value */
335 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
336 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
338 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
342 scm_remember_upto_here_1 (b
);
347 scm_i_normbig (SCM b
)
349 /* convert a big back to a fixnum if it'll fit */
350 /* presume b is a bignum */
351 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
353 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
354 if (SCM_FIXABLE (val
))
355 b
= SCM_I_MAKINUM (val
);
360 static SCM_C_INLINE_KEYWORD SCM
361 scm_i_mpz2num (mpz_t b
)
363 /* convert a mpz number to a SCM number. */
364 if (mpz_fits_slong_p (b
))
366 long val
= mpz_get_si (b
);
367 if (SCM_FIXABLE (val
))
368 return SCM_I_MAKINUM (val
);
372 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
373 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
378 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
379 static SCM
scm_divide2real (SCM x
, SCM y
);
382 scm_i_make_ratio (SCM numerator
, SCM denominator
)
383 #define FUNC_NAME "make-ratio"
385 /* First make sure the arguments are proper.
387 if (SCM_I_INUMP (denominator
))
389 if (scm_is_eq (denominator
, SCM_INUM0
))
390 scm_num_overflow ("make-ratio");
391 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
396 if (!(SCM_BIGP(denominator
)))
397 SCM_WRONG_TYPE_ARG (2, denominator
);
399 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
400 SCM_WRONG_TYPE_ARG (1, numerator
);
402 /* Then flip signs so that the denominator is positive.
404 if (scm_is_true (scm_negative_p (denominator
)))
406 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
407 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
410 /* Now consider for each of the four fixnum/bignum combinations
411 whether the rational number is really an integer.
413 if (SCM_I_INUMP (numerator
))
415 long x
= SCM_I_INUM (numerator
);
416 if (scm_is_eq (numerator
, SCM_INUM0
))
418 if (SCM_I_INUMP (denominator
))
421 y
= SCM_I_INUM (denominator
);
423 return SCM_I_MAKINUM(1);
425 return SCM_I_MAKINUM (x
/ y
);
429 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
430 of that value for the denominator, as a bignum. Apart from
431 that case, abs(bignum) > abs(inum) so inum/bignum is not an
433 if (x
== SCM_MOST_NEGATIVE_FIXNUM
434 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
435 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
436 return SCM_I_MAKINUM(-1);
439 else if (SCM_BIGP (numerator
))
441 if (SCM_I_INUMP (denominator
))
443 long yy
= SCM_I_INUM (denominator
);
444 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
445 return scm_divide (numerator
, denominator
);
449 if (scm_is_eq (numerator
, denominator
))
450 return SCM_I_MAKINUM(1);
451 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
452 SCM_I_BIG_MPZ (denominator
)))
453 return scm_divide(numerator
, denominator
);
457 /* No, it's a proper fraction.
460 SCM divisor
= scm_gcd (numerator
, denominator
);
461 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
463 numerator
= scm_divide (numerator
, divisor
);
464 denominator
= scm_divide (denominator
, divisor
);
467 return scm_double_cell (scm_tc16_fraction
,
468 SCM_UNPACK (numerator
),
469 SCM_UNPACK (denominator
), 0);
475 scm_i_fraction2double (SCM z
)
477 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
478 SCM_FRACTION_DENOMINATOR (z
)));
481 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
483 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
485 #define FUNC_NAME s_scm_exact_p
491 if (SCM_FRACTIONP (x
))
495 SCM_WRONG_TYPE_ARG (1, x
);
500 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
502 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
504 #define FUNC_NAME s_scm_odd_p
508 long val
= SCM_I_INUM (n
);
509 return scm_from_bool ((val
& 1L) != 0);
511 else if (SCM_BIGP (n
))
513 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
514 scm_remember_upto_here_1 (n
);
515 return scm_from_bool (odd_p
);
517 else if (scm_is_true (scm_inf_p (n
)))
519 else if (SCM_REALP (n
))
521 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
527 SCM_WRONG_TYPE_ARG (1, n
);
530 SCM_WRONG_TYPE_ARG (1, n
);
535 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
537 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
539 #define FUNC_NAME s_scm_even_p
543 long val
= SCM_I_INUM (n
);
544 return scm_from_bool ((val
& 1L) == 0);
546 else if (SCM_BIGP (n
))
548 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
549 scm_remember_upto_here_1 (n
);
550 return scm_from_bool (even_p
);
552 else if (scm_is_true (scm_inf_p (n
)))
554 else if (SCM_REALP (n
))
556 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
562 SCM_WRONG_TYPE_ARG (1, n
);
565 SCM_WRONG_TYPE_ARG (1, n
);
569 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
571 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
572 "or @samp{-inf.0}, @code{#f} otherwise.")
573 #define FUNC_NAME s_scm_inf_p
576 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
577 else if (SCM_COMPLEXP (x
))
578 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
579 || xisinf (SCM_COMPLEX_IMAG (x
)));
585 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
587 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
589 #define FUNC_NAME s_scm_nan_p
592 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
593 else if (SCM_COMPLEXP (n
))
594 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
595 || xisnan (SCM_COMPLEX_IMAG (n
)));
601 /* Guile's idea of infinity. */
602 static double guile_Inf
;
604 /* Guile's idea of not a number. */
605 static double guile_NaN
;
608 guile_ieee_init (void)
610 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
612 /* Some version of gcc on some old version of Linux used to crash when
613 trying to make Inf and NaN. */
616 /* C99 INFINITY, when available.
617 FIXME: The standard allows for INFINITY to be something that overflows
618 at compile time. We ought to have a configure test to check for that
619 before trying to use it. (But in practice we believe this is not a
620 problem on any system guile is likely to target.) */
621 guile_Inf
= INFINITY
;
624 extern unsigned int DINFINITY
[2];
625 guile_Inf
= (*((double *) (DINFINITY
)));
632 if (guile_Inf
== tmp
)
640 #if defined (HAVE_ISNAN)
643 /* C99 NAN, when available */
648 extern unsigned int DQNAN
[2];
649 guile_NaN
= (*((double *)(DQNAN
)));
652 guile_NaN
= guile_Inf
/ guile_Inf
;
658 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
661 #define FUNC_NAME s_scm_inf
663 static int initialized
= 0;
669 return scm_from_double (guile_Inf
);
673 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
676 #define FUNC_NAME s_scm_nan
678 static int initialized
= 0;
684 return scm_from_double (guile_NaN
);
689 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
691 "Return the absolute value of @var{x}.")
696 long int xx
= SCM_I_INUM (x
);
699 else if (SCM_POSFIXABLE (-xx
))
700 return SCM_I_MAKINUM (-xx
);
702 return scm_i_long2big (-xx
);
704 else if (SCM_BIGP (x
))
706 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
708 return scm_i_clonebig (x
, 0);
712 else if (SCM_REALP (x
))
714 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
715 double xx
= SCM_REAL_VALUE (x
);
717 return scm_from_double (-xx
);
721 else if (SCM_FRACTIONP (x
))
723 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
725 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
726 SCM_FRACTION_DENOMINATOR (x
));
729 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
734 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
735 /* "Return the quotient of the numbers @var{x} and @var{y}."
738 scm_quotient (SCM x
, SCM y
)
742 long xx
= SCM_I_INUM (x
);
745 long yy
= SCM_I_INUM (y
);
747 scm_num_overflow (s_quotient
);
752 return SCM_I_MAKINUM (z
);
754 return scm_i_long2big (z
);
757 else if (SCM_BIGP (y
))
759 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
760 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
761 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
763 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
764 scm_remember_upto_here_1 (y
);
765 return SCM_I_MAKINUM (-1);
768 return SCM_I_MAKINUM (0);
771 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
773 else if (SCM_BIGP (x
))
777 long yy
= SCM_I_INUM (y
);
779 scm_num_overflow (s_quotient
);
784 SCM result
= scm_i_mkbig ();
787 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
790 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
793 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
794 scm_remember_upto_here_1 (x
);
795 return scm_i_normbig (result
);
798 else if (SCM_BIGP (y
))
800 SCM result
= scm_i_mkbig ();
801 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
804 scm_remember_upto_here_2 (x
, y
);
805 return scm_i_normbig (result
);
808 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
811 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
814 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
815 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
817 * "(remainder 13 4) @result{} 1\n"
818 * "(remainder -13 4) @result{} -1\n"
822 scm_remainder (SCM x
, SCM y
)
828 long yy
= SCM_I_INUM (y
);
830 scm_num_overflow (s_remainder
);
833 long z
= SCM_I_INUM (x
) % yy
;
834 return SCM_I_MAKINUM (z
);
837 else if (SCM_BIGP (y
))
839 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
840 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
841 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
843 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
844 scm_remember_upto_here_1 (y
);
845 return SCM_I_MAKINUM (0);
851 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
853 else if (SCM_BIGP (x
))
857 long yy
= SCM_I_INUM (y
);
859 scm_num_overflow (s_remainder
);
862 SCM result
= scm_i_mkbig ();
865 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
866 scm_remember_upto_here_1 (x
);
867 return scm_i_normbig (result
);
870 else if (SCM_BIGP (y
))
872 SCM result
= scm_i_mkbig ();
873 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
876 scm_remember_upto_here_2 (x
, y
);
877 return scm_i_normbig (result
);
880 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
883 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
887 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
888 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
890 * "(modulo 13 4) @result{} 1\n"
891 * "(modulo -13 4) @result{} 3\n"
895 scm_modulo (SCM x
, SCM y
)
899 long xx
= SCM_I_INUM (x
);
902 long yy
= SCM_I_INUM (y
);
904 scm_num_overflow (s_modulo
);
907 /* C99 specifies that "%" is the remainder corresponding to a
908 quotient rounded towards zero, and that's also traditional
909 for machine division, so z here should be well defined. */
927 return SCM_I_MAKINUM (result
);
930 else if (SCM_BIGP (y
))
932 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
939 SCM pos_y
= scm_i_clonebig (y
, 0);
940 /* do this after the last scm_op */
941 mpz_init_set_si (z_x
, xx
);
942 result
= pos_y
; /* re-use this bignum */
943 mpz_mod (SCM_I_BIG_MPZ (result
),
945 SCM_I_BIG_MPZ (pos_y
));
946 scm_remember_upto_here_1 (pos_y
);
950 result
= scm_i_mkbig ();
951 /* do this after the last scm_op */
952 mpz_init_set_si (z_x
, xx
);
953 mpz_mod (SCM_I_BIG_MPZ (result
),
956 scm_remember_upto_here_1 (y
);
959 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
960 mpz_add (SCM_I_BIG_MPZ (result
),
962 SCM_I_BIG_MPZ (result
));
963 scm_remember_upto_here_1 (y
);
964 /* and do this before the next one */
966 return scm_i_normbig (result
);
970 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
972 else if (SCM_BIGP (x
))
976 long yy
= SCM_I_INUM (y
);
978 scm_num_overflow (s_modulo
);
981 SCM result
= scm_i_mkbig ();
982 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
984 (yy
< 0) ? - yy
: yy
);
985 scm_remember_upto_here_1 (x
);
986 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
987 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
988 SCM_I_BIG_MPZ (result
),
990 return scm_i_normbig (result
);
993 else if (SCM_BIGP (y
))
996 SCM result
= scm_i_mkbig ();
997 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
998 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
999 mpz_mod (SCM_I_BIG_MPZ (result
),
1001 SCM_I_BIG_MPZ (pos_y
));
1003 scm_remember_upto_here_1 (x
);
1004 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1005 mpz_add (SCM_I_BIG_MPZ (result
),
1007 SCM_I_BIG_MPZ (result
));
1008 scm_remember_upto_here_2 (y
, pos_y
);
1009 return scm_i_normbig (result
);
1013 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1016 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1019 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
1020 /* "Return the greatest common divisor of all arguments.\n"
1021 * "If called without arguments, 0 is returned."
1024 scm_gcd (SCM x
, SCM y
)
1027 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1029 if (SCM_I_INUMP (x
))
1031 if (SCM_I_INUMP (y
))
1033 long xx
= SCM_I_INUM (x
);
1034 long yy
= SCM_I_INUM (y
);
1035 long u
= xx
< 0 ? -xx
: xx
;
1036 long v
= yy
< 0 ? -yy
: yy
;
1046 /* Determine a common factor 2^k */
1047 while (!(1 & (u
| v
)))
1053 /* Now, any factor 2^n can be eliminated */
1073 return (SCM_POSFIXABLE (result
)
1074 ? SCM_I_MAKINUM (result
)
1075 : scm_i_long2big (result
));
1077 else if (SCM_BIGP (y
))
1083 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1085 else if (SCM_BIGP (x
))
1087 if (SCM_I_INUMP (y
))
1089 unsigned long result
;
1092 yy
= SCM_I_INUM (y
);
1097 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1098 scm_remember_upto_here_1 (x
);
1099 return (SCM_POSFIXABLE (result
)
1100 ? SCM_I_MAKINUM (result
)
1101 : scm_from_ulong (result
));
1103 else if (SCM_BIGP (y
))
1105 SCM result
= scm_i_mkbig ();
1106 mpz_gcd (SCM_I_BIG_MPZ (result
),
1109 scm_remember_upto_here_2 (x
, y
);
1110 return scm_i_normbig (result
);
1113 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1116 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1119 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1120 /* "Return the least common multiple of the arguments.\n"
1121 * "If called without arguments, 1 is returned."
1124 scm_lcm (SCM n1
, SCM n2
)
1126 if (SCM_UNBNDP (n2
))
1128 if (SCM_UNBNDP (n1
))
1129 return SCM_I_MAKINUM (1L);
1130 n2
= SCM_I_MAKINUM (1L);
1133 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1134 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1135 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1136 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1138 if (SCM_I_INUMP (n1
))
1140 if (SCM_I_INUMP (n2
))
1142 SCM d
= scm_gcd (n1
, n2
);
1143 if (scm_is_eq (d
, SCM_INUM0
))
1146 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1150 /* inum n1, big n2 */
1153 SCM result
= scm_i_mkbig ();
1154 long nn1
= SCM_I_INUM (n1
);
1155 if (nn1
== 0) return SCM_INUM0
;
1156 if (nn1
< 0) nn1
= - nn1
;
1157 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1158 scm_remember_upto_here_1 (n2
);
1166 if (SCM_I_INUMP (n2
))
1173 SCM result
= scm_i_mkbig ();
1174 mpz_lcm(SCM_I_BIG_MPZ (result
),
1176 SCM_I_BIG_MPZ (n2
));
1177 scm_remember_upto_here_2(n1
, n2
);
1178 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1184 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1189 + + + x (map digit:logand X Y)
1190 + - + x (map digit:logand X (lognot (+ -1 Y)))
1191 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1192 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1197 + + + (map digit:logior X Y)
1198 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1199 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1200 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1205 + + + (map digit:logxor X Y)
1206 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1207 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1208 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1213 + + (any digit:logand X Y)
1214 + - (any digit:logand X (lognot (+ -1 Y)))
1215 - + (any digit:logand (lognot (+ -1 X)) Y)
1220 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1222 "Return the bitwise AND of the integer arguments.\n\n"
1224 "(logand) @result{} -1\n"
1225 "(logand 7) @result{} 7\n"
1226 "(logand #b111 #b011 #b001) @result{} 1\n"
1228 #define FUNC_NAME s_scm_logand
1232 if (SCM_UNBNDP (n2
))
1234 if (SCM_UNBNDP (n1
))
1235 return SCM_I_MAKINUM (-1);
1236 else if (!SCM_NUMBERP (n1
))
1237 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1238 else if (SCM_NUMBERP (n1
))
1241 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1244 if (SCM_I_INUMP (n1
))
1246 nn1
= SCM_I_INUM (n1
);
1247 if (SCM_I_INUMP (n2
))
1249 long nn2
= SCM_I_INUM (n2
);
1250 return SCM_I_MAKINUM (nn1
& nn2
);
1252 else if SCM_BIGP (n2
)
1258 SCM result_z
= scm_i_mkbig ();
1260 mpz_init_set_si (nn1_z
, nn1
);
1261 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1262 scm_remember_upto_here_1 (n2
);
1264 return scm_i_normbig (result_z
);
1268 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1270 else if (SCM_BIGP (n1
))
1272 if (SCM_I_INUMP (n2
))
1275 nn1
= SCM_I_INUM (n1
);
1278 else if (SCM_BIGP (n2
))
1280 SCM result_z
= scm_i_mkbig ();
1281 mpz_and (SCM_I_BIG_MPZ (result_z
),
1283 SCM_I_BIG_MPZ (n2
));
1284 scm_remember_upto_here_2 (n1
, n2
);
1285 return scm_i_normbig (result_z
);
1288 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1291 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1296 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1298 "Return the bitwise OR of the integer arguments.\n\n"
1300 "(logior) @result{} 0\n"
1301 "(logior 7) @result{} 7\n"
1302 "(logior #b000 #b001 #b011) @result{} 3\n"
1304 #define FUNC_NAME s_scm_logior
1308 if (SCM_UNBNDP (n2
))
1310 if (SCM_UNBNDP (n1
))
1312 else if (SCM_NUMBERP (n1
))
1315 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1318 if (SCM_I_INUMP (n1
))
1320 nn1
= SCM_I_INUM (n1
);
1321 if (SCM_I_INUMP (n2
))
1323 long nn2
= SCM_I_INUM (n2
);
1324 return SCM_I_MAKINUM (nn1
| nn2
);
1326 else if (SCM_BIGP (n2
))
1332 SCM result_z
= scm_i_mkbig ();
1334 mpz_init_set_si (nn1_z
, nn1
);
1335 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1336 scm_remember_upto_here_1 (n2
);
1338 return scm_i_normbig (result_z
);
1342 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1344 else if (SCM_BIGP (n1
))
1346 if (SCM_I_INUMP (n2
))
1349 nn1
= SCM_I_INUM (n1
);
1352 else if (SCM_BIGP (n2
))
1354 SCM result_z
= scm_i_mkbig ();
1355 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1357 SCM_I_BIG_MPZ (n2
));
1358 scm_remember_upto_here_2 (n1
, n2
);
1359 return scm_i_normbig (result_z
);
1362 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1365 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1370 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1372 "Return the bitwise XOR of the integer arguments. A bit is\n"
1373 "set in the result if it is set in an odd number of arguments.\n"
1375 "(logxor) @result{} 0\n"
1376 "(logxor 7) @result{} 7\n"
1377 "(logxor #b000 #b001 #b011) @result{} 2\n"
1378 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1380 #define FUNC_NAME s_scm_logxor
1384 if (SCM_UNBNDP (n2
))
1386 if (SCM_UNBNDP (n1
))
1388 else if (SCM_NUMBERP (n1
))
1391 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1394 if (SCM_I_INUMP (n1
))
1396 nn1
= SCM_I_INUM (n1
);
1397 if (SCM_I_INUMP (n2
))
1399 long nn2
= SCM_I_INUM (n2
);
1400 return SCM_I_MAKINUM (nn1
^ nn2
);
1402 else if (SCM_BIGP (n2
))
1406 SCM result_z
= scm_i_mkbig ();
1408 mpz_init_set_si (nn1_z
, nn1
);
1409 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1410 scm_remember_upto_here_1 (n2
);
1412 return scm_i_normbig (result_z
);
1416 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1418 else if (SCM_BIGP (n1
))
1420 if (SCM_I_INUMP (n2
))
1423 nn1
= SCM_I_INUM (n1
);
1426 else if (SCM_BIGP (n2
))
1428 SCM result_z
= scm_i_mkbig ();
1429 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1431 SCM_I_BIG_MPZ (n2
));
1432 scm_remember_upto_here_2 (n1
, n2
);
1433 return scm_i_normbig (result_z
);
1436 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1439 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1444 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1446 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1447 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1448 "without actually calculating the @code{logand}, just testing\n"
1452 "(logtest #b0100 #b1011) @result{} #f\n"
1453 "(logtest #b0100 #b0111) @result{} #t\n"
1455 #define FUNC_NAME s_scm_logtest
1459 if (SCM_I_INUMP (j
))
1461 nj
= SCM_I_INUM (j
);
1462 if (SCM_I_INUMP (k
))
1464 long nk
= SCM_I_INUM (k
);
1465 return scm_from_bool (nj
& nk
);
1467 else if (SCM_BIGP (k
))
1475 mpz_init_set_si (nj_z
, nj
);
1476 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1477 scm_remember_upto_here_1 (k
);
1478 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1484 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1486 else if (SCM_BIGP (j
))
1488 if (SCM_I_INUMP (k
))
1491 nj
= SCM_I_INUM (j
);
1494 else if (SCM_BIGP (k
))
1498 mpz_init (result_z
);
1502 scm_remember_upto_here_2 (j
, k
);
1503 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1504 mpz_clear (result_z
);
1508 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1511 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1516 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1518 "Test whether bit number @var{index} in @var{j} is set.\n"
1519 "@var{index} starts from 0 for the least significant bit.\n"
1522 "(logbit? 0 #b1101) @result{} #t\n"
1523 "(logbit? 1 #b1101) @result{} #f\n"
1524 "(logbit? 2 #b1101) @result{} #t\n"
1525 "(logbit? 3 #b1101) @result{} #t\n"
1526 "(logbit? 4 #b1101) @result{} #f\n"
1528 #define FUNC_NAME s_scm_logbit_p
1530 unsigned long int iindex
;
1531 iindex
= scm_to_ulong (index
);
1533 if (SCM_I_INUMP (j
))
1535 /* bits above what's in an inum follow the sign bit */
1536 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1537 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1539 else if (SCM_BIGP (j
))
1541 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1542 scm_remember_upto_here_1 (j
);
1543 return scm_from_bool (val
);
1546 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1551 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1553 "Return the integer which is the ones-complement of the integer\n"
1557 "(number->string (lognot #b10000000) 2)\n"
1558 " @result{} \"-10000001\"\n"
1559 "(number->string (lognot #b0) 2)\n"
1560 " @result{} \"-1\"\n"
1562 #define FUNC_NAME s_scm_lognot
1564 if (SCM_I_INUMP (n
)) {
1565 /* No overflow here, just need to toggle all the bits making up the inum.
1566 Enhancement: No need to strip the tag and add it back, could just xor
1567 a block of 1 bits, if that worked with the various debug versions of
1569 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1571 } else if (SCM_BIGP (n
)) {
1572 SCM result
= scm_i_mkbig ();
1573 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1574 scm_remember_upto_here_1 (n
);
1578 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1583 /* returns 0 if IN is not an integer. OUT must already be
1586 coerce_to_big (SCM in
, mpz_t out
)
1589 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1590 else if (SCM_I_INUMP (in
))
1591 mpz_set_si (out
, SCM_I_INUM (in
));
1598 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1599 (SCM n
, SCM k
, SCM m
),
1600 "Return @var{n} raised to the integer exponent\n"
1601 "@var{k}, modulo @var{m}.\n"
1604 "(modulo-expt 2 3 5)\n"
1607 #define FUNC_NAME s_scm_modulo_expt
1613 /* There are two classes of error we might encounter --
1614 1) Math errors, which we'll report by calling scm_num_overflow,
1616 2) wrong-type errors, which of course we'll report by calling
1618 We don't report those errors immediately, however; instead we do
1619 some cleanup first. These variables tell us which error (if
1620 any) we should report after cleaning up.
1622 int report_overflow
= 0;
1624 int position_of_wrong_type
= 0;
1625 SCM value_of_wrong_type
= SCM_INUM0
;
1627 SCM result
= SCM_UNDEFINED
;
1633 if (scm_is_eq (m
, SCM_INUM0
))
1635 report_overflow
= 1;
1639 if (!coerce_to_big (n
, n_tmp
))
1641 value_of_wrong_type
= n
;
1642 position_of_wrong_type
= 1;
1646 if (!coerce_to_big (k
, k_tmp
))
1648 value_of_wrong_type
= k
;
1649 position_of_wrong_type
= 2;
1653 if (!coerce_to_big (m
, m_tmp
))
1655 value_of_wrong_type
= m
;
1656 position_of_wrong_type
= 3;
1660 /* if the exponent K is negative, and we simply call mpz_powm, we
1661 will get a divide-by-zero exception when an inverse 1/n mod m
1662 doesn't exist (or is not unique). Since exceptions are hard to
1663 handle, we'll attempt the inversion "by hand" -- that way, we get
1664 a simple failure code, which is easy to handle. */
1666 if (-1 == mpz_sgn (k_tmp
))
1668 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1670 report_overflow
= 1;
1673 mpz_neg (k_tmp
, k_tmp
);
1676 result
= scm_i_mkbig ();
1677 mpz_powm (SCM_I_BIG_MPZ (result
),
1682 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1683 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1690 if (report_overflow
)
1691 scm_num_overflow (FUNC_NAME
);
1693 if (position_of_wrong_type
)
1694 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1695 value_of_wrong_type
);
1697 return scm_i_normbig (result
);
1701 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1703 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1704 "exact integer, @var{n} can be any number.\n"
1706 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1707 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1708 "includes @math{0^0} is 1.\n"
1711 "(integer-expt 2 5) @result{} 32\n"
1712 "(integer-expt -3 3) @result{} -27\n"
1713 "(integer-expt 5 -3) @result{} 1/125\n"
1714 "(integer-expt 0 0) @result{} 1\n"
1716 #define FUNC_NAME s_scm_integer_expt
1719 SCM z_i2
= SCM_BOOL_F
;
1721 SCM acc
= SCM_I_MAKINUM (1L);
1723 /* 0^0 == 1 according to R5RS */
1724 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1725 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1726 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1727 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1729 if (SCM_I_INUMP (k
))
1730 i2
= SCM_I_INUM (k
);
1731 else if (SCM_BIGP (k
))
1733 z_i2
= scm_i_clonebig (k
, 1);
1734 scm_remember_upto_here_1 (k
);
1738 SCM_WRONG_TYPE_ARG (2, k
);
1742 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1744 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1745 n
= scm_divide (n
, SCM_UNDEFINED
);
1749 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1753 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1755 return scm_product (acc
, n
);
1757 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1758 acc
= scm_product (acc
, n
);
1759 n
= scm_product (n
, n
);
1760 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1768 n
= scm_divide (n
, SCM_UNDEFINED
);
1775 return scm_product (acc
, n
);
1777 acc
= scm_product (acc
, n
);
1778 n
= scm_product (n
, n
);
1785 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1787 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1788 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1790 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1791 "@var{cnt} is negative it's a division, rounded towards negative\n"
1792 "infinity. (Note that this is not the same rounding as\n"
1793 "@code{quotient} does.)\n"
1795 "With @var{n} viewed as an infinite precision twos complement,\n"
1796 "@code{ash} means a left shift introducing zero bits, or a right\n"
1797 "shift dropping bits.\n"
1800 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1801 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1803 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1804 "(ash -23 -2) @result{} -6\n"
1806 #define FUNC_NAME s_scm_ash
1809 bits_to_shift
= scm_to_long (cnt
);
1811 if (SCM_I_INUMP (n
))
1813 long nn
= SCM_I_INUM (n
);
1815 if (bits_to_shift
> 0)
1817 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1818 overflow a non-zero fixnum. For smaller shifts we check the
1819 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1820 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1821 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1827 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1829 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1832 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1836 SCM result
= scm_i_long2big (nn
);
1837 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1844 bits_to_shift
= -bits_to_shift
;
1845 if (bits_to_shift
>= SCM_LONG_BIT
)
1846 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1848 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1852 else if (SCM_BIGP (n
))
1856 if (bits_to_shift
== 0)
1859 result
= scm_i_mkbig ();
1860 if (bits_to_shift
>= 0)
1862 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1868 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1869 we have to allocate a bignum even if the result is going to be a
1871 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1873 return scm_i_normbig (result
);
1879 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1885 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1886 (SCM n
, SCM start
, SCM end
),
1887 "Return the integer composed of the @var{start} (inclusive)\n"
1888 "through @var{end} (exclusive) bits of @var{n}. The\n"
1889 "@var{start}th bit becomes the 0-th bit in the result.\n"
1892 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1893 " @result{} \"1010\"\n"
1894 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1895 " @result{} \"10110\"\n"
1897 #define FUNC_NAME s_scm_bit_extract
1899 unsigned long int istart
, iend
, bits
;
1900 istart
= scm_to_ulong (start
);
1901 iend
= scm_to_ulong (end
);
1902 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1904 /* how many bits to keep */
1905 bits
= iend
- istart
;
1907 if (SCM_I_INUMP (n
))
1909 long int in
= SCM_I_INUM (n
);
1911 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1912 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1913 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1915 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1917 /* Since we emulate two's complement encoded numbers, this
1918 * special case requires us to produce a result that has
1919 * more bits than can be stored in a fixnum.
1921 SCM result
= scm_i_long2big (in
);
1922 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1927 /* mask down to requisite bits */
1928 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1929 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1931 else if (SCM_BIGP (n
))
1936 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1940 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1941 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1942 such bits into a ulong. */
1943 result
= scm_i_mkbig ();
1944 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1945 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1946 result
= scm_i_normbig (result
);
1948 scm_remember_upto_here_1 (n
);
1952 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1957 static const char scm_logtab
[] = {
1958 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1961 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1963 "Return the number of bits in integer @var{n}. If integer is\n"
1964 "positive, the 1-bits in its binary representation are counted.\n"
1965 "If negative, the 0-bits in its two's-complement binary\n"
1966 "representation are counted. If 0, 0 is returned.\n"
1969 "(logcount #b10101010)\n"
1976 #define FUNC_NAME s_scm_logcount
1978 if (SCM_I_INUMP (n
))
1980 unsigned long int c
= 0;
1981 long int nn
= SCM_I_INUM (n
);
1986 c
+= scm_logtab
[15 & nn
];
1989 return SCM_I_MAKINUM (c
);
1991 else if (SCM_BIGP (n
))
1993 unsigned long count
;
1994 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1995 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1997 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1998 scm_remember_upto_here_1 (n
);
1999 return SCM_I_MAKINUM (count
);
2002 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2007 static const char scm_ilentab
[] = {
2008 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2012 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2014 "Return the number of bits necessary to represent @var{n}.\n"
2017 "(integer-length #b10101010)\n"
2019 "(integer-length 0)\n"
2021 "(integer-length #b1111)\n"
2024 #define FUNC_NAME s_scm_integer_length
2026 if (SCM_I_INUMP (n
))
2028 unsigned long int c
= 0;
2030 long int nn
= SCM_I_INUM (n
);
2036 l
= scm_ilentab
[15 & nn
];
2039 return SCM_I_MAKINUM (c
- 4 + l
);
2041 else if (SCM_BIGP (n
))
2043 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2044 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2045 1 too big, so check for that and adjust. */
2046 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2047 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2048 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2049 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2051 scm_remember_upto_here_1 (n
);
2052 return SCM_I_MAKINUM (size
);
2055 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2059 /*** NUMBERS -> STRINGS ***/
2060 #define SCM_MAX_DBL_PREC 60
2061 #define SCM_MAX_DBL_RADIX 36
2063 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2064 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2065 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2068 void init_dblprec(int *prec
, int radix
) {
2069 /* determine floating point precision by adding successively
2070 smaller increments to 1.0 until it is considered == 1.0 */
2071 double f
= ((double)1.0)/radix
;
2072 double fsum
= 1.0 + f
;
2077 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2089 void init_fx_radix(double *fx_list
, int radix
)
2091 /* initialize a per-radix list of tolerances. When added
2092 to a number < 1.0, we can determine if we should raund
2093 up and quit converting a number to a string. */
2097 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2098 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2101 /* use this array as a way to generate a single digit */
2102 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2105 idbl2str (double f
, char *a
, int radix
)
2107 int efmt
, dpt
, d
, i
, wp
;
2109 #ifdef DBL_MIN_10_EXP
2112 #endif /* DBL_MIN_10_EXP */
2117 radix
> SCM_MAX_DBL_RADIX
)
2119 /* revert to existing behavior */
2123 wp
= scm_dblprec
[radix
-2];
2124 fx
= fx_per_radix
[radix
-2];
2128 #ifdef HAVE_COPYSIGN
2129 double sgn
= copysign (1.0, f
);
2134 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2140 strcpy (a
, "-inf.0");
2142 strcpy (a
, "+inf.0");
2145 else if (xisnan (f
))
2147 strcpy (a
, "+nan.0");
2157 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2158 make-uniform-vector, from causing infinite loops. */
2159 /* just do the checking...if it passes, we do the conversion for our
2160 radix again below */
2167 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2175 while (f_cpy
> 10.0)
2178 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2199 if (f
+ fx
[wp
] >= radix
)
2206 /* adding 9999 makes this equivalent to abs(x) % 3 */
2207 dpt
= (exp
+ 9999) % 3;
2211 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2233 a
[ch
++] = number_chars
[d
];
2236 if (f
+ fx
[wp
] >= 1.0)
2238 a
[ch
- 1] = number_chars
[d
+1];
2250 if ((dpt
> 4) && (exp
> 6))
2252 d
= (a
[0] == '-' ? 2 : 1);
2253 for (i
= ch
++; i
> d
; i
--)
2266 if (a
[ch
- 1] == '.')
2267 a
[ch
++] = '0'; /* trailing zero */
2276 for (i
= radix
; i
<= exp
; i
*= radix
);
2277 for (i
/= radix
; i
; i
/= radix
)
2279 a
[ch
++] = number_chars
[exp
/ i
];
2288 icmplx2str (double real
, double imag
, char *str
, int radix
)
2292 i
= idbl2str (real
, str
, radix
);
2295 /* Don't output a '+' for negative numbers or for Inf and
2296 NaN. They will provide their own sign. */
2297 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2299 i
+= idbl2str (imag
, &str
[i
], radix
);
2306 iflo2str (SCM flt
, char *str
, int radix
)
2309 if (SCM_REALP (flt
))
2310 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2312 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2317 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2318 characters in the result.
2320 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2322 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2327 return scm_iuint2str (-num
, rad
, p
) + 1;
2330 return scm_iuint2str (num
, rad
, p
);
2333 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2334 characters in the result.
2336 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2338 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2342 scm_t_uintmax n
= num
;
2344 for (n
/= rad
; n
> 0; n
/= rad
)
2354 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2359 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2361 "Return a string holding the external representation of the\n"
2362 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2363 "inexact, a radix of 10 will be used.")
2364 #define FUNC_NAME s_scm_number_to_string
2368 if (SCM_UNBNDP (radix
))
2371 base
= scm_to_signed_integer (radix
, 2, 36);
2373 if (SCM_I_INUMP (n
))
2375 char num_buf
[SCM_INTBUFLEN
];
2376 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2377 return scm_from_locale_stringn (num_buf
, length
);
2379 else if (SCM_BIGP (n
))
2381 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2382 scm_remember_upto_here_1 (n
);
2383 return scm_take_locale_string (str
);
2385 else if (SCM_FRACTIONP (n
))
2387 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2388 scm_from_locale_string ("/"),
2389 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2391 else if (SCM_INEXACTP (n
))
2393 char num_buf
[FLOBUFLEN
];
2394 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2397 SCM_WRONG_TYPE_ARG (1, n
);
2402 /* These print routines used to be stubbed here so that scm_repl.c
2403 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2406 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2408 char num_buf
[FLOBUFLEN
];
2409 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2414 scm_i_print_double (double val
, SCM port
)
2416 char num_buf
[FLOBUFLEN
];
2417 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2421 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2424 char num_buf
[FLOBUFLEN
];
2425 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2430 scm_i_print_complex (double real
, double imag
, SCM port
)
2432 char num_buf
[FLOBUFLEN
];
2433 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2437 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2440 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2441 scm_lfwrite_str (str
, port
);
2442 scm_remember_upto_here_1 (str
);
2447 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2449 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2450 scm_remember_upto_here_1 (exp
);
2451 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2455 /*** END nums->strs ***/
2458 /*** STRINGS -> NUMBERS ***/
2460 /* The following functions implement the conversion from strings to numbers.
2461 * The implementation somehow follows the grammar for numbers as it is given
2462 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2463 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2464 * points should be noted about the implementation:
2465 * * Each function keeps a local index variable 'idx' that points at the
2466 * current position within the parsed string. The global index is only
2467 * updated if the function could parse the corresponding syntactic unit
2469 * * Similarly, the functions keep track of indicators of inexactness ('#',
2470 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2471 * global exactness information is only updated after each part has been
2472 * successfully parsed.
2473 * * Sequences of digits are parsed into temporary variables holding fixnums.
2474 * Only if these fixnums would overflow, the result variables are updated
2475 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2476 * the temporary variables holding the fixnums are cleared, and the process
2477 * starts over again. If for example fixnums were able to store five decimal
2478 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2479 * and the result was computed as 12345 * 100000 + 67890. In other words,
2480 * only every five digits two bignum operations were performed.
2483 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2485 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2487 /* In non ASCII-style encodings the following macro might not work. */
2488 #define XDIGIT2UINT(d) \
2489 (uc_is_property_decimal_digit ((int) (unsigned char) d) \
2491 : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
2494 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2495 unsigned int radix
, enum t_exactness
*p_exactness
)
2497 unsigned int idx
= *p_idx
;
2498 unsigned int hash_seen
= 0;
2499 scm_t_bits shift
= 1;
2501 unsigned int digit_value
;
2504 size_t len
= scm_i_string_length (mem
);
2509 c
= scm_i_string_ref (mem
, idx
);
2510 if (!uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2512 digit_value
= XDIGIT2UINT (c
);
2513 if (digit_value
>= radix
)
2517 result
= SCM_I_MAKINUM (digit_value
);
2520 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2521 if (uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2525 digit_value
= XDIGIT2UINT (c
);
2526 if (digit_value
>= radix
)
2538 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2540 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2542 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2549 shift
= shift
* radix
;
2550 add
= add
* radix
+ digit_value
;
2555 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2557 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2561 *p_exactness
= INEXACT
;
2567 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2568 * covers the parts of the rules that start at a potential point. The value
2569 * of the digits up to the point have been parsed by the caller and are given
2570 * in variable result. The content of *p_exactness indicates, whether a hash
2571 * has already been seen in the digits before the point.
2574 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2577 mem2decimal_from_point (SCM result
, SCM mem
,
2578 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2580 unsigned int idx
= *p_idx
;
2581 enum t_exactness x
= *p_exactness
;
2582 size_t len
= scm_i_string_length (mem
);
2587 if (scm_i_string_ref (mem
, idx
) == '.')
2589 scm_t_bits shift
= 1;
2591 unsigned int digit_value
;
2592 SCM big_shift
= SCM_I_MAKINUM (1);
2597 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2598 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2603 digit_value
= DIGIT2UINT (c
);
2614 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2616 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2617 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2619 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2627 add
= add
* 10 + digit_value
;
2633 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2634 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2635 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2638 result
= scm_divide (result
, big_shift
);
2640 /* We've seen a decimal point, thus the value is implicitly inexact. */
2652 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2654 switch (scm_i_string_ref (mem
, idx
))
2666 c
= scm_i_string_ref (mem
, idx
);
2674 c
= scm_i_string_ref (mem
, idx
);
2683 c
= scm_i_string_ref (mem
, idx
);
2688 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2692 exponent
= DIGIT2UINT (c
);
2695 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2696 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2699 if (exponent
<= SCM_MAXEXP
)
2700 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2706 if (exponent
> SCM_MAXEXP
)
2708 size_t exp_len
= idx
- start
;
2709 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2710 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2711 scm_out_of_range ("string->number", exp_num
);
2714 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2716 result
= scm_product (result
, e
);
2718 result
= scm_divide2real (result
, e
);
2720 /* We've seen an exponent, thus the value is implicitly inexact. */
2738 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2741 mem2ureal (SCM mem
, unsigned int *p_idx
,
2742 unsigned int radix
, enum t_exactness
*p_exactness
)
2744 unsigned int idx
= *p_idx
;
2746 size_t len
= scm_i_string_length (mem
);
2748 /* Start off believing that the number will be exact. This changes
2749 to INEXACT if we see a decimal point or a hash. */
2750 enum t_exactness x
= EXACT
;
2755 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2761 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2763 /* Cobble up the fractional part. We might want to set the
2764 NaN's mantissa from it. */
2766 mem2uinteger (mem
, &idx
, 10, &x
);
2771 if (scm_i_string_ref (mem
, idx
) == '.')
2775 else if (idx
+ 1 == len
)
2777 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2780 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2787 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2788 if (scm_is_false (uinteger
))
2793 else if (scm_i_string_ref (mem
, idx
) == '/')
2801 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2802 if (scm_is_false (divisor
))
2805 /* both are int/big here, I assume */
2806 result
= scm_i_make_ratio (uinteger
, divisor
);
2808 else if (radix
== 10)
2810 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2811 if (scm_is_false (result
))
2820 /* Update *p_exactness if the number just read was inexact. This is
2821 important for complex numbers, so that a complex number is
2822 treated as inexact overall if either its real or imaginary part
2828 /* When returning an inexact zero, make sure it is represented as a
2829 floating point value so that we can change its sign.
2831 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2832 result
= scm_from_double (0.0);
2838 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2841 mem2complex (SCM mem
, unsigned int idx
,
2842 unsigned int radix
, enum t_exactness
*p_exactness
)
2847 size_t len
= scm_i_string_length (mem
);
2852 c
= scm_i_string_ref (mem
, idx
);
2867 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2868 if (scm_is_false (ureal
))
2870 /* input must be either +i or -i */
2875 if (scm_i_string_ref (mem
, idx
) == 'i'
2876 || scm_i_string_ref (mem
, idx
) == 'I')
2882 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2889 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2890 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2895 c
= scm_i_string_ref (mem
, idx
);
2899 /* either +<ureal>i or -<ureal>i */
2906 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2909 /* polar input: <real>@<real>. */
2920 c
= scm_i_string_ref (mem
, idx
);
2938 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2939 if (scm_is_false (angle
))
2944 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2945 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2947 result
= scm_make_polar (ureal
, angle
);
2952 /* expecting input matching <real>[+-]<ureal>?i */
2959 int sign
= (c
== '+') ? 1 : -1;
2960 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2962 if (scm_is_false (imag
))
2963 imag
= SCM_I_MAKINUM (sign
);
2964 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2965 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2969 if (scm_i_string_ref (mem
, idx
) != 'i'
2970 && scm_i_string_ref (mem
, idx
) != 'I')
2977 return scm_make_rectangular (ureal
, imag
);
2986 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2988 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2991 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
2993 unsigned int idx
= 0;
2994 unsigned int radix
= NO_RADIX
;
2995 enum t_exactness forced_x
= NO_EXACTNESS
;
2996 enum t_exactness implicit_x
= EXACT
;
2998 size_t len
= scm_i_string_length (mem
);
3000 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3001 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3003 switch (scm_i_string_ref (mem
, idx
+ 1))
3006 if (radix
!= NO_RADIX
)
3011 if (radix
!= NO_RADIX
)
3016 if (forced_x
!= NO_EXACTNESS
)
3021 if (forced_x
!= NO_EXACTNESS
)
3026 if (radix
!= NO_RADIX
)
3031 if (radix
!= NO_RADIX
)
3041 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3042 if (radix
== NO_RADIX
)
3043 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3045 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3047 if (scm_is_false (result
))
3053 if (SCM_INEXACTP (result
))
3054 return scm_inexact_to_exact (result
);
3058 if (SCM_INEXACTP (result
))
3061 return scm_exact_to_inexact (result
);
3064 if (implicit_x
== INEXACT
)
3066 if (SCM_INEXACTP (result
))
3069 return scm_exact_to_inexact (result
);
3077 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3078 unsigned int default_radix
)
3080 SCM str
= scm_from_locale_stringn (mem
, len
);
3082 return scm_i_string_to_number (str
, default_radix
);
3086 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3087 (SCM string
, SCM radix
),
3088 "Return a number of the maximally precise representation\n"
3089 "expressed by the given @var{string}. @var{radix} must be an\n"
3090 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3091 "is a default radix that may be overridden by an explicit radix\n"
3092 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3093 "supplied, then the default radix is 10. If string is not a\n"
3094 "syntactically valid notation for a number, then\n"
3095 "@code{string->number} returns @code{#f}.")
3096 #define FUNC_NAME s_scm_string_to_number
3100 SCM_VALIDATE_STRING (1, string
);
3102 if (SCM_UNBNDP (radix
))
3105 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3107 answer
= scm_i_string_to_number (string
, base
);
3108 scm_remember_upto_here_1 (string
);
3114 /*** END strs->nums ***/
3118 scm_bigequal (SCM x
, SCM y
)
3120 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3121 scm_remember_upto_here_2 (x
, y
);
3122 return scm_from_bool (0 == result
);
3126 scm_real_equalp (SCM x
, SCM y
)
3128 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3132 scm_complex_equalp (SCM x
, SCM y
)
3134 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3135 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3139 scm_i_fraction_equalp (SCM x
, SCM y
)
3141 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3142 SCM_FRACTION_NUMERATOR (y
)))
3143 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3144 SCM_FRACTION_DENOMINATOR (y
))))
3151 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3153 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3155 #define FUNC_NAME s_scm_number_p
3157 return scm_from_bool (SCM_NUMBERP (x
));
3161 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3163 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3164 "otherwise. Note that the sets of real, rational and integer\n"
3165 "values form subsets of the set of complex numbers, i. e. the\n"
3166 "predicate will also be fulfilled if @var{x} is a real,\n"
3167 "rational or integer number.")
3168 #define FUNC_NAME s_scm_complex_p
3170 /* all numbers are complex. */
3171 return scm_number_p (x
);
3175 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3177 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3178 "otherwise. Note that the set of integer values forms a subset of\n"
3179 "the set of real numbers, i. e. the predicate will also be\n"
3180 "fulfilled if @var{x} is an integer number.")
3181 #define FUNC_NAME s_scm_real_p
3183 /* we can't represent irrational numbers. */
3184 return scm_rational_p (x
);
3188 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3190 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3191 "otherwise. Note that the set of integer values forms a subset of\n"
3192 "the set of rational numbers, i. e. the predicate will also be\n"
3193 "fulfilled if @var{x} is an integer number.")
3194 #define FUNC_NAME s_scm_rational_p
3196 if (SCM_I_INUMP (x
))
3198 else if (SCM_IMP (x
))
3200 else if (SCM_BIGP (x
))
3202 else if (SCM_FRACTIONP (x
))
3204 else if (SCM_REALP (x
))
3205 /* due to their limited precision, all floating point numbers are
3206 rational as well. */
3213 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3215 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3217 #define FUNC_NAME s_scm_integer_p
3220 if (SCM_I_INUMP (x
))
3226 if (!SCM_INEXACTP (x
))
3228 if (SCM_COMPLEXP (x
))
3230 r
= SCM_REAL_VALUE (x
);
3231 /* +/-inf passes r==floor(r), making those #t */
3239 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3241 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3243 #define FUNC_NAME s_scm_inexact_p
3245 if (SCM_INEXACTP (x
))
3247 if (SCM_NUMBERP (x
))
3249 SCM_WRONG_TYPE_ARG (1, x
);
3254 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3255 /* "Return @code{#t} if all parameters are numerically equal." */
3257 scm_num_eq_p (SCM x
, SCM y
)
3260 if (SCM_I_INUMP (x
))
3262 long xx
= SCM_I_INUM (x
);
3263 if (SCM_I_INUMP (y
))
3265 long yy
= SCM_I_INUM (y
);
3266 return scm_from_bool (xx
== yy
);
3268 else if (SCM_BIGP (y
))
3270 else if (SCM_REALP (y
))
3272 /* On a 32-bit system an inum fits a double, we can cast the inum
3273 to a double and compare.
3275 But on a 64-bit system an inum is bigger than a double and
3276 casting it to a double (call that dxx) will round. dxx is at
3277 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3278 an integer and fits a long. So we cast yy to a long and
3279 compare with plain xx.
3281 An alternative (for any size system actually) would be to check
3282 yy is an integer (with floor) and is in range of an inum
3283 (compare against appropriate powers of 2) then test
3284 xx==(long)yy. It's just a matter of which casts/comparisons
3285 might be fastest or easiest for the cpu. */
3287 double yy
= SCM_REAL_VALUE (y
);
3288 return scm_from_bool ((double) xx
== yy
3289 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3290 || xx
== (long) yy
));
3292 else if (SCM_COMPLEXP (y
))
3293 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3294 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3295 else if (SCM_FRACTIONP (y
))
3298 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3300 else if (SCM_BIGP (x
))
3302 if (SCM_I_INUMP (y
))
3304 else if (SCM_BIGP (y
))
3306 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3307 scm_remember_upto_here_2 (x
, y
);
3308 return scm_from_bool (0 == cmp
);
3310 else if (SCM_REALP (y
))
3313 if (xisnan (SCM_REAL_VALUE (y
)))
3315 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3316 scm_remember_upto_here_1 (x
);
3317 return scm_from_bool (0 == cmp
);
3319 else if (SCM_COMPLEXP (y
))
3322 if (0.0 != SCM_COMPLEX_IMAG (y
))
3324 if (xisnan (SCM_COMPLEX_REAL (y
)))
3326 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3327 scm_remember_upto_here_1 (x
);
3328 return scm_from_bool (0 == cmp
);
3330 else if (SCM_FRACTIONP (y
))
3333 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3335 else if (SCM_REALP (x
))
3337 double xx
= SCM_REAL_VALUE (x
);
3338 if (SCM_I_INUMP (y
))
3340 /* see comments with inum/real above */
3341 long yy
= SCM_I_INUM (y
);
3342 return scm_from_bool (xx
== (double) yy
3343 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3344 || (long) xx
== yy
));
3346 else if (SCM_BIGP (y
))
3349 if (xisnan (SCM_REAL_VALUE (x
)))
3351 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3352 scm_remember_upto_here_1 (y
);
3353 return scm_from_bool (0 == cmp
);
3355 else if (SCM_REALP (y
))
3356 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3357 else if (SCM_COMPLEXP (y
))
3358 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3359 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3360 else if (SCM_FRACTIONP (y
))
3362 double xx
= SCM_REAL_VALUE (x
);
3366 return scm_from_bool (xx
< 0.0);
3367 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3371 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3373 else if (SCM_COMPLEXP (x
))
3375 if (SCM_I_INUMP (y
))
3376 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3377 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3378 else if (SCM_BIGP (y
))
3381 if (0.0 != SCM_COMPLEX_IMAG (x
))
3383 if (xisnan (SCM_COMPLEX_REAL (x
)))
3385 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3386 scm_remember_upto_here_1 (y
);
3387 return scm_from_bool (0 == cmp
);
3389 else if (SCM_REALP (y
))
3390 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3391 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3392 else if (SCM_COMPLEXP (y
))
3393 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3394 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3395 else if (SCM_FRACTIONP (y
))
3398 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3400 xx
= SCM_COMPLEX_REAL (x
);
3404 return scm_from_bool (xx
< 0.0);
3405 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3409 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3411 else if (SCM_FRACTIONP (x
))
3413 if (SCM_I_INUMP (y
))
3415 else if (SCM_BIGP (y
))
3417 else if (SCM_REALP (y
))
3419 double yy
= SCM_REAL_VALUE (y
);
3423 return scm_from_bool (0.0 < yy
);
3424 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3427 else if (SCM_COMPLEXP (y
))
3430 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3432 yy
= SCM_COMPLEX_REAL (y
);
3436 return scm_from_bool (0.0 < yy
);
3437 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3440 else if (SCM_FRACTIONP (y
))
3441 return scm_i_fraction_equalp (x
, y
);
3443 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3446 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3450 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3451 done are good for inums, but for bignums an answer can almost always be
3452 had by just examining a few high bits of the operands, as done by GMP in
3453 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3454 of the float exponent to take into account. */
3456 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3457 /* "Return @code{#t} if the list of parameters is monotonically\n"
3461 scm_less_p (SCM x
, SCM y
)
3464 if (SCM_I_INUMP (x
))
3466 long xx
= SCM_I_INUM (x
);
3467 if (SCM_I_INUMP (y
))
3469 long yy
= SCM_I_INUM (y
);
3470 return scm_from_bool (xx
< yy
);
3472 else if (SCM_BIGP (y
))
3474 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3475 scm_remember_upto_here_1 (y
);
3476 return scm_from_bool (sgn
> 0);
3478 else if (SCM_REALP (y
))
3479 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3480 else if (SCM_FRACTIONP (y
))
3482 /* "x < a/b" becomes "x*b < a" */
3484 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3485 y
= SCM_FRACTION_NUMERATOR (y
);
3489 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3491 else if (SCM_BIGP (x
))
3493 if (SCM_I_INUMP (y
))
3495 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3496 scm_remember_upto_here_1 (x
);
3497 return scm_from_bool (sgn
< 0);
3499 else if (SCM_BIGP (y
))
3501 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3502 scm_remember_upto_here_2 (x
, y
);
3503 return scm_from_bool (cmp
< 0);
3505 else if (SCM_REALP (y
))
3508 if (xisnan (SCM_REAL_VALUE (y
)))
3510 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3511 scm_remember_upto_here_1 (x
);
3512 return scm_from_bool (cmp
< 0);
3514 else if (SCM_FRACTIONP (y
))
3517 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3519 else if (SCM_REALP (x
))
3521 if (SCM_I_INUMP (y
))
3522 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3523 else if (SCM_BIGP (y
))
3526 if (xisnan (SCM_REAL_VALUE (x
)))
3528 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3529 scm_remember_upto_here_1 (y
);
3530 return scm_from_bool (cmp
> 0);
3532 else if (SCM_REALP (y
))
3533 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3534 else if (SCM_FRACTIONP (y
))
3536 double xx
= SCM_REAL_VALUE (x
);
3540 return scm_from_bool (xx
< 0.0);
3541 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3545 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3547 else if (SCM_FRACTIONP (x
))
3549 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3551 /* "a/b < y" becomes "a < y*b" */
3552 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3553 x
= SCM_FRACTION_NUMERATOR (x
);
3556 else if (SCM_REALP (y
))
3558 double yy
= SCM_REAL_VALUE (y
);
3562 return scm_from_bool (0.0 < yy
);
3563 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3566 else if (SCM_FRACTIONP (y
))
3568 /* "a/b < c/d" becomes "a*d < c*b" */
3569 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3570 SCM_FRACTION_DENOMINATOR (y
));
3571 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3572 SCM_FRACTION_DENOMINATOR (x
));
3578 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3581 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3585 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3586 /* "Return @code{#t} if the list of parameters is monotonically\n"
3589 #define FUNC_NAME s_scm_gr_p
3591 scm_gr_p (SCM x
, SCM y
)
3593 if (!SCM_NUMBERP (x
))
3594 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3595 else if (!SCM_NUMBERP (y
))
3596 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3598 return scm_less_p (y
, x
);
3603 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3604 /* "Return @code{#t} if the list of parameters is monotonically\n"
3607 #define FUNC_NAME s_scm_leq_p
3609 scm_leq_p (SCM x
, SCM y
)
3611 if (!SCM_NUMBERP (x
))
3612 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3613 else if (!SCM_NUMBERP (y
))
3614 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3615 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3618 return scm_not (scm_less_p (y
, x
));
3623 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3624 /* "Return @code{#t} if the list of parameters is monotonically\n"
3627 #define FUNC_NAME s_scm_geq_p
3629 scm_geq_p (SCM x
, SCM y
)
3631 if (!SCM_NUMBERP (x
))
3632 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3633 else if (!SCM_NUMBERP (y
))
3634 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3635 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3638 return scm_not (scm_less_p (x
, y
));
3643 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3644 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3650 if (SCM_I_INUMP (z
))
3651 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3652 else if (SCM_BIGP (z
))
3654 else if (SCM_REALP (z
))
3655 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3656 else if (SCM_COMPLEXP (z
))
3657 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3658 && SCM_COMPLEX_IMAG (z
) == 0.0);
3659 else if (SCM_FRACTIONP (z
))
3662 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3666 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3667 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3671 scm_positive_p (SCM x
)
3673 if (SCM_I_INUMP (x
))
3674 return scm_from_bool (SCM_I_INUM (x
) > 0);
3675 else if (SCM_BIGP (x
))
3677 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3678 scm_remember_upto_here_1 (x
);
3679 return scm_from_bool (sgn
> 0);
3681 else if (SCM_REALP (x
))
3682 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3683 else if (SCM_FRACTIONP (x
))
3684 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3686 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3690 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3691 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3695 scm_negative_p (SCM x
)
3697 if (SCM_I_INUMP (x
))
3698 return scm_from_bool (SCM_I_INUM (x
) < 0);
3699 else if (SCM_BIGP (x
))
3701 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3702 scm_remember_upto_here_1 (x
);
3703 return scm_from_bool (sgn
< 0);
3705 else if (SCM_REALP (x
))
3706 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3707 else if (SCM_FRACTIONP (x
))
3708 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3710 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3714 /* scm_min and scm_max return an inexact when either argument is inexact, as
3715 required by r5rs. On that basis, for exact/inexact combinations the
3716 exact is converted to inexact to compare and possibly return. This is
3717 unlike scm_less_p above which takes some trouble to preserve all bits in
3718 its test, such trouble is not required for min and max. */
3720 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3721 /* "Return the maximum of all parameter values."
3724 scm_max (SCM x
, SCM y
)
3729 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3730 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3733 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3736 if (SCM_I_INUMP (x
))
3738 long xx
= SCM_I_INUM (x
);
3739 if (SCM_I_INUMP (y
))
3741 long yy
= SCM_I_INUM (y
);
3742 return (xx
< yy
) ? y
: x
;
3744 else if (SCM_BIGP (y
))
3746 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3747 scm_remember_upto_here_1 (y
);
3748 return (sgn
< 0) ? x
: y
;
3750 else if (SCM_REALP (y
))
3753 /* if y==NaN then ">" is false and we return NaN */
3754 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3756 else if (SCM_FRACTIONP (y
))
3759 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3762 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3764 else if (SCM_BIGP (x
))
3766 if (SCM_I_INUMP (y
))
3768 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3769 scm_remember_upto_here_1 (x
);
3770 return (sgn
< 0) ? y
: x
;
3772 else if (SCM_BIGP (y
))
3774 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3775 scm_remember_upto_here_2 (x
, y
);
3776 return (cmp
> 0) ? x
: y
;
3778 else if (SCM_REALP (y
))
3780 /* if y==NaN then xx>yy is false, so we return the NaN y */
3783 xx
= scm_i_big2dbl (x
);
3784 yy
= SCM_REAL_VALUE (y
);
3785 return (xx
> yy
? scm_from_double (xx
) : y
);
3787 else if (SCM_FRACTIONP (y
))
3792 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3794 else if (SCM_REALP (x
))
3796 if (SCM_I_INUMP (y
))
3798 double z
= SCM_I_INUM (y
);
3799 /* if x==NaN then "<" is false and we return NaN */
3800 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3802 else if (SCM_BIGP (y
))
3807 else if (SCM_REALP (y
))
3809 /* if x==NaN then our explicit check means we return NaN
3810 if y==NaN then ">" is false and we return NaN
3811 calling isnan is unavoidable, since it's the only way to know
3812 which of x or y causes any compares to be false */
3813 double xx
= SCM_REAL_VALUE (x
);
3814 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3816 else if (SCM_FRACTIONP (y
))
3818 double yy
= scm_i_fraction2double (y
);
3819 double xx
= SCM_REAL_VALUE (x
);
3820 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3823 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3825 else if (SCM_FRACTIONP (x
))
3827 if (SCM_I_INUMP (y
))
3831 else if (SCM_BIGP (y
))
3835 else if (SCM_REALP (y
))
3837 double xx
= scm_i_fraction2double (x
);
3838 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3840 else if (SCM_FRACTIONP (y
))
3845 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3848 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3852 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3853 /* "Return the minium of all parameter values."
3856 scm_min (SCM x
, SCM y
)
3861 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3862 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3865 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3868 if (SCM_I_INUMP (x
))
3870 long xx
= SCM_I_INUM (x
);
3871 if (SCM_I_INUMP (y
))
3873 long yy
= SCM_I_INUM (y
);
3874 return (xx
< yy
) ? x
: y
;
3876 else if (SCM_BIGP (y
))
3878 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3879 scm_remember_upto_here_1 (y
);
3880 return (sgn
< 0) ? y
: x
;
3882 else if (SCM_REALP (y
))
3885 /* if y==NaN then "<" is false and we return NaN */
3886 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3888 else if (SCM_FRACTIONP (y
))
3891 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3894 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3896 else if (SCM_BIGP (x
))
3898 if (SCM_I_INUMP (y
))
3900 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3901 scm_remember_upto_here_1 (x
);
3902 return (sgn
< 0) ? x
: y
;
3904 else if (SCM_BIGP (y
))
3906 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3907 scm_remember_upto_here_2 (x
, y
);
3908 return (cmp
> 0) ? y
: x
;
3910 else if (SCM_REALP (y
))
3912 /* if y==NaN then xx<yy is false, so we return the NaN y */
3915 xx
= scm_i_big2dbl (x
);
3916 yy
= SCM_REAL_VALUE (y
);
3917 return (xx
< yy
? scm_from_double (xx
) : y
);
3919 else if (SCM_FRACTIONP (y
))
3924 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3926 else if (SCM_REALP (x
))
3928 if (SCM_I_INUMP (y
))
3930 double z
= SCM_I_INUM (y
);
3931 /* if x==NaN then "<" is false and we return NaN */
3932 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3934 else if (SCM_BIGP (y
))
3939 else if (SCM_REALP (y
))
3941 /* if x==NaN then our explicit check means we return NaN
3942 if y==NaN then "<" is false and we return NaN
3943 calling isnan is unavoidable, since it's the only way to know
3944 which of x or y causes any compares to be false */
3945 double xx
= SCM_REAL_VALUE (x
);
3946 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3948 else if (SCM_FRACTIONP (y
))
3950 double yy
= scm_i_fraction2double (y
);
3951 double xx
= SCM_REAL_VALUE (x
);
3952 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3955 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3957 else if (SCM_FRACTIONP (x
))
3959 if (SCM_I_INUMP (y
))
3963 else if (SCM_BIGP (y
))
3967 else if (SCM_REALP (y
))
3969 double xx
= scm_i_fraction2double (x
);
3970 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3972 else if (SCM_FRACTIONP (y
))
3977 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3980 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3984 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3985 /* "Return the sum of all parameter values. Return 0 if called without\n"
3989 scm_sum (SCM x
, SCM y
)
3991 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3993 if (SCM_NUMBERP (x
)) return x
;
3994 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3995 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3998 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4000 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4002 long xx
= SCM_I_INUM (x
);
4003 long yy
= SCM_I_INUM (y
);
4004 long int z
= xx
+ yy
;
4005 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4007 else if (SCM_BIGP (y
))
4012 else if (SCM_REALP (y
))
4014 long int xx
= SCM_I_INUM (x
);
4015 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4017 else if (SCM_COMPLEXP (y
))
4019 long int xx
= SCM_I_INUM (x
);
4020 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4021 SCM_COMPLEX_IMAG (y
));
4023 else if (SCM_FRACTIONP (y
))
4024 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4025 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4026 SCM_FRACTION_DENOMINATOR (y
));
4028 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4029 } else if (SCM_BIGP (x
))
4031 if (SCM_I_INUMP (y
))
4036 inum
= SCM_I_INUM (y
);
4039 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4042 SCM result
= scm_i_mkbig ();
4043 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4044 scm_remember_upto_here_1 (x
);
4045 /* we know the result will have to be a bignum */
4048 return scm_i_normbig (result
);
4052 SCM result
= scm_i_mkbig ();
4053 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4054 scm_remember_upto_here_1 (x
);
4055 /* we know the result will have to be a bignum */
4058 return scm_i_normbig (result
);
4061 else if (SCM_BIGP (y
))
4063 SCM result
= scm_i_mkbig ();
4064 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4065 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4066 mpz_add (SCM_I_BIG_MPZ (result
),
4069 scm_remember_upto_here_2 (x
, y
);
4070 /* we know the result will have to be a bignum */
4073 return scm_i_normbig (result
);
4075 else if (SCM_REALP (y
))
4077 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4078 scm_remember_upto_here_1 (x
);
4079 return scm_from_double (result
);
4081 else if (SCM_COMPLEXP (y
))
4083 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4084 + SCM_COMPLEX_REAL (y
));
4085 scm_remember_upto_here_1 (x
);
4086 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4088 else if (SCM_FRACTIONP (y
))
4089 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4090 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4091 SCM_FRACTION_DENOMINATOR (y
));
4093 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4095 else if (SCM_REALP (x
))
4097 if (SCM_I_INUMP (y
))
4098 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4099 else if (SCM_BIGP (y
))
4101 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4102 scm_remember_upto_here_1 (y
);
4103 return scm_from_double (result
);
4105 else if (SCM_REALP (y
))
4106 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4107 else if (SCM_COMPLEXP (y
))
4108 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4109 SCM_COMPLEX_IMAG (y
));
4110 else if (SCM_FRACTIONP (y
))
4111 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4113 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4115 else if (SCM_COMPLEXP (x
))
4117 if (SCM_I_INUMP (y
))
4118 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4119 SCM_COMPLEX_IMAG (x
));
4120 else if (SCM_BIGP (y
))
4122 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4123 + SCM_COMPLEX_REAL (x
));
4124 scm_remember_upto_here_1 (y
);
4125 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4127 else if (SCM_REALP (y
))
4128 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4129 SCM_COMPLEX_IMAG (x
));
4130 else if (SCM_COMPLEXP (y
))
4131 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4132 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4133 else if (SCM_FRACTIONP (y
))
4134 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4135 SCM_COMPLEX_IMAG (x
));
4137 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4139 else if (SCM_FRACTIONP (x
))
4141 if (SCM_I_INUMP (y
))
4142 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4143 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4144 SCM_FRACTION_DENOMINATOR (x
));
4145 else if (SCM_BIGP (y
))
4146 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4147 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4148 SCM_FRACTION_DENOMINATOR (x
));
4149 else if (SCM_REALP (y
))
4150 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4151 else if (SCM_COMPLEXP (y
))
4152 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4153 SCM_COMPLEX_IMAG (y
));
4154 else if (SCM_FRACTIONP (y
))
4155 /* a/b + c/d = (ad + bc) / bd */
4156 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4157 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4158 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4160 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4163 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4167 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4169 "Return @math{@var{x}+1}.")
4170 #define FUNC_NAME s_scm_oneplus
4172 return scm_sum (x
, SCM_I_MAKINUM (1));
4177 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4178 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4179 * the sum of all but the first argument are subtracted from the first
4181 #define FUNC_NAME s_difference
4183 scm_difference (SCM x
, SCM y
)
4185 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4188 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4190 if (SCM_I_INUMP (x
))
4192 long xx
= -SCM_I_INUM (x
);
4193 if (SCM_FIXABLE (xx
))
4194 return SCM_I_MAKINUM (xx
);
4196 return scm_i_long2big (xx
);
4198 else if (SCM_BIGP (x
))
4199 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4200 bignum, but negating that gives a fixnum. */
4201 return scm_i_normbig (scm_i_clonebig (x
, 0));
4202 else if (SCM_REALP (x
))
4203 return scm_from_double (-SCM_REAL_VALUE (x
));
4204 else if (SCM_COMPLEXP (x
))
4205 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4206 -SCM_COMPLEX_IMAG (x
));
4207 else if (SCM_FRACTIONP (x
))
4208 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4209 SCM_FRACTION_DENOMINATOR (x
));
4211 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4214 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4216 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4218 long int xx
= SCM_I_INUM (x
);
4219 long int yy
= SCM_I_INUM (y
);
4220 long int z
= xx
- yy
;
4221 if (SCM_FIXABLE (z
))
4222 return SCM_I_MAKINUM (z
);
4224 return scm_i_long2big (z
);
4226 else if (SCM_BIGP (y
))
4228 /* inum-x - big-y */
4229 long xx
= SCM_I_INUM (x
);
4232 return scm_i_clonebig (y
, 0);
4235 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4236 SCM result
= scm_i_mkbig ();
4239 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4242 /* x - y == -(y + -x) */
4243 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4244 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4246 scm_remember_upto_here_1 (y
);
4248 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4249 /* we know the result will have to be a bignum */
4252 return scm_i_normbig (result
);
4255 else if (SCM_REALP (y
))
4257 long int xx
= SCM_I_INUM (x
);
4258 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4260 else if (SCM_COMPLEXP (y
))
4262 long int xx
= SCM_I_INUM (x
);
4263 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4264 - SCM_COMPLEX_IMAG (y
));
4266 else if (SCM_FRACTIONP (y
))
4267 /* a - b/c = (ac - b) / c */
4268 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4269 SCM_FRACTION_NUMERATOR (y
)),
4270 SCM_FRACTION_DENOMINATOR (y
));
4272 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4274 else if (SCM_BIGP (x
))
4276 if (SCM_I_INUMP (y
))
4278 /* big-x - inum-y */
4279 long yy
= SCM_I_INUM (y
);
4280 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4282 scm_remember_upto_here_1 (x
);
4284 return (SCM_FIXABLE (-yy
) ?
4285 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4288 SCM result
= scm_i_mkbig ();
4291 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4293 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4294 scm_remember_upto_here_1 (x
);
4296 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4297 /* we know the result will have to be a bignum */
4300 return scm_i_normbig (result
);
4303 else if (SCM_BIGP (y
))
4305 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4306 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4307 SCM result
= scm_i_mkbig ();
4308 mpz_sub (SCM_I_BIG_MPZ (result
),
4311 scm_remember_upto_here_2 (x
, y
);
4312 /* we know the result will have to be a bignum */
4313 if ((sgn_x
== 1) && (sgn_y
== -1))
4315 if ((sgn_x
== -1) && (sgn_y
== 1))
4317 return scm_i_normbig (result
);
4319 else if (SCM_REALP (y
))
4321 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4322 scm_remember_upto_here_1 (x
);
4323 return scm_from_double (result
);
4325 else if (SCM_COMPLEXP (y
))
4327 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4328 - SCM_COMPLEX_REAL (y
));
4329 scm_remember_upto_here_1 (x
);
4330 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4332 else if (SCM_FRACTIONP (y
))
4333 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4334 SCM_FRACTION_NUMERATOR (y
)),
4335 SCM_FRACTION_DENOMINATOR (y
));
4336 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4338 else if (SCM_REALP (x
))
4340 if (SCM_I_INUMP (y
))
4341 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4342 else if (SCM_BIGP (y
))
4344 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4345 scm_remember_upto_here_1 (x
);
4346 return scm_from_double (result
);
4348 else if (SCM_REALP (y
))
4349 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4350 else if (SCM_COMPLEXP (y
))
4351 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4352 -SCM_COMPLEX_IMAG (y
));
4353 else if (SCM_FRACTIONP (y
))
4354 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4356 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4358 else if (SCM_COMPLEXP (x
))
4360 if (SCM_I_INUMP (y
))
4361 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4362 SCM_COMPLEX_IMAG (x
));
4363 else if (SCM_BIGP (y
))
4365 double real_part
= (SCM_COMPLEX_REAL (x
)
4366 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4367 scm_remember_upto_here_1 (x
);
4368 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4370 else if (SCM_REALP (y
))
4371 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4372 SCM_COMPLEX_IMAG (x
));
4373 else if (SCM_COMPLEXP (y
))
4374 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4375 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4376 else if (SCM_FRACTIONP (y
))
4377 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4378 SCM_COMPLEX_IMAG (x
));
4380 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4382 else if (SCM_FRACTIONP (x
))
4384 if (SCM_I_INUMP (y
))
4385 /* a/b - c = (a - cb) / b */
4386 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4387 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4388 SCM_FRACTION_DENOMINATOR (x
));
4389 else if (SCM_BIGP (y
))
4390 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4391 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4392 SCM_FRACTION_DENOMINATOR (x
));
4393 else if (SCM_REALP (y
))
4394 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4395 else if (SCM_COMPLEXP (y
))
4396 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4397 -SCM_COMPLEX_IMAG (y
));
4398 else if (SCM_FRACTIONP (y
))
4399 /* a/b - c/d = (ad - bc) / bd */
4400 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4401 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4402 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4404 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4407 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4412 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4414 "Return @math{@var{x}-1}.")
4415 #define FUNC_NAME s_scm_oneminus
4417 return scm_difference (x
, SCM_I_MAKINUM (1));
4422 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4423 /* "Return the product of all arguments. If called without arguments,\n"
4427 scm_product (SCM x
, SCM y
)
4429 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4432 return SCM_I_MAKINUM (1L);
4433 else if (SCM_NUMBERP (x
))
4436 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4439 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4444 xx
= SCM_I_INUM (x
);
4448 case 0: return x
; break;
4449 case 1: return y
; break;
4452 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4454 long yy
= SCM_I_INUM (y
);
4456 SCM k
= SCM_I_MAKINUM (kk
);
4457 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4461 SCM result
= scm_i_long2big (xx
);
4462 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4463 return scm_i_normbig (result
);
4466 else if (SCM_BIGP (y
))
4468 SCM result
= scm_i_mkbig ();
4469 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4470 scm_remember_upto_here_1 (y
);
4473 else if (SCM_REALP (y
))
4474 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4475 else if (SCM_COMPLEXP (y
))
4476 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4477 xx
* SCM_COMPLEX_IMAG (y
));
4478 else if (SCM_FRACTIONP (y
))
4479 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4480 SCM_FRACTION_DENOMINATOR (y
));
4482 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4484 else if (SCM_BIGP (x
))
4486 if (SCM_I_INUMP (y
))
4491 else if (SCM_BIGP (y
))
4493 SCM result
= scm_i_mkbig ();
4494 mpz_mul (SCM_I_BIG_MPZ (result
),
4497 scm_remember_upto_here_2 (x
, y
);
4500 else if (SCM_REALP (y
))
4502 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4503 scm_remember_upto_here_1 (x
);
4504 return scm_from_double (result
);
4506 else if (SCM_COMPLEXP (y
))
4508 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4509 scm_remember_upto_here_1 (x
);
4510 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4511 z
* SCM_COMPLEX_IMAG (y
));
4513 else if (SCM_FRACTIONP (y
))
4514 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4515 SCM_FRACTION_DENOMINATOR (y
));
4517 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4519 else if (SCM_REALP (x
))
4521 if (SCM_I_INUMP (y
))
4523 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4524 if (scm_is_eq (y
, SCM_INUM0
))
4526 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4528 else if (SCM_BIGP (y
))
4530 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4531 scm_remember_upto_here_1 (y
);
4532 return scm_from_double (result
);
4534 else if (SCM_REALP (y
))
4535 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4536 else if (SCM_COMPLEXP (y
))
4537 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4538 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4539 else if (SCM_FRACTIONP (y
))
4540 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4542 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4544 else if (SCM_COMPLEXP (x
))
4546 if (SCM_I_INUMP (y
))
4548 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4549 if (scm_is_eq (y
, SCM_INUM0
))
4551 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4552 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4554 else if (SCM_BIGP (y
))
4556 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4557 scm_remember_upto_here_1 (y
);
4558 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4559 z
* SCM_COMPLEX_IMAG (x
));
4561 else if (SCM_REALP (y
))
4562 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4563 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4564 else if (SCM_COMPLEXP (y
))
4566 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4567 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4568 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4569 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4571 else if (SCM_FRACTIONP (y
))
4573 double yy
= scm_i_fraction2double (y
);
4574 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4575 yy
* SCM_COMPLEX_IMAG (x
));
4578 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4580 else if (SCM_FRACTIONP (x
))
4582 if (SCM_I_INUMP (y
))
4583 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4584 SCM_FRACTION_DENOMINATOR (x
));
4585 else if (SCM_BIGP (y
))
4586 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4587 SCM_FRACTION_DENOMINATOR (x
));
4588 else if (SCM_REALP (y
))
4589 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4590 else if (SCM_COMPLEXP (y
))
4592 double xx
= scm_i_fraction2double (x
);
4593 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4594 xx
* SCM_COMPLEX_IMAG (y
));
4596 else if (SCM_FRACTIONP (y
))
4597 /* a/b * c/d = ac / bd */
4598 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4599 SCM_FRACTION_NUMERATOR (y
)),
4600 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4601 SCM_FRACTION_DENOMINATOR (y
)));
4603 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4606 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4609 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4610 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4611 #define ALLOW_DIVIDE_BY_ZERO
4612 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4615 /* The code below for complex division is adapted from the GNU
4616 libstdc++, which adapted it from f2c's libF77, and is subject to
4619 /****************************************************************
4620 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4622 Permission to use, copy, modify, and distribute this software
4623 and its documentation for any purpose and without fee is hereby
4624 granted, provided that the above copyright notice appear in all
4625 copies and that both that the copyright notice and this
4626 permission notice and warranty disclaimer appear in supporting
4627 documentation, and that the names of AT&T Bell Laboratories or
4628 Bellcore or any of their entities not be used in advertising or
4629 publicity pertaining to distribution of the software without
4630 specific, written prior permission.
4632 AT&T and Bellcore disclaim all warranties with regard to this
4633 software, including all implied warranties of merchantability
4634 and fitness. In no event shall AT&T or Bellcore be liable for
4635 any special, indirect or consequential damages or any damages
4636 whatsoever resulting from loss of use, data or profits, whether
4637 in an action of contract, negligence or other tortious action,
4638 arising out of or in connection with the use or performance of
4640 ****************************************************************/
4642 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4643 /* Divide the first argument by the product of the remaining
4644 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4646 #define FUNC_NAME s_divide
4648 scm_i_divide (SCM x
, SCM y
, int inexact
)
4652 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4655 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4656 else if (SCM_I_INUMP (x
))
4658 long xx
= SCM_I_INUM (x
);
4659 if (xx
== 1 || xx
== -1)
4661 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4663 scm_num_overflow (s_divide
);
4668 return scm_from_double (1.0 / (double) xx
);
4669 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4672 else if (SCM_BIGP (x
))
4675 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4676 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4678 else if (SCM_REALP (x
))
4680 double xx
= SCM_REAL_VALUE (x
);
4681 #ifndef ALLOW_DIVIDE_BY_ZERO
4683 scm_num_overflow (s_divide
);
4686 return scm_from_double (1.0 / xx
);
4688 else if (SCM_COMPLEXP (x
))
4690 double r
= SCM_COMPLEX_REAL (x
);
4691 double i
= SCM_COMPLEX_IMAG (x
);
4692 if (fabs(r
) <= fabs(i
))
4695 double d
= i
* (1.0 + t
* t
);
4696 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4701 double d
= r
* (1.0 + t
* t
);
4702 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4705 else if (SCM_FRACTIONP (x
))
4706 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4707 SCM_FRACTION_NUMERATOR (x
));
4709 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4712 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4714 long xx
= SCM_I_INUM (x
);
4715 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4717 long yy
= SCM_I_INUM (y
);
4720 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4721 scm_num_overflow (s_divide
);
4723 return scm_from_double ((double) xx
/ (double) yy
);
4726 else if (xx
% yy
!= 0)
4729 return scm_from_double ((double) xx
/ (double) yy
);
4730 else return scm_i_make_ratio (x
, y
);
4735 if (SCM_FIXABLE (z
))
4736 return SCM_I_MAKINUM (z
);
4738 return scm_i_long2big (z
);
4741 else if (SCM_BIGP (y
))
4744 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4745 else return scm_i_make_ratio (x
, y
);
4747 else if (SCM_REALP (y
))
4749 double yy
= SCM_REAL_VALUE (y
);
4750 #ifndef ALLOW_DIVIDE_BY_ZERO
4752 scm_num_overflow (s_divide
);
4755 return scm_from_double ((double) xx
/ yy
);
4757 else if (SCM_COMPLEXP (y
))
4760 complex_div
: /* y _must_ be a complex number */
4762 double r
= SCM_COMPLEX_REAL (y
);
4763 double i
= SCM_COMPLEX_IMAG (y
);
4764 if (fabs(r
) <= fabs(i
))
4767 double d
= i
* (1.0 + t
* t
);
4768 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4773 double d
= r
* (1.0 + t
* t
);
4774 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4778 else if (SCM_FRACTIONP (y
))
4779 /* a / b/c = ac / b */
4780 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4781 SCM_FRACTION_NUMERATOR (y
));
4783 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4785 else if (SCM_BIGP (x
))
4787 if (SCM_I_INUMP (y
))
4789 long int yy
= SCM_I_INUM (y
);
4792 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4793 scm_num_overflow (s_divide
);
4795 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4796 scm_remember_upto_here_1 (x
);
4797 return (sgn
== 0) ? scm_nan () : scm_inf ();
4804 /* FIXME: HMM, what are the relative performance issues here?
4805 We need to test. Is it faster on average to test
4806 divisible_p, then perform whichever operation, or is it
4807 faster to perform the integer div opportunistically and
4808 switch to real if there's a remainder? For now we take the
4809 middle ground: test, then if divisible, use the faster div
4812 long abs_yy
= yy
< 0 ? -yy
: yy
;
4813 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4817 SCM result
= scm_i_mkbig ();
4818 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4819 scm_remember_upto_here_1 (x
);
4821 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4822 return scm_i_normbig (result
);
4827 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4828 else return scm_i_make_ratio (x
, y
);
4832 else if (SCM_BIGP (y
))
4834 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4837 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4838 scm_num_overflow (s_divide
);
4840 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4841 scm_remember_upto_here_1 (x
);
4842 return (sgn
== 0) ? scm_nan () : scm_inf ();
4850 /* It's easily possible for the ratio x/y to fit a double
4851 but one or both x and y be too big to fit a double,
4852 hence the use of mpq_get_d rather than converting and
4855 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4856 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4857 return scm_from_double (mpq_get_d (q
));
4861 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4865 SCM result
= scm_i_mkbig ();
4866 mpz_divexact (SCM_I_BIG_MPZ (result
),
4869 scm_remember_upto_here_2 (x
, y
);
4870 return scm_i_normbig (result
);
4873 return scm_i_make_ratio (x
, y
);
4877 else if (SCM_REALP (y
))
4879 double yy
= SCM_REAL_VALUE (y
);
4880 #ifndef ALLOW_DIVIDE_BY_ZERO
4882 scm_num_overflow (s_divide
);
4885 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4887 else if (SCM_COMPLEXP (y
))
4889 a
= scm_i_big2dbl (x
);
4892 else if (SCM_FRACTIONP (y
))
4893 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4894 SCM_FRACTION_NUMERATOR (y
));
4896 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4898 else if (SCM_REALP (x
))
4900 double rx
= SCM_REAL_VALUE (x
);
4901 if (SCM_I_INUMP (y
))
4903 long int yy
= SCM_I_INUM (y
);
4904 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4906 scm_num_overflow (s_divide
);
4909 return scm_from_double (rx
/ (double) yy
);
4911 else if (SCM_BIGP (y
))
4913 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4914 scm_remember_upto_here_1 (y
);
4915 return scm_from_double (rx
/ dby
);
4917 else if (SCM_REALP (y
))
4919 double yy
= SCM_REAL_VALUE (y
);
4920 #ifndef ALLOW_DIVIDE_BY_ZERO
4922 scm_num_overflow (s_divide
);
4925 return scm_from_double (rx
/ yy
);
4927 else if (SCM_COMPLEXP (y
))
4932 else if (SCM_FRACTIONP (y
))
4933 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4935 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4937 else if (SCM_COMPLEXP (x
))
4939 double rx
= SCM_COMPLEX_REAL (x
);
4940 double ix
= SCM_COMPLEX_IMAG (x
);
4941 if (SCM_I_INUMP (y
))
4943 long int yy
= SCM_I_INUM (y
);
4944 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4946 scm_num_overflow (s_divide
);
4951 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4954 else if (SCM_BIGP (y
))
4956 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4957 scm_remember_upto_here_1 (y
);
4958 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4960 else if (SCM_REALP (y
))
4962 double yy
= SCM_REAL_VALUE (y
);
4963 #ifndef ALLOW_DIVIDE_BY_ZERO
4965 scm_num_overflow (s_divide
);
4968 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4970 else if (SCM_COMPLEXP (y
))
4972 double ry
= SCM_COMPLEX_REAL (y
);
4973 double iy
= SCM_COMPLEX_IMAG (y
);
4974 if (fabs(ry
) <= fabs(iy
))
4977 double d
= iy
* (1.0 + t
* t
);
4978 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4983 double d
= ry
* (1.0 + t
* t
);
4984 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4987 else if (SCM_FRACTIONP (y
))
4989 double yy
= scm_i_fraction2double (y
);
4990 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4993 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4995 else if (SCM_FRACTIONP (x
))
4997 if (SCM_I_INUMP (y
))
4999 long int yy
= SCM_I_INUM (y
);
5000 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5002 scm_num_overflow (s_divide
);
5005 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5006 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5008 else if (SCM_BIGP (y
))
5010 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5011 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5013 else if (SCM_REALP (y
))
5015 double yy
= SCM_REAL_VALUE (y
);
5016 #ifndef ALLOW_DIVIDE_BY_ZERO
5018 scm_num_overflow (s_divide
);
5021 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5023 else if (SCM_COMPLEXP (y
))
5025 a
= scm_i_fraction2double (x
);
5028 else if (SCM_FRACTIONP (y
))
5029 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5030 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5032 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5035 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5039 scm_divide (SCM x
, SCM y
)
5041 return scm_i_divide (x
, y
, 0);
5044 static SCM
scm_divide2real (SCM x
, SCM y
)
5046 return scm_i_divide (x
, y
, 1);
5052 scm_asinh (double x
)
5057 #define asinh scm_asinh
5058 return log (x
+ sqrt (x
* x
+ 1));
5061 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5062 /* "Return the inverse hyperbolic sine of @var{x}."
5067 scm_acosh (double x
)
5072 #define acosh scm_acosh
5073 return log (x
+ sqrt (x
* x
- 1));
5076 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5077 /* "Return the inverse hyperbolic cosine of @var{x}."
5082 scm_atanh (double x
)
5087 #define atanh scm_atanh
5088 return 0.5 * log ((1 + x
) / (1 - x
));
5091 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5092 /* "Return the inverse hyperbolic tangent of @var{x}."
5097 scm_c_truncate (double x
)
5108 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5109 half-way case (ie. when x is an integer plus 0.5) going upwards.
5110 Then half-way cases are identified and adjusted down if the
5111 round-upwards didn't give the desired even integer.
5113 "plus_half == result" identifies a half-way case. If plus_half, which is
5114 x + 0.5, is an integer then x must be an integer plus 0.5.
5116 An odd "result" value is identified with result/2 != floor(result/2).
5117 This is done with plus_half, since that value is ready for use sooner in
5118 a pipelined cpu, and we're already requiring plus_half == result.
5120 Note however that we need to be careful when x is big and already an
5121 integer. In that case "x+0.5" may round to an adjacent integer, causing
5122 us to return such a value, incorrectly. For instance if the hardware is
5123 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5124 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5125 returned. Or if the hardware is in round-upwards mode, then other bigger
5126 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5127 representable value, 2^128+2^76 (or whatever), again incorrect.
5129 These bad roundings of x+0.5 are avoided by testing at the start whether
5130 x is already an integer. If it is then clearly that's the desired result
5131 already. And if it's not then the exponent must be small enough to allow
5132 an 0.5 to be represented, and hence added without a bad rounding. */
5135 scm_c_round (double x
)
5137 double plus_half
, result
;
5142 plus_half
= x
+ 0.5;
5143 result
= floor (plus_half
);
5144 /* Adjust so that the rounding is towards even. */
5145 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5150 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5152 "Round the number @var{x} towards zero.")
5153 #define FUNC_NAME s_scm_truncate_number
5155 if (scm_is_false (scm_negative_p (x
)))
5156 return scm_floor (x
);
5158 return scm_ceiling (x
);
5162 static SCM exactly_one_half
;
5164 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5166 "Round the number @var{x} towards the nearest integer. "
5167 "When it is exactly halfway between two integers, "
5168 "round towards the even one.")
5169 #define FUNC_NAME s_scm_round_number
5171 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5173 else if (SCM_REALP (x
))
5174 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5177 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5178 single quotient+remainder division then examining to see which way
5179 the rounding should go. */
5180 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5181 SCM result
= scm_floor (plus_half
);
5182 /* Adjust so that the rounding is towards even. */
5183 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5184 && scm_is_true (scm_odd_p (result
)))
5185 return scm_difference (result
, SCM_I_MAKINUM (1));
5192 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5194 "Round the number @var{x} towards minus infinity.")
5195 #define FUNC_NAME s_scm_floor
5197 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5199 else if (SCM_REALP (x
))
5200 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5201 else if (SCM_FRACTIONP (x
))
5203 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5204 SCM_FRACTION_DENOMINATOR (x
));
5205 if (scm_is_false (scm_negative_p (x
)))
5207 /* For positive x, rounding towards zero is correct. */
5212 /* For negative x, we need to return q-1 unless x is an
5213 integer. But fractions are never integer, per our
5215 return scm_difference (q
, SCM_I_MAKINUM (1));
5219 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5223 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5225 "Round the number @var{x} towards infinity.")
5226 #define FUNC_NAME s_scm_ceiling
5228 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5230 else if (SCM_REALP (x
))
5231 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5232 else if (SCM_FRACTIONP (x
))
5234 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5235 SCM_FRACTION_DENOMINATOR (x
));
5236 if (scm_is_false (scm_positive_p (x
)))
5238 /* For negative x, rounding towards zero is correct. */
5243 /* For positive x, we need to return q+1 unless x is an
5244 integer. But fractions are never integer, per our
5246 return scm_sum (q
, SCM_I_MAKINUM (1));
5250 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5254 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5255 /* "Return the square root of the real number @var{x}."
5257 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5258 /* "Return the absolute value of the real number @var{x}."
5260 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5261 /* "Return the @var{x}th power of e."
5263 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5264 /* "Return the natural logarithm of the real number @var{x}."
5266 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5267 /* "Return the sine of the real number @var{x}."
5269 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5270 /* "Return the cosine of the real number @var{x}."
5272 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5273 /* "Return the tangent of the real number @var{x}."
5275 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5276 /* "Return the arc sine of the real number @var{x}."
5278 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5279 /* "Return the arc cosine of the real number @var{x}."
5281 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5282 /* "Return the arc tangent of the real number @var{x}."
5284 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5285 /* "Return the hyperbolic sine of the real number @var{x}."
5287 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5288 /* "Return the hyperbolic cosine of the real number @var{x}."
5290 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5291 /* "Return the hyperbolic tangent of the real number @var{x}."
5299 static void scm_two_doubles (SCM x
,
5301 const char *sstring
,
5305 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5307 if (SCM_I_INUMP (x
))
5308 xy
->x
= SCM_I_INUM (x
);
5309 else if (SCM_BIGP (x
))
5310 xy
->x
= scm_i_big2dbl (x
);
5311 else if (SCM_REALP (x
))
5312 xy
->x
= SCM_REAL_VALUE (x
);
5313 else if (SCM_FRACTIONP (x
))
5314 xy
->x
= scm_i_fraction2double (x
);
5316 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5318 if (SCM_I_INUMP (y
))
5319 xy
->y
= SCM_I_INUM (y
);
5320 else if (SCM_BIGP (y
))
5321 xy
->y
= scm_i_big2dbl (y
);
5322 else if (SCM_REALP (y
))
5323 xy
->y
= SCM_REAL_VALUE (y
);
5324 else if (SCM_FRACTIONP (y
))
5325 xy
->y
= scm_i_fraction2double (y
);
5327 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5331 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5333 "Return @var{x} raised to the power of @var{y}. This\n"
5334 "procedure does not accept complex arguments.")
5335 #define FUNC_NAME s_scm_sys_expt
5338 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5339 return scm_from_double (pow (xy
.x
, xy
.y
));
5344 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5346 "Return the arc tangent of the two arguments @var{x} and\n"
5347 "@var{y}. This is similar to calculating the arc tangent of\n"
5348 "@var{x} / @var{y}, except that the signs of both arguments\n"
5349 "are used to determine the quadrant of the result. This\n"
5350 "procedure does not accept complex arguments.")
5351 #define FUNC_NAME s_scm_sys_atan2
5354 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5355 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5360 scm_c_make_rectangular (double re
, double im
)
5363 return scm_from_double (re
);
5367 SCM_NEWSMOB (z
, scm_tc16_complex
,
5368 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5370 SCM_COMPLEX_REAL (z
) = re
;
5371 SCM_COMPLEX_IMAG (z
) = im
;
5376 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5377 (SCM real_part
, SCM imaginary_part
),
5378 "Return a complex number constructed of the given @var{real-part} "
5379 "and @var{imaginary-part} parts.")
5380 #define FUNC_NAME s_scm_make_rectangular
5383 scm_two_doubles (real_part
, imaginary_part
, FUNC_NAME
, &xy
);
5384 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5389 scm_c_make_polar (double mag
, double ang
)
5393 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5394 use it on Glibc-based systems that have it (it's a GNU extension). See
5395 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5397 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5398 sincos (ang
, &s
, &c
);
5403 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5406 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5408 "Return the complex number @var{x} * e^(i * @var{y}).")
5409 #define FUNC_NAME s_scm_make_polar
5412 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5413 return scm_c_make_polar (xy
.x
, xy
.y
);
5418 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5419 /* "Return the real part of the number @var{z}."
5422 scm_real_part (SCM z
)
5424 if (SCM_I_INUMP (z
))
5426 else if (SCM_BIGP (z
))
5428 else if (SCM_REALP (z
))
5430 else if (SCM_COMPLEXP (z
))
5431 return scm_from_double (SCM_COMPLEX_REAL (z
));
5432 else if (SCM_FRACTIONP (z
))
5435 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5439 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5440 /* "Return the imaginary part of the number @var{z}."
5443 scm_imag_part (SCM z
)
5445 if (SCM_I_INUMP (z
))
5447 else if (SCM_BIGP (z
))
5449 else if (SCM_REALP (z
))
5451 else if (SCM_COMPLEXP (z
))
5452 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5453 else if (SCM_FRACTIONP (z
))
5456 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5459 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5460 /* "Return the numerator of the number @var{z}."
5463 scm_numerator (SCM z
)
5465 if (SCM_I_INUMP (z
))
5467 else if (SCM_BIGP (z
))
5469 else if (SCM_FRACTIONP (z
))
5470 return SCM_FRACTION_NUMERATOR (z
);
5471 else if (SCM_REALP (z
))
5472 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5474 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5478 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5479 /* "Return the denominator of the number @var{z}."
5482 scm_denominator (SCM z
)
5484 if (SCM_I_INUMP (z
))
5485 return SCM_I_MAKINUM (1);
5486 else if (SCM_BIGP (z
))
5487 return SCM_I_MAKINUM (1);
5488 else if (SCM_FRACTIONP (z
))
5489 return SCM_FRACTION_DENOMINATOR (z
);
5490 else if (SCM_REALP (z
))
5491 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5493 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5496 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5497 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5498 * "@code{abs} for real arguments, but also allows complex numbers."
5501 scm_magnitude (SCM z
)
5503 if (SCM_I_INUMP (z
))
5505 long int zz
= SCM_I_INUM (z
);
5508 else if (SCM_POSFIXABLE (-zz
))
5509 return SCM_I_MAKINUM (-zz
);
5511 return scm_i_long2big (-zz
);
5513 else if (SCM_BIGP (z
))
5515 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5516 scm_remember_upto_here_1 (z
);
5518 return scm_i_clonebig (z
, 0);
5522 else if (SCM_REALP (z
))
5523 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5524 else if (SCM_COMPLEXP (z
))
5525 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5526 else if (SCM_FRACTIONP (z
))
5528 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5530 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5531 SCM_FRACTION_DENOMINATOR (z
));
5534 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5538 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5539 /* "Return the angle of the complex number @var{z}."
5544 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5545 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5546 But if atan2 follows the floating point rounding mode, then the value
5547 is not a constant. Maybe it'd be close enough though. */
5548 if (SCM_I_INUMP (z
))
5550 if (SCM_I_INUM (z
) >= 0)
5553 return scm_from_double (atan2 (0.0, -1.0));
5555 else if (SCM_BIGP (z
))
5557 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5558 scm_remember_upto_here_1 (z
);
5560 return scm_from_double (atan2 (0.0, -1.0));
5564 else if (SCM_REALP (z
))
5566 if (SCM_REAL_VALUE (z
) >= 0)
5569 return scm_from_double (atan2 (0.0, -1.0));
5571 else if (SCM_COMPLEXP (z
))
5572 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5573 else if (SCM_FRACTIONP (z
))
5575 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5577 else return scm_from_double (atan2 (0.0, -1.0));
5580 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5584 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5585 /* Convert the number @var{x} to its inexact representation.\n"
5588 scm_exact_to_inexact (SCM z
)
5590 if (SCM_I_INUMP (z
))
5591 return scm_from_double ((double) SCM_I_INUM (z
));
5592 else if (SCM_BIGP (z
))
5593 return scm_from_double (scm_i_big2dbl (z
));
5594 else if (SCM_FRACTIONP (z
))
5595 return scm_from_double (scm_i_fraction2double (z
));
5596 else if (SCM_INEXACTP (z
))
5599 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5603 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5605 "Return an exact number that is numerically closest to @var{z}.")
5606 #define FUNC_NAME s_scm_inexact_to_exact
5608 if (SCM_I_INUMP (z
))
5610 else if (SCM_BIGP (z
))
5612 else if (SCM_REALP (z
))
5614 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5615 SCM_OUT_OF_RANGE (1, z
);
5622 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5623 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5624 scm_i_mpz2num (mpq_denref (frac
)));
5626 /* When scm_i_make_ratio throws, we leak the memory allocated
5633 else if (SCM_FRACTIONP (z
))
5636 SCM_WRONG_TYPE_ARG (1, z
);
5640 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5642 "Returns the @emph{simplest} rational number differing\n"
5643 "from @var{x} by no more than @var{eps}.\n"
5645 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5646 "exact result when both its arguments are exact. Thus, you might need\n"
5647 "to use @code{inexact->exact} on the arguments.\n"
5650 "(rationalize (inexact->exact 1.2) 1/100)\n"
5653 #define FUNC_NAME s_scm_rationalize
5655 if (SCM_I_INUMP (x
))
5657 else if (SCM_BIGP (x
))
5659 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5661 /* Use continued fractions to find closest ratio. All
5662 arithmetic is done with exact numbers.
5665 SCM ex
= scm_inexact_to_exact (x
);
5666 SCM int_part
= scm_floor (ex
);
5667 SCM tt
= SCM_I_MAKINUM (1);
5668 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5669 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5673 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5676 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5677 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5679 /* We stop after a million iterations just to be absolutely sure
5680 that we don't go into an infinite loop. The process normally
5681 converges after less than a dozen iterations.
5684 eps
= scm_abs (eps
);
5685 while (++i
< 1000000)
5687 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5688 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5689 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5691 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5692 eps
))) /* abs(x-a/b) <= eps */
5694 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5695 if (scm_is_false (scm_exact_p (x
))
5696 || scm_is_false (scm_exact_p (eps
)))
5697 return scm_exact_to_inexact (res
);
5701 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5703 tt
= scm_floor (rx
); /* tt = floor (rx) */
5709 scm_num_overflow (s_scm_rationalize
);
5712 SCM_WRONG_TYPE_ARG (1, x
);
5716 /* conversion functions */
5719 scm_is_integer (SCM val
)
5721 return scm_is_true (scm_integer_p (val
));
5725 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5727 if (SCM_I_INUMP (val
))
5729 scm_t_signed_bits n
= SCM_I_INUM (val
);
5730 return n
>= min
&& n
<= max
;
5732 else if (SCM_BIGP (val
))
5734 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5736 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5738 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5740 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5741 return n
>= min
&& n
<= max
;
5751 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5752 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5755 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5756 SCM_I_BIG_MPZ (val
));
5758 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5770 return n
>= min
&& n
<= max
;
5778 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5780 if (SCM_I_INUMP (val
))
5782 scm_t_signed_bits n
= SCM_I_INUM (val
);
5783 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5785 else if (SCM_BIGP (val
))
5787 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5789 else if (max
<= ULONG_MAX
)
5791 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5793 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5794 return n
>= min
&& n
<= max
;
5804 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5807 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5808 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5811 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5812 SCM_I_BIG_MPZ (val
));
5814 return n
>= min
&& n
<= max
;
5822 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5824 scm_error (scm_out_of_range_key
,
5826 "Value out of range ~S to ~S: ~S",
5827 scm_list_3 (min
, max
, bad_val
),
5828 scm_list_1 (bad_val
));
5831 #define TYPE scm_t_intmax
5832 #define TYPE_MIN min
5833 #define TYPE_MAX max
5834 #define SIZEOF_TYPE 0
5835 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5836 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5837 #include "libguile/conv-integer.i.c"
5839 #define TYPE scm_t_uintmax
5840 #define TYPE_MIN min
5841 #define TYPE_MAX max
5842 #define SIZEOF_TYPE 0
5843 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5844 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5845 #include "libguile/conv-uinteger.i.c"
5847 #define TYPE scm_t_int8
5848 #define TYPE_MIN SCM_T_INT8_MIN
5849 #define TYPE_MAX SCM_T_INT8_MAX
5850 #define SIZEOF_TYPE 1
5851 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5852 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5853 #include "libguile/conv-integer.i.c"
5855 #define TYPE scm_t_uint8
5857 #define TYPE_MAX SCM_T_UINT8_MAX
5858 #define SIZEOF_TYPE 1
5859 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5860 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5861 #include "libguile/conv-uinteger.i.c"
5863 #define TYPE scm_t_int16
5864 #define TYPE_MIN SCM_T_INT16_MIN
5865 #define TYPE_MAX SCM_T_INT16_MAX
5866 #define SIZEOF_TYPE 2
5867 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5868 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5869 #include "libguile/conv-integer.i.c"
5871 #define TYPE scm_t_uint16
5873 #define TYPE_MAX SCM_T_UINT16_MAX
5874 #define SIZEOF_TYPE 2
5875 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5876 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5877 #include "libguile/conv-uinteger.i.c"
5879 #define TYPE scm_t_int32
5880 #define TYPE_MIN SCM_T_INT32_MIN
5881 #define TYPE_MAX SCM_T_INT32_MAX
5882 #define SIZEOF_TYPE 4
5883 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5884 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5885 #include "libguile/conv-integer.i.c"
5887 #define TYPE scm_t_uint32
5889 #define TYPE_MAX SCM_T_UINT32_MAX
5890 #define SIZEOF_TYPE 4
5891 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5892 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5893 #include "libguile/conv-uinteger.i.c"
5895 #define TYPE scm_t_wchar
5896 #define TYPE_MIN (scm_t_int32)-1
5897 #define TYPE_MAX (scm_t_int32)0x10ffff
5898 #define SIZEOF_TYPE 4
5899 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
5900 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
5901 #include "libguile/conv-integer.i.c"
5903 #if SCM_HAVE_T_INT64
5905 #define TYPE scm_t_int64
5906 #define TYPE_MIN SCM_T_INT64_MIN
5907 #define TYPE_MAX SCM_T_INT64_MAX
5908 #define SIZEOF_TYPE 8
5909 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5910 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5911 #include "libguile/conv-integer.i.c"
5913 #define TYPE scm_t_uint64
5915 #define TYPE_MAX SCM_T_UINT64_MAX
5916 #define SIZEOF_TYPE 8
5917 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5918 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5919 #include "libguile/conv-uinteger.i.c"
5924 scm_to_mpz (SCM val
, mpz_t rop
)
5926 if (SCM_I_INUMP (val
))
5927 mpz_set_si (rop
, SCM_I_INUM (val
));
5928 else if (SCM_BIGP (val
))
5929 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5931 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5935 scm_from_mpz (mpz_t val
)
5937 return scm_i_mpz2num (val
);
5941 scm_is_real (SCM val
)
5943 return scm_is_true (scm_real_p (val
));
5947 scm_is_rational (SCM val
)
5949 return scm_is_true (scm_rational_p (val
));
5953 scm_to_double (SCM val
)
5955 if (SCM_I_INUMP (val
))
5956 return SCM_I_INUM (val
);
5957 else if (SCM_BIGP (val
))
5958 return scm_i_big2dbl (val
);
5959 else if (SCM_FRACTIONP (val
))
5960 return scm_i_fraction2double (val
);
5961 else if (SCM_REALP (val
))
5962 return SCM_REAL_VALUE (val
);
5964 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5968 scm_from_double (double val
)
5970 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5971 SCM_REAL_VALUE (z
) = val
;
5975 #if SCM_ENABLE_DISCOURAGED == 1
5978 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5982 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5986 scm_out_of_range (NULL
, num
);
5989 return scm_to_double (num
);
5993 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5997 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6001 scm_out_of_range (NULL
, num
);
6004 return scm_to_double (num
);
6010 scm_is_complex (SCM val
)
6012 return scm_is_true (scm_complex_p (val
));
6016 scm_c_real_part (SCM z
)
6018 if (SCM_COMPLEXP (z
))
6019 return SCM_COMPLEX_REAL (z
);
6022 /* Use the scm_real_part to get proper error checking and
6025 return scm_to_double (scm_real_part (z
));
6030 scm_c_imag_part (SCM z
)
6032 if (SCM_COMPLEXP (z
))
6033 return SCM_COMPLEX_IMAG (z
);
6036 /* Use the scm_imag_part to get proper error checking and
6037 dispatching. The result will almost always be 0.0, but not
6040 return scm_to_double (scm_imag_part (z
));
6045 scm_c_magnitude (SCM z
)
6047 return scm_to_double (scm_magnitude (z
));
6053 return scm_to_double (scm_angle (z
));
6057 scm_is_number (SCM z
)
6059 return scm_is_true (scm_number_p (z
));
6063 /* In the following functions we dispatch to the real-arg funcs like log()
6064 when we know the arg is real, instead of just handing everything to
6065 clog() for instance. This is in case clog() doesn't optimize for a
6066 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6067 well use it to go straight to the applicable C func. */
6069 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6071 "Return the natural logarithm of @var{z}.")
6072 #define FUNC_NAME s_scm_log
6074 if (SCM_COMPLEXP (z
))
6076 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6077 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6079 double re
= SCM_COMPLEX_REAL (z
);
6080 double im
= SCM_COMPLEX_IMAG (z
);
6081 return scm_c_make_rectangular (log (hypot (re
, im
)),
6087 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6088 although the value itself overflows. */
6089 double re
= scm_to_double (z
);
6090 double l
= log (fabs (re
));
6092 return scm_from_double (l
);
6094 return scm_c_make_rectangular (l
, M_PI
);
6100 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6102 "Return the base 10 logarithm of @var{z}.")
6103 #define FUNC_NAME s_scm_log10
6105 if (SCM_COMPLEXP (z
))
6107 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6108 clog() and a multiply by M_LOG10E, rather than the fallback
6109 log10+hypot+atan2.) */
6110 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6111 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6113 double re
= SCM_COMPLEX_REAL (z
);
6114 double im
= SCM_COMPLEX_IMAG (z
);
6115 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6116 M_LOG10E
* atan2 (im
, re
));
6121 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6122 although the value itself overflows. */
6123 double re
= scm_to_double (z
);
6124 double l
= log10 (fabs (re
));
6126 return scm_from_double (l
);
6128 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6134 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6136 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6137 "base of natural logarithms (2.71828@dots{}).")
6138 #define FUNC_NAME s_scm_exp
6140 if (SCM_COMPLEXP (z
))
6142 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6143 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6145 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6146 SCM_COMPLEX_IMAG (z
));
6151 /* When z is a negative bignum the conversion to double overflows,
6152 giving -infinity, but that's ok, the exp is still 0.0. */
6153 return scm_from_double (exp (scm_to_double (z
)));
6159 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6161 "Return the square root of @var{z}. Of the two possible roots\n"
6162 "(positive and negative), the one with the a positive real part\n"
6163 "is returned, or if that's zero then a positive imaginary part.\n"
6167 "(sqrt 9.0) @result{} 3.0\n"
6168 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6169 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6170 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6172 #define FUNC_NAME s_scm_sqrt
6174 if (SCM_COMPLEXP (x
))
6176 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6177 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6179 double re
= SCM_COMPLEX_REAL (x
);
6180 double im
= SCM_COMPLEX_IMAG (x
);
6181 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6182 0.5 * atan2 (im
, re
));
6187 double xx
= scm_to_double (x
);
6189 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6191 return scm_from_double (sqrt (xx
));
6203 mpz_init_set_si (z_negative_one
, -1);
6205 /* It may be possible to tune the performance of some algorithms by using
6206 * the following constants to avoid the creation of bignums. Please, before
6207 * using these values, remember the two rules of program optimization:
6208 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6209 scm_c_define ("most-positive-fixnum",
6210 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6211 scm_c_define ("most-negative-fixnum",
6212 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6214 scm_add_feature ("complex");
6215 scm_add_feature ("inexact");
6216 scm_flo0
= scm_from_double (0.0);
6218 /* determine floating point precision */
6219 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6221 init_dblprec(&scm_dblprec
[i
-2],i
);
6222 init_fx_radix(fx_per_radix
[i
-2],i
);
6225 /* hard code precision for base 10 if the preprocessor tells us to... */
6226 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6229 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6230 SCM_I_MAKINUM (2)));
6231 #include "libguile/numbers.x"