1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc(), csqrt(), etc */
58 #include "libguile/_scm.h"
59 #include "libguile/feature.h"
60 #include "libguile/ports.h"
61 #include "libguile/root.h"
62 #include "libguile/smob.h"
63 #include "libguile/strings.h"
65 #include "libguile/validate.h"
66 #include "libguile/numbers.h"
67 #include "libguile/deprecation.h"
69 #include "libguile/eq.h"
71 #include "libguile/discouraged.h"
73 /* values per glibc, if not already defined */
75 #define M_LOG10E 0.43429448190325182765
78 #define M_PI 3.14159265358979323846
84 Wonder if this might be faster for some of our code? A switch on
85 the numtag would jump directly to the right case, and the
86 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
88 #define SCM_I_NUMTAG_NOTNUM 0
89 #define SCM_I_NUMTAG_INUM 1
90 #define SCM_I_NUMTAG_BIG scm_tc16_big
91 #define SCM_I_NUMTAG_REAL scm_tc16_real
92 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
93 #define SCM_I_NUMTAG(x) \
94 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
95 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
96 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
97 : SCM_I_NUMTAG_NOTNUM)))
99 /* the macro above will not work as is with fractions */
102 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
104 /* FLOBUFLEN is the maximum number of characters neccessary for the
105 * printed or scm_string representation of an inexact number.
107 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
110 #if ! defined (HAVE_ISNAN)
115 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
118 #if ! defined (HAVE_ISINF)
123 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
130 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
131 an explicit check. In some future gmp (don't know what version number),
132 mpz_cmp_d is supposed to do this itself. */
134 #define xmpz_cmp_d(z, d) \
135 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
137 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
140 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
141 isinf. It does have finite and isnan though, hence the use of those.
142 fpclass would be a possibility on that system too. */
146 #if defined (HAVE_ISINF)
148 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
149 return (! (finite (x
) || isnan (x
)));
158 #if defined (HAVE_ISNAN)
165 #if defined (GUILE_I)
166 /* For an SCM object Z which is a complex number (ie. satisfies
167 SCM_COMPLEXP), return its value as a C level "complex double". */
168 #define SCM_COMPLEX_VALUE(z) \
169 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
172 /* Convert a C "complex double" to an SCM value. */
173 #if HAVE_COMPLEX_DOUBLE
175 scm_from_complex_double (complex double z
)
177 return scm_c_make_rectangular (creal (z
), cimag (z
));
179 #endif /* HAVE_COMPLEX_DOUBLE */
183 static mpz_t z_negative_one
;
190 /* Return a newly created bignum. */
191 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
192 mpz_init (SCM_I_BIG_MPZ (z
));
197 scm_i_long2big (long x
)
199 /* Return a newly created bignum initialized to X. */
200 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
201 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
206 scm_i_ulong2big (unsigned long x
)
208 /* Return a newly created bignum initialized to X. */
209 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
210 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
215 scm_i_clonebig (SCM src_big
, int same_sign_p
)
217 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
218 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
219 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
221 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
226 scm_i_bigcmp (SCM x
, SCM y
)
228 /* Return neg if x < y, pos if x > y, and 0 if x == y */
229 /* presume we already know x and y are bignums */
230 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
231 scm_remember_upto_here_2 (x
, y
);
236 scm_i_dbl2big (double d
)
238 /* results are only defined if d is an integer */
239 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
240 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
244 /* Convert a integer in double representation to a SCM number. */
247 scm_i_dbl2num (double u
)
249 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
250 powers of 2, so there's no rounding when making "double" values
251 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
252 get rounded on a 64-bit machine, hence the "+1".
254 The use of floor() to force to an integer value ensures we get a
255 "numerically closest" value without depending on how a
256 double->long cast or how mpz_set_d will round. For reference,
257 double->long probably follows the hardware rounding mode,
258 mpz_set_d truncates towards zero. */
260 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
261 representable as a double? */
263 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
264 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
265 return SCM_I_MAKINUM ((long) u
);
267 return scm_i_dbl2big (u
);
270 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
271 with R5RS exact->inexact.
273 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
274 (ie. truncate towards zero), then adjust to get the closest double by
275 examining the next lower bit and adding 1 (to the absolute value) if
278 Bignums exactly half way between representable doubles are rounded to the
279 next higher absolute value (ie. away from zero). This seems like an
280 adequate interpretation of R5RS "numerically closest", and it's easier
281 and faster than a full "nearest-even" style.
283 The bit test must be done on the absolute value of the mpz_t, which means
284 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
285 negatives as twos complement.
287 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
288 following the hardware rounding mode, but applied to the absolute value
289 of the mpz_t operand. This is not what we want so we put the high
290 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
291 mpz_get_d is supposed to always truncate towards zero.
293 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
294 is a slowdown. It'd be faster to pick out the relevant high bits with
295 mpz_getlimbn if we could be bothered coding that, and if the new
296 truncating gmp doesn't come out. */
299 scm_i_big2dbl (SCM b
)
304 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
308 /* Current GMP, eg. 4.1.3, force truncation towards zero */
310 if (bits
> DBL_MANT_DIG
)
312 size_t shift
= bits
- DBL_MANT_DIG
;
313 mpz_init2 (tmp
, DBL_MANT_DIG
);
314 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
315 result
= ldexp (mpz_get_d (tmp
), shift
);
320 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
325 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
328 if (bits
> DBL_MANT_DIG
)
330 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
331 /* test bit number "pos" in absolute value */
332 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
333 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
335 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
339 scm_remember_upto_here_1 (b
);
344 scm_i_normbig (SCM b
)
346 /* convert a big back to a fixnum if it'll fit */
347 /* presume b is a bignum */
348 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
350 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
351 if (SCM_FIXABLE (val
))
352 b
= SCM_I_MAKINUM (val
);
357 static SCM_C_INLINE_KEYWORD SCM
358 scm_i_mpz2num (mpz_t b
)
360 /* convert a mpz number to a SCM number. */
361 if (mpz_fits_slong_p (b
))
363 long val
= mpz_get_si (b
);
364 if (SCM_FIXABLE (val
))
365 return SCM_I_MAKINUM (val
);
369 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
370 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
375 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
376 static SCM
scm_divide2real (SCM x
, SCM y
);
379 scm_i_make_ratio (SCM numerator
, SCM denominator
)
380 #define FUNC_NAME "make-ratio"
382 /* First make sure the arguments are proper.
384 if (SCM_I_INUMP (denominator
))
386 if (scm_is_eq (denominator
, SCM_INUM0
))
387 scm_num_overflow ("make-ratio");
388 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
393 if (!(SCM_BIGP(denominator
)))
394 SCM_WRONG_TYPE_ARG (2, denominator
);
396 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
397 SCM_WRONG_TYPE_ARG (1, numerator
);
399 /* Then flip signs so that the denominator is positive.
401 if (scm_is_true (scm_negative_p (denominator
)))
403 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
404 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
407 /* Now consider for each of the four fixnum/bignum combinations
408 whether the rational number is really an integer.
410 if (SCM_I_INUMP (numerator
))
412 long x
= SCM_I_INUM (numerator
);
413 if (scm_is_eq (numerator
, SCM_INUM0
))
415 if (SCM_I_INUMP (denominator
))
418 y
= SCM_I_INUM (denominator
);
420 return SCM_I_MAKINUM(1);
422 return SCM_I_MAKINUM (x
/ y
);
426 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
427 of that value for the denominator, as a bignum. Apart from
428 that case, abs(bignum) > abs(inum) so inum/bignum is not an
430 if (x
== SCM_MOST_NEGATIVE_FIXNUM
431 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
432 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
433 return SCM_I_MAKINUM(-1);
436 else if (SCM_BIGP (numerator
))
438 if (SCM_I_INUMP (denominator
))
440 long yy
= SCM_I_INUM (denominator
);
441 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
442 return scm_divide (numerator
, denominator
);
446 if (scm_is_eq (numerator
, denominator
))
447 return SCM_I_MAKINUM(1);
448 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
449 SCM_I_BIG_MPZ (denominator
)))
450 return scm_divide(numerator
, denominator
);
454 /* No, it's a proper fraction.
457 SCM divisor
= scm_gcd (numerator
, denominator
);
458 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
460 numerator
= scm_divide (numerator
, divisor
);
461 denominator
= scm_divide (denominator
, divisor
);
464 return scm_double_cell (scm_tc16_fraction
,
465 SCM_UNPACK (numerator
),
466 SCM_UNPACK (denominator
), 0);
472 scm_i_fraction2double (SCM z
)
474 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
475 SCM_FRACTION_DENOMINATOR (z
)));
478 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
480 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
482 #define FUNC_NAME s_scm_exact_p
488 if (SCM_FRACTIONP (x
))
492 SCM_WRONG_TYPE_ARG (1, x
);
497 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
499 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
501 #define FUNC_NAME s_scm_odd_p
505 long val
= SCM_I_INUM (n
);
506 return scm_from_bool ((val
& 1L) != 0);
508 else if (SCM_BIGP (n
))
510 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
511 scm_remember_upto_here_1 (n
);
512 return scm_from_bool (odd_p
);
514 else if (scm_is_true (scm_inf_p (n
)))
516 else if (SCM_REALP (n
))
518 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
524 SCM_WRONG_TYPE_ARG (1, n
);
527 SCM_WRONG_TYPE_ARG (1, n
);
532 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
534 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
536 #define FUNC_NAME s_scm_even_p
540 long val
= SCM_I_INUM (n
);
541 return scm_from_bool ((val
& 1L) == 0);
543 else if (SCM_BIGP (n
))
545 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
546 scm_remember_upto_here_1 (n
);
547 return scm_from_bool (even_p
);
549 else if (scm_is_true (scm_inf_p (n
)))
551 else if (SCM_REALP (n
))
553 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
559 SCM_WRONG_TYPE_ARG (1, n
);
562 SCM_WRONG_TYPE_ARG (1, n
);
566 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
568 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
569 "or @samp{-inf.0}, @code{#f} otherwise.")
570 #define FUNC_NAME s_scm_inf_p
573 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
574 else if (SCM_COMPLEXP (x
))
575 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
576 || xisinf (SCM_COMPLEX_IMAG (x
)));
582 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
584 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
586 #define FUNC_NAME s_scm_nan_p
589 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
590 else if (SCM_COMPLEXP (n
))
591 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
592 || xisnan (SCM_COMPLEX_IMAG (n
)));
598 /* Guile's idea of infinity. */
599 static double guile_Inf
;
601 /* Guile's idea of not a number. */
602 static double guile_NaN
;
605 guile_ieee_init (void)
607 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
609 /* Some version of gcc on some old version of Linux used to crash when
610 trying to make Inf and NaN. */
613 /* C99 INFINITY, when available.
614 FIXME: The standard allows for INFINITY to be something that overflows
615 at compile time. We ought to have a configure test to check for that
616 before trying to use it. (But in practice we believe this is not a
617 problem on any system guile is likely to target.) */
618 guile_Inf
= INFINITY
;
621 extern unsigned int DINFINITY
[2];
622 guile_Inf
= (*((double *) (DINFINITY
)));
629 if (guile_Inf
== tmp
)
637 #if defined (HAVE_ISNAN)
640 /* C99 NAN, when available */
645 extern unsigned int DQNAN
[2];
646 guile_NaN
= (*((double *)(DQNAN
)));
649 guile_NaN
= guile_Inf
/ guile_Inf
;
655 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
658 #define FUNC_NAME s_scm_inf
660 static int initialized
= 0;
666 return scm_from_double (guile_Inf
);
670 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
673 #define FUNC_NAME s_scm_nan
675 static int initialized
= 0;
681 return scm_from_double (guile_NaN
);
686 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
688 "Return the absolute value of @var{x}.")
693 long int xx
= SCM_I_INUM (x
);
696 else if (SCM_POSFIXABLE (-xx
))
697 return SCM_I_MAKINUM (-xx
);
699 return scm_i_long2big (-xx
);
701 else if (SCM_BIGP (x
))
703 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
705 return scm_i_clonebig (x
, 0);
709 else if (SCM_REALP (x
))
711 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
712 double xx
= SCM_REAL_VALUE (x
);
714 return scm_from_double (-xx
);
718 else if (SCM_FRACTIONP (x
))
720 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
722 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
723 SCM_FRACTION_DENOMINATOR (x
));
726 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
731 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
732 /* "Return the quotient of the numbers @var{x} and @var{y}."
735 scm_quotient (SCM x
, SCM y
)
739 long xx
= SCM_I_INUM (x
);
742 long yy
= SCM_I_INUM (y
);
744 scm_num_overflow (s_quotient
);
749 return SCM_I_MAKINUM (z
);
751 return scm_i_long2big (z
);
754 else if (SCM_BIGP (y
))
756 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
757 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
758 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
760 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
761 scm_remember_upto_here_1 (y
);
762 return SCM_I_MAKINUM (-1);
765 return SCM_I_MAKINUM (0);
768 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
770 else if (SCM_BIGP (x
))
774 long yy
= SCM_I_INUM (y
);
776 scm_num_overflow (s_quotient
);
781 SCM result
= scm_i_mkbig ();
784 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
787 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
790 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
791 scm_remember_upto_here_1 (x
);
792 return scm_i_normbig (result
);
795 else if (SCM_BIGP (y
))
797 SCM result
= scm_i_mkbig ();
798 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
801 scm_remember_upto_here_2 (x
, y
);
802 return scm_i_normbig (result
);
805 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
808 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
811 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
812 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
814 * "(remainder 13 4) @result{} 1\n"
815 * "(remainder -13 4) @result{} -1\n"
819 scm_remainder (SCM x
, SCM y
)
825 long yy
= SCM_I_INUM (y
);
827 scm_num_overflow (s_remainder
);
830 long z
= SCM_I_INUM (x
) % yy
;
831 return SCM_I_MAKINUM (z
);
834 else if (SCM_BIGP (y
))
836 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
837 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
838 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
840 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
841 scm_remember_upto_here_1 (y
);
842 return SCM_I_MAKINUM (0);
848 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
850 else if (SCM_BIGP (x
))
854 long yy
= SCM_I_INUM (y
);
856 scm_num_overflow (s_remainder
);
859 SCM result
= scm_i_mkbig ();
862 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
863 scm_remember_upto_here_1 (x
);
864 return scm_i_normbig (result
);
867 else if (SCM_BIGP (y
))
869 SCM result
= scm_i_mkbig ();
870 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
873 scm_remember_upto_here_2 (x
, y
);
874 return scm_i_normbig (result
);
877 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
880 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
884 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
885 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
887 * "(modulo 13 4) @result{} 1\n"
888 * "(modulo -13 4) @result{} 3\n"
892 scm_modulo (SCM x
, SCM y
)
896 long xx
= SCM_I_INUM (x
);
899 long yy
= SCM_I_INUM (y
);
901 scm_num_overflow (s_modulo
);
904 /* C99 specifies that "%" is the remainder corresponding to a
905 quotient rounded towards zero, and that's also traditional
906 for machine division, so z here should be well defined. */
924 return SCM_I_MAKINUM (result
);
927 else if (SCM_BIGP (y
))
929 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
936 SCM pos_y
= scm_i_clonebig (y
, 0);
937 /* do this after the last scm_op */
938 mpz_init_set_si (z_x
, xx
);
939 result
= pos_y
; /* re-use this bignum */
940 mpz_mod (SCM_I_BIG_MPZ (result
),
942 SCM_I_BIG_MPZ (pos_y
));
943 scm_remember_upto_here_1 (pos_y
);
947 result
= scm_i_mkbig ();
948 /* do this after the last scm_op */
949 mpz_init_set_si (z_x
, xx
);
950 mpz_mod (SCM_I_BIG_MPZ (result
),
953 scm_remember_upto_here_1 (y
);
956 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
957 mpz_add (SCM_I_BIG_MPZ (result
),
959 SCM_I_BIG_MPZ (result
));
960 scm_remember_upto_here_1 (y
);
961 /* and do this before the next one */
963 return scm_i_normbig (result
);
967 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
969 else if (SCM_BIGP (x
))
973 long yy
= SCM_I_INUM (y
);
975 scm_num_overflow (s_modulo
);
978 SCM result
= scm_i_mkbig ();
979 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
981 (yy
< 0) ? - yy
: yy
);
982 scm_remember_upto_here_1 (x
);
983 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
984 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
985 SCM_I_BIG_MPZ (result
),
987 return scm_i_normbig (result
);
990 else if (SCM_BIGP (y
))
993 SCM result
= scm_i_mkbig ();
994 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
995 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
996 mpz_mod (SCM_I_BIG_MPZ (result
),
998 SCM_I_BIG_MPZ (pos_y
));
1000 scm_remember_upto_here_1 (x
);
1001 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1002 mpz_add (SCM_I_BIG_MPZ (result
),
1004 SCM_I_BIG_MPZ (result
));
1005 scm_remember_upto_here_2 (y
, pos_y
);
1006 return scm_i_normbig (result
);
1010 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1013 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1016 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
1017 /* "Return the greatest common divisor of all arguments.\n"
1018 * "If called without arguments, 0 is returned."
1021 scm_gcd (SCM x
, SCM y
)
1024 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1026 if (SCM_I_INUMP (x
))
1028 if (SCM_I_INUMP (y
))
1030 long xx
= SCM_I_INUM (x
);
1031 long yy
= SCM_I_INUM (y
);
1032 long u
= xx
< 0 ? -xx
: xx
;
1033 long v
= yy
< 0 ? -yy
: yy
;
1043 /* Determine a common factor 2^k */
1044 while (!(1 & (u
| v
)))
1050 /* Now, any factor 2^n can be eliminated */
1070 return (SCM_POSFIXABLE (result
)
1071 ? SCM_I_MAKINUM (result
)
1072 : scm_i_long2big (result
));
1074 else if (SCM_BIGP (y
))
1080 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1082 else if (SCM_BIGP (x
))
1084 if (SCM_I_INUMP (y
))
1086 unsigned long result
;
1089 yy
= SCM_I_INUM (y
);
1094 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1095 scm_remember_upto_here_1 (x
);
1096 return (SCM_POSFIXABLE (result
)
1097 ? SCM_I_MAKINUM (result
)
1098 : scm_from_ulong (result
));
1100 else if (SCM_BIGP (y
))
1102 SCM result
= scm_i_mkbig ();
1103 mpz_gcd (SCM_I_BIG_MPZ (result
),
1106 scm_remember_upto_here_2 (x
, y
);
1107 return scm_i_normbig (result
);
1110 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1113 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1116 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1117 /* "Return the least common multiple of the arguments.\n"
1118 * "If called without arguments, 1 is returned."
1121 scm_lcm (SCM n1
, SCM n2
)
1123 if (SCM_UNBNDP (n2
))
1125 if (SCM_UNBNDP (n1
))
1126 return SCM_I_MAKINUM (1L);
1127 n2
= SCM_I_MAKINUM (1L);
1130 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1131 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1132 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1133 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1135 if (SCM_I_INUMP (n1
))
1137 if (SCM_I_INUMP (n2
))
1139 SCM d
= scm_gcd (n1
, n2
);
1140 if (scm_is_eq (d
, SCM_INUM0
))
1143 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1147 /* inum n1, big n2 */
1150 SCM result
= scm_i_mkbig ();
1151 long nn1
= SCM_I_INUM (n1
);
1152 if (nn1
== 0) return SCM_INUM0
;
1153 if (nn1
< 0) nn1
= - nn1
;
1154 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1155 scm_remember_upto_here_1 (n2
);
1163 if (SCM_I_INUMP (n2
))
1170 SCM result
= scm_i_mkbig ();
1171 mpz_lcm(SCM_I_BIG_MPZ (result
),
1173 SCM_I_BIG_MPZ (n2
));
1174 scm_remember_upto_here_2(n1
, n2
);
1175 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1181 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1186 + + + x (map digit:logand X Y)
1187 + - + x (map digit:logand X (lognot (+ -1 Y)))
1188 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1189 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1194 + + + (map digit:logior X Y)
1195 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1196 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1197 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1202 + + + (map digit:logxor X Y)
1203 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1204 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1205 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1210 + + (any digit:logand X Y)
1211 + - (any digit:logand X (lognot (+ -1 Y)))
1212 - + (any digit:logand (lognot (+ -1 X)) Y)
1217 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1219 "Return the bitwise AND of the integer arguments.\n\n"
1221 "(logand) @result{} -1\n"
1222 "(logand 7) @result{} 7\n"
1223 "(logand #b111 #b011 #b001) @result{} 1\n"
1225 #define FUNC_NAME s_scm_logand
1229 if (SCM_UNBNDP (n2
))
1231 if (SCM_UNBNDP (n1
))
1232 return SCM_I_MAKINUM (-1);
1233 else if (!SCM_NUMBERP (n1
))
1234 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1235 else if (SCM_NUMBERP (n1
))
1238 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1241 if (SCM_I_INUMP (n1
))
1243 nn1
= SCM_I_INUM (n1
);
1244 if (SCM_I_INUMP (n2
))
1246 long nn2
= SCM_I_INUM (n2
);
1247 return SCM_I_MAKINUM (nn1
& nn2
);
1249 else if SCM_BIGP (n2
)
1255 SCM result_z
= scm_i_mkbig ();
1257 mpz_init_set_si (nn1_z
, nn1
);
1258 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1259 scm_remember_upto_here_1 (n2
);
1261 return scm_i_normbig (result_z
);
1265 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1267 else if (SCM_BIGP (n1
))
1269 if (SCM_I_INUMP (n2
))
1272 nn1
= SCM_I_INUM (n1
);
1275 else if (SCM_BIGP (n2
))
1277 SCM result_z
= scm_i_mkbig ();
1278 mpz_and (SCM_I_BIG_MPZ (result_z
),
1280 SCM_I_BIG_MPZ (n2
));
1281 scm_remember_upto_here_2 (n1
, n2
);
1282 return scm_i_normbig (result_z
);
1285 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1288 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1293 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1295 "Return the bitwise OR of the integer arguments.\n\n"
1297 "(logior) @result{} 0\n"
1298 "(logior 7) @result{} 7\n"
1299 "(logior #b000 #b001 #b011) @result{} 3\n"
1301 #define FUNC_NAME s_scm_logior
1305 if (SCM_UNBNDP (n2
))
1307 if (SCM_UNBNDP (n1
))
1309 else if (SCM_NUMBERP (n1
))
1312 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1315 if (SCM_I_INUMP (n1
))
1317 nn1
= SCM_I_INUM (n1
);
1318 if (SCM_I_INUMP (n2
))
1320 long nn2
= SCM_I_INUM (n2
);
1321 return SCM_I_MAKINUM (nn1
| nn2
);
1323 else if (SCM_BIGP (n2
))
1329 SCM result_z
= scm_i_mkbig ();
1331 mpz_init_set_si (nn1_z
, nn1
);
1332 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1333 scm_remember_upto_here_1 (n2
);
1335 return scm_i_normbig (result_z
);
1339 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1341 else if (SCM_BIGP (n1
))
1343 if (SCM_I_INUMP (n2
))
1346 nn1
= SCM_I_INUM (n1
);
1349 else if (SCM_BIGP (n2
))
1351 SCM result_z
= scm_i_mkbig ();
1352 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1354 SCM_I_BIG_MPZ (n2
));
1355 scm_remember_upto_here_2 (n1
, n2
);
1356 return scm_i_normbig (result_z
);
1359 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1362 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1367 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1369 "Return the bitwise XOR of the integer arguments. A bit is\n"
1370 "set in the result if it is set in an odd number of arguments.\n"
1372 "(logxor) @result{} 0\n"
1373 "(logxor 7) @result{} 7\n"
1374 "(logxor #b000 #b001 #b011) @result{} 2\n"
1375 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1377 #define FUNC_NAME s_scm_logxor
1381 if (SCM_UNBNDP (n2
))
1383 if (SCM_UNBNDP (n1
))
1385 else if (SCM_NUMBERP (n1
))
1388 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1391 if (SCM_I_INUMP (n1
))
1393 nn1
= SCM_I_INUM (n1
);
1394 if (SCM_I_INUMP (n2
))
1396 long nn2
= SCM_I_INUM (n2
);
1397 return SCM_I_MAKINUM (nn1
^ nn2
);
1399 else if (SCM_BIGP (n2
))
1403 SCM result_z
= scm_i_mkbig ();
1405 mpz_init_set_si (nn1_z
, nn1
);
1406 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1407 scm_remember_upto_here_1 (n2
);
1409 return scm_i_normbig (result_z
);
1413 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1415 else if (SCM_BIGP (n1
))
1417 if (SCM_I_INUMP (n2
))
1420 nn1
= SCM_I_INUM (n1
);
1423 else if (SCM_BIGP (n2
))
1425 SCM result_z
= scm_i_mkbig ();
1426 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1428 SCM_I_BIG_MPZ (n2
));
1429 scm_remember_upto_here_2 (n1
, n2
);
1430 return scm_i_normbig (result_z
);
1433 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1436 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1441 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1443 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1444 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1445 "without actually calculating the @code{logand}, just testing\n"
1449 "(logtest #b0100 #b1011) @result{} #f\n"
1450 "(logtest #b0100 #b0111) @result{} #t\n"
1452 #define FUNC_NAME s_scm_logtest
1456 if (SCM_I_INUMP (j
))
1458 nj
= SCM_I_INUM (j
);
1459 if (SCM_I_INUMP (k
))
1461 long nk
= SCM_I_INUM (k
);
1462 return scm_from_bool (nj
& nk
);
1464 else if (SCM_BIGP (k
))
1472 mpz_init_set_si (nj_z
, nj
);
1473 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1474 scm_remember_upto_here_1 (k
);
1475 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1481 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1483 else if (SCM_BIGP (j
))
1485 if (SCM_I_INUMP (k
))
1488 nj
= SCM_I_INUM (j
);
1491 else if (SCM_BIGP (k
))
1495 mpz_init (result_z
);
1499 scm_remember_upto_here_2 (j
, k
);
1500 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1501 mpz_clear (result_z
);
1505 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1508 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1513 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1515 "Test whether bit number @var{index} in @var{j} is set.\n"
1516 "@var{index} starts from 0 for the least significant bit.\n"
1519 "(logbit? 0 #b1101) @result{} #t\n"
1520 "(logbit? 1 #b1101) @result{} #f\n"
1521 "(logbit? 2 #b1101) @result{} #t\n"
1522 "(logbit? 3 #b1101) @result{} #t\n"
1523 "(logbit? 4 #b1101) @result{} #f\n"
1525 #define FUNC_NAME s_scm_logbit_p
1527 unsigned long int iindex
;
1528 iindex
= scm_to_ulong (index
);
1530 if (SCM_I_INUMP (j
))
1532 /* bits above what's in an inum follow the sign bit */
1533 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1534 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1536 else if (SCM_BIGP (j
))
1538 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1539 scm_remember_upto_here_1 (j
);
1540 return scm_from_bool (val
);
1543 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1548 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1550 "Return the integer which is the ones-complement of the integer\n"
1554 "(number->string (lognot #b10000000) 2)\n"
1555 " @result{} \"-10000001\"\n"
1556 "(number->string (lognot #b0) 2)\n"
1557 " @result{} \"-1\"\n"
1559 #define FUNC_NAME s_scm_lognot
1561 if (SCM_I_INUMP (n
)) {
1562 /* No overflow here, just need to toggle all the bits making up the inum.
1563 Enhancement: No need to strip the tag and add it back, could just xor
1564 a block of 1 bits, if that worked with the various debug versions of
1566 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1568 } else if (SCM_BIGP (n
)) {
1569 SCM result
= scm_i_mkbig ();
1570 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1571 scm_remember_upto_here_1 (n
);
1575 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1580 /* returns 0 if IN is not an integer. OUT must already be
1583 coerce_to_big (SCM in
, mpz_t out
)
1586 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1587 else if (SCM_I_INUMP (in
))
1588 mpz_set_si (out
, SCM_I_INUM (in
));
1595 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1596 (SCM n
, SCM k
, SCM m
),
1597 "Return @var{n} raised to the integer exponent\n"
1598 "@var{k}, modulo @var{m}.\n"
1601 "(modulo-expt 2 3 5)\n"
1604 #define FUNC_NAME s_scm_modulo_expt
1610 /* There are two classes of error we might encounter --
1611 1) Math errors, which we'll report by calling scm_num_overflow,
1613 2) wrong-type errors, which of course we'll report by calling
1615 We don't report those errors immediately, however; instead we do
1616 some cleanup first. These variables tell us which error (if
1617 any) we should report after cleaning up.
1619 int report_overflow
= 0;
1621 int position_of_wrong_type
= 0;
1622 SCM value_of_wrong_type
= SCM_INUM0
;
1624 SCM result
= SCM_UNDEFINED
;
1630 if (scm_is_eq (m
, SCM_INUM0
))
1632 report_overflow
= 1;
1636 if (!coerce_to_big (n
, n_tmp
))
1638 value_of_wrong_type
= n
;
1639 position_of_wrong_type
= 1;
1643 if (!coerce_to_big (k
, k_tmp
))
1645 value_of_wrong_type
= k
;
1646 position_of_wrong_type
= 2;
1650 if (!coerce_to_big (m
, m_tmp
))
1652 value_of_wrong_type
= m
;
1653 position_of_wrong_type
= 3;
1657 /* if the exponent K is negative, and we simply call mpz_powm, we
1658 will get a divide-by-zero exception when an inverse 1/n mod m
1659 doesn't exist (or is not unique). Since exceptions are hard to
1660 handle, we'll attempt the inversion "by hand" -- that way, we get
1661 a simple failure code, which is easy to handle. */
1663 if (-1 == mpz_sgn (k_tmp
))
1665 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1667 report_overflow
= 1;
1670 mpz_neg (k_tmp
, k_tmp
);
1673 result
= scm_i_mkbig ();
1674 mpz_powm (SCM_I_BIG_MPZ (result
),
1679 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1680 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1687 if (report_overflow
)
1688 scm_num_overflow (FUNC_NAME
);
1690 if (position_of_wrong_type
)
1691 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1692 value_of_wrong_type
);
1694 return scm_i_normbig (result
);
1698 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1700 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1701 "exact integer, @var{n} can be any number.\n"
1703 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1704 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1705 "includes @math{0^0} is 1.\n"
1708 "(integer-expt 2 5) @result{} 32\n"
1709 "(integer-expt -3 3) @result{} -27\n"
1710 "(integer-expt 5 -3) @result{} 1/125\n"
1711 "(integer-expt 0 0) @result{} 1\n"
1713 #define FUNC_NAME s_scm_integer_expt
1716 SCM z_i2
= SCM_BOOL_F
;
1718 SCM acc
= SCM_I_MAKINUM (1L);
1720 /* 0^0 == 1 according to R5RS */
1721 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1722 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1723 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1724 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1726 if (SCM_I_INUMP (k
))
1727 i2
= SCM_I_INUM (k
);
1728 else if (SCM_BIGP (k
))
1730 z_i2
= scm_i_clonebig (k
, 1);
1731 scm_remember_upto_here_1 (k
);
1735 SCM_WRONG_TYPE_ARG (2, k
);
1739 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1741 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1742 n
= scm_divide (n
, SCM_UNDEFINED
);
1746 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1750 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1752 return scm_product (acc
, n
);
1754 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1755 acc
= scm_product (acc
, n
);
1756 n
= scm_product (n
, n
);
1757 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1765 n
= scm_divide (n
, SCM_UNDEFINED
);
1772 return scm_product (acc
, n
);
1774 acc
= scm_product (acc
, n
);
1775 n
= scm_product (n
, n
);
1782 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1784 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1785 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1787 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1788 "@var{cnt} is negative it's a division, rounded towards negative\n"
1789 "infinity. (Note that this is not the same rounding as\n"
1790 "@code{quotient} does.)\n"
1792 "With @var{n} viewed as an infinite precision twos complement,\n"
1793 "@code{ash} means a left shift introducing zero bits, or a right\n"
1794 "shift dropping bits.\n"
1797 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1798 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1800 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1801 "(ash -23 -2) @result{} -6\n"
1803 #define FUNC_NAME s_scm_ash
1806 bits_to_shift
= scm_to_long (cnt
);
1808 if (SCM_I_INUMP (n
))
1810 long nn
= SCM_I_INUM (n
);
1812 if (bits_to_shift
> 0)
1814 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1815 overflow a non-zero fixnum. For smaller shifts we check the
1816 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1817 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1818 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1824 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1826 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1829 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1833 SCM result
= scm_i_long2big (nn
);
1834 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1841 bits_to_shift
= -bits_to_shift
;
1842 if (bits_to_shift
>= SCM_LONG_BIT
)
1843 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1845 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1849 else if (SCM_BIGP (n
))
1853 if (bits_to_shift
== 0)
1856 result
= scm_i_mkbig ();
1857 if (bits_to_shift
>= 0)
1859 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1865 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1866 we have to allocate a bignum even if the result is going to be a
1868 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1870 return scm_i_normbig (result
);
1876 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1882 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1883 (SCM n
, SCM start
, SCM end
),
1884 "Return the integer composed of the @var{start} (inclusive)\n"
1885 "through @var{end} (exclusive) bits of @var{n}. The\n"
1886 "@var{start}th bit becomes the 0-th bit in the result.\n"
1889 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1890 " @result{} \"1010\"\n"
1891 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1892 " @result{} \"10110\"\n"
1894 #define FUNC_NAME s_scm_bit_extract
1896 unsigned long int istart
, iend
, bits
;
1897 istart
= scm_to_ulong (start
);
1898 iend
= scm_to_ulong (end
);
1899 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1901 /* how many bits to keep */
1902 bits
= iend
- istart
;
1904 if (SCM_I_INUMP (n
))
1906 long int in
= SCM_I_INUM (n
);
1908 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1909 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1910 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1912 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1914 /* Since we emulate two's complement encoded numbers, this
1915 * special case requires us to produce a result that has
1916 * more bits than can be stored in a fixnum.
1918 SCM result
= scm_i_long2big (in
);
1919 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1924 /* mask down to requisite bits */
1925 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1926 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1928 else if (SCM_BIGP (n
))
1933 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1937 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1938 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1939 such bits into a ulong. */
1940 result
= scm_i_mkbig ();
1941 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1942 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1943 result
= scm_i_normbig (result
);
1945 scm_remember_upto_here_1 (n
);
1949 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1954 static const char scm_logtab
[] = {
1955 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1958 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1960 "Return the number of bits in integer @var{n}. If integer is\n"
1961 "positive, the 1-bits in its binary representation are counted.\n"
1962 "If negative, the 0-bits in its two's-complement binary\n"
1963 "representation are counted. If 0, 0 is returned.\n"
1966 "(logcount #b10101010)\n"
1973 #define FUNC_NAME s_scm_logcount
1975 if (SCM_I_INUMP (n
))
1977 unsigned long int c
= 0;
1978 long int nn
= SCM_I_INUM (n
);
1983 c
+= scm_logtab
[15 & nn
];
1986 return SCM_I_MAKINUM (c
);
1988 else if (SCM_BIGP (n
))
1990 unsigned long count
;
1991 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1992 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1994 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1995 scm_remember_upto_here_1 (n
);
1996 return SCM_I_MAKINUM (count
);
1999 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2004 static const char scm_ilentab
[] = {
2005 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2009 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2011 "Return the number of bits necessary to represent @var{n}.\n"
2014 "(integer-length #b10101010)\n"
2016 "(integer-length 0)\n"
2018 "(integer-length #b1111)\n"
2021 #define FUNC_NAME s_scm_integer_length
2023 if (SCM_I_INUMP (n
))
2025 unsigned long int c
= 0;
2027 long int nn
= SCM_I_INUM (n
);
2033 l
= scm_ilentab
[15 & nn
];
2036 return SCM_I_MAKINUM (c
- 4 + l
);
2038 else if (SCM_BIGP (n
))
2040 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2041 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2042 1 too big, so check for that and adjust. */
2043 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2044 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2045 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2046 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2048 scm_remember_upto_here_1 (n
);
2049 return SCM_I_MAKINUM (size
);
2052 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2056 /*** NUMBERS -> STRINGS ***/
2057 #define SCM_MAX_DBL_PREC 60
2058 #define SCM_MAX_DBL_RADIX 36
2060 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2061 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2062 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2065 void init_dblprec(int *prec
, int radix
) {
2066 /* determine floating point precision by adding successively
2067 smaller increments to 1.0 until it is considered == 1.0 */
2068 double f
= ((double)1.0)/radix
;
2069 double fsum
= 1.0 + f
;
2074 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2086 void init_fx_radix(double *fx_list
, int radix
)
2088 /* initialize a per-radix list of tolerances. When added
2089 to a number < 1.0, we can determine if we should raund
2090 up and quit converting a number to a string. */
2094 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2095 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2098 /* use this array as a way to generate a single digit */
2099 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2102 idbl2str (double f
, char *a
, int radix
)
2104 int efmt
, dpt
, d
, i
, wp
;
2106 #ifdef DBL_MIN_10_EXP
2109 #endif /* DBL_MIN_10_EXP */
2114 radix
> SCM_MAX_DBL_RADIX
)
2116 /* revert to existing behavior */
2120 wp
= scm_dblprec
[radix
-2];
2121 fx
= fx_per_radix
[radix
-2];
2125 #ifdef HAVE_COPYSIGN
2126 double sgn
= copysign (1.0, f
);
2131 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2137 strcpy (a
, "-inf.0");
2139 strcpy (a
, "+inf.0");
2142 else if (xisnan (f
))
2144 strcpy (a
, "+nan.0");
2154 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2155 make-uniform-vector, from causing infinite loops. */
2156 /* just do the checking...if it passes, we do the conversion for our
2157 radix again below */
2164 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2172 while (f_cpy
> 10.0)
2175 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2196 if (f
+ fx
[wp
] >= radix
)
2203 /* adding 9999 makes this equivalent to abs(x) % 3 */
2204 dpt
= (exp
+ 9999) % 3;
2208 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2230 a
[ch
++] = number_chars
[d
];
2233 if (f
+ fx
[wp
] >= 1.0)
2235 a
[ch
- 1] = number_chars
[d
+1];
2247 if ((dpt
> 4) && (exp
> 6))
2249 d
= (a
[0] == '-' ? 2 : 1);
2250 for (i
= ch
++; i
> d
; i
--)
2263 if (a
[ch
- 1] == '.')
2264 a
[ch
++] = '0'; /* trailing zero */
2273 for (i
= radix
; i
<= exp
; i
*= radix
);
2274 for (i
/= radix
; i
; i
/= radix
)
2276 a
[ch
++] = number_chars
[exp
/ i
];
2285 icmplx2str (double real
, double imag
, char *str
, int radix
)
2289 i
= idbl2str (real
, str
, radix
);
2292 /* Don't output a '+' for negative numbers or for Inf and
2293 NaN. They will provide their own sign. */
2294 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2296 i
+= idbl2str (imag
, &str
[i
], radix
);
2303 iflo2str (SCM flt
, char *str
, int radix
)
2306 if (SCM_REALP (flt
))
2307 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2309 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2314 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2315 characters in the result.
2317 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2319 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2324 return scm_iuint2str (-num
, rad
, p
) + 1;
2327 return scm_iuint2str (num
, rad
, p
);
2330 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2331 characters in the result.
2333 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2335 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2339 scm_t_uintmax n
= num
;
2341 for (n
/= rad
; n
> 0; n
/= rad
)
2351 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2356 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2358 "Return a string holding the external representation of the\n"
2359 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2360 "inexact, a radix of 10 will be used.")
2361 #define FUNC_NAME s_scm_number_to_string
2365 if (SCM_UNBNDP (radix
))
2368 base
= scm_to_signed_integer (radix
, 2, 36);
2370 if (SCM_I_INUMP (n
))
2372 char num_buf
[SCM_INTBUFLEN
];
2373 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2374 return scm_from_locale_stringn (num_buf
, length
);
2376 else if (SCM_BIGP (n
))
2378 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2379 scm_remember_upto_here_1 (n
);
2380 return scm_take_locale_string (str
);
2382 else if (SCM_FRACTIONP (n
))
2384 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2385 scm_from_locale_string ("/"),
2386 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2388 else if (SCM_INEXACTP (n
))
2390 char num_buf
[FLOBUFLEN
];
2391 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2394 SCM_WRONG_TYPE_ARG (1, n
);
2399 /* These print routines used to be stubbed here so that scm_repl.c
2400 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2403 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2405 char num_buf
[FLOBUFLEN
];
2406 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2411 scm_i_print_double (double val
, SCM port
)
2413 char num_buf
[FLOBUFLEN
];
2414 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2418 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2421 char num_buf
[FLOBUFLEN
];
2422 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2427 scm_i_print_complex (double real
, double imag
, SCM port
)
2429 char num_buf
[FLOBUFLEN
];
2430 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2434 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2437 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2438 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2439 scm_remember_upto_here_1 (str
);
2444 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2446 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2447 scm_remember_upto_here_1 (exp
);
2448 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2452 /*** END nums->strs ***/
2455 /*** STRINGS -> NUMBERS ***/
2457 /* The following functions implement the conversion from strings to numbers.
2458 * The implementation somehow follows the grammar for numbers as it is given
2459 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2460 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2461 * points should be noted about the implementation:
2462 * * Each function keeps a local index variable 'idx' that points at the
2463 * current position within the parsed string. The global index is only
2464 * updated if the function could parse the corresponding syntactic unit
2466 * * Similarly, the functions keep track of indicators of inexactness ('#',
2467 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2468 * global exactness information is only updated after each part has been
2469 * successfully parsed.
2470 * * Sequences of digits are parsed into temporary variables holding fixnums.
2471 * Only if these fixnums would overflow, the result variables are updated
2472 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2473 * the temporary variables holding the fixnums are cleared, and the process
2474 * starts over again. If for example fixnums were able to store five decimal
2475 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2476 * and the result was computed as 12345 * 100000 + 67890. In other words,
2477 * only every five digits two bignum operations were performed.
2480 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2482 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2484 /* In non ASCII-style encodings the following macro might not work. */
2485 #define XDIGIT2UINT(d) \
2486 (isdigit ((int) (unsigned char) d) \
2488 : tolower ((int) (unsigned char) d) - 'a' + 10)
2491 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2492 unsigned int radix
, enum t_exactness
*p_exactness
)
2494 unsigned int idx
= *p_idx
;
2495 unsigned int hash_seen
= 0;
2496 scm_t_bits shift
= 1;
2498 unsigned int digit_value
;
2506 if (!isxdigit ((int) (unsigned char) c
))
2508 digit_value
= XDIGIT2UINT (c
);
2509 if (digit_value
>= radix
)
2513 result
= SCM_I_MAKINUM (digit_value
);
2517 if (isxdigit ((int) (unsigned char) c
))
2521 digit_value
= XDIGIT2UINT (c
);
2522 if (digit_value
>= radix
)
2534 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2536 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2538 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2545 shift
= shift
* radix
;
2546 add
= add
* radix
+ digit_value
;
2551 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2553 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2557 *p_exactness
= INEXACT
;
2563 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2564 * covers the parts of the rules that start at a potential point. The value
2565 * of the digits up to the point have been parsed by the caller and are given
2566 * in variable result. The content of *p_exactness indicates, whether a hash
2567 * has already been seen in the digits before the point.
2570 /* In non ASCII-style encodings the following macro might not work. */
2571 #define DIGIT2UINT(d) ((d) - '0')
2574 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2575 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2577 unsigned int idx
= *p_idx
;
2578 enum t_exactness x
= *p_exactness
;
2583 if (mem
[idx
] == '.')
2585 scm_t_bits shift
= 1;
2587 unsigned int digit_value
;
2588 SCM big_shift
= SCM_I_MAKINUM (1);
2594 if (isdigit ((int) (unsigned char) c
))
2599 digit_value
= DIGIT2UINT (c
);
2610 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2612 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2613 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2615 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2623 add
= add
* 10 + digit_value
;
2629 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2630 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2631 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2634 result
= scm_divide (result
, big_shift
);
2636 /* We've seen a decimal point, thus the value is implicitly inexact. */
2648 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2675 if (!isdigit ((int) (unsigned char) c
))
2679 exponent
= DIGIT2UINT (c
);
2683 if (isdigit ((int) (unsigned char) c
))
2686 if (exponent
<= SCM_MAXEXP
)
2687 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2693 if (exponent
> SCM_MAXEXP
)
2695 size_t exp_len
= idx
- start
;
2696 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2697 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2698 scm_out_of_range ("string->number", exp_num
);
2701 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2703 result
= scm_product (result
, e
);
2705 result
= scm_divide2real (result
, e
);
2707 /* We've seen an exponent, thus the value is implicitly inexact. */
2725 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2728 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2729 unsigned int radix
, enum t_exactness
*p_exactness
)
2731 unsigned int idx
= *p_idx
;
2737 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2743 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2745 enum t_exactness x
= EXACT
;
2747 /* Cobble up the fractional part. We might want to set the
2748 NaN's mantissa from it. */
2750 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2755 if (mem
[idx
] == '.')
2759 else if (idx
+ 1 == len
)
2761 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2764 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2765 p_idx
, p_exactness
);
2769 enum t_exactness x
= EXACT
;
2772 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2773 if (scm_is_false (uinteger
))
2778 else if (mem
[idx
] == '/')
2784 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2785 if (scm_is_false (divisor
))
2788 /* both are int/big here, I assume */
2789 result
= scm_i_make_ratio (uinteger
, divisor
);
2791 else if (radix
== 10)
2793 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2794 if (scm_is_false (result
))
2805 /* When returning an inexact zero, make sure it is represented as a
2806 floating point value so that we can change its sign.
2808 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2809 result
= scm_from_double (0.0);
2815 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2818 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2819 unsigned int radix
, enum t_exactness
*p_exactness
)
2843 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2844 if (scm_is_false (ureal
))
2846 /* input must be either +i or -i */
2851 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2857 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2864 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2865 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2874 /* either +<ureal>i or -<ureal>i */
2881 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2884 /* polar input: <real>@<real>. */
2909 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2910 if (scm_is_false (angle
))
2915 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2916 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2918 result
= scm_make_polar (ureal
, angle
);
2923 /* expecting input matching <real>[+-]<ureal>?i */
2930 int sign
= (c
== '+') ? 1 : -1;
2931 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2933 if (scm_is_false (imag
))
2934 imag
= SCM_I_MAKINUM (sign
);
2935 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2936 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2940 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2947 return scm_make_rectangular (ureal
, imag
);
2956 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2958 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2961 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2962 unsigned int default_radix
)
2964 unsigned int idx
= 0;
2965 unsigned int radix
= NO_RADIX
;
2966 enum t_exactness forced_x
= NO_EXACTNESS
;
2967 enum t_exactness implicit_x
= EXACT
;
2970 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2971 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2973 switch (mem
[idx
+ 1])
2976 if (radix
!= NO_RADIX
)
2981 if (radix
!= NO_RADIX
)
2986 if (forced_x
!= NO_EXACTNESS
)
2991 if (forced_x
!= NO_EXACTNESS
)
2996 if (radix
!= NO_RADIX
)
3001 if (radix
!= NO_RADIX
)
3011 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3012 if (radix
== NO_RADIX
)
3013 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
3015 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
3017 if (scm_is_false (result
))
3023 if (SCM_INEXACTP (result
))
3024 return scm_inexact_to_exact (result
);
3028 if (SCM_INEXACTP (result
))
3031 return scm_exact_to_inexact (result
);
3034 if (implicit_x
== INEXACT
)
3036 if (SCM_INEXACTP (result
))
3039 return scm_exact_to_inexact (result
);
3047 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3048 (SCM string
, SCM radix
),
3049 "Return a number of the maximally precise representation\n"
3050 "expressed by the given @var{string}. @var{radix} must be an\n"
3051 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3052 "is a default radix that may be overridden by an explicit radix\n"
3053 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3054 "supplied, then the default radix is 10. If string is not a\n"
3055 "syntactically valid notation for a number, then\n"
3056 "@code{string->number} returns @code{#f}.")
3057 #define FUNC_NAME s_scm_string_to_number
3061 SCM_VALIDATE_STRING (1, string
);
3063 if (SCM_UNBNDP (radix
))
3066 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3068 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3069 scm_i_string_length (string
),
3071 scm_remember_upto_here_1 (string
);
3077 /*** END strs->nums ***/
3081 scm_bigequal (SCM x
, SCM y
)
3083 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3084 scm_remember_upto_here_2 (x
, y
);
3085 return scm_from_bool (0 == result
);
3089 scm_real_equalp (SCM x
, SCM y
)
3091 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3095 scm_complex_equalp (SCM x
, SCM y
)
3097 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3098 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3102 scm_i_fraction_equalp (SCM x
, SCM y
)
3104 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3105 SCM_FRACTION_NUMERATOR (y
)))
3106 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3107 SCM_FRACTION_DENOMINATOR (y
))))
3114 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3116 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3118 #define FUNC_NAME s_scm_number_p
3120 return scm_from_bool (SCM_NUMBERP (x
));
3124 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3126 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3127 "otherwise. Note that the sets of real, rational and integer\n"
3128 "values form subsets of the set of complex numbers, i. e. the\n"
3129 "predicate will also be fulfilled if @var{x} is a real,\n"
3130 "rational or integer number.")
3131 #define FUNC_NAME s_scm_complex_p
3133 /* all numbers are complex. */
3134 return scm_number_p (x
);
3138 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3140 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3141 "otherwise. Note that the set of integer values forms a subset of\n"
3142 "the set of real numbers, i. e. the predicate will also be\n"
3143 "fulfilled if @var{x} is an integer number.")
3144 #define FUNC_NAME s_scm_real_p
3146 /* we can't represent irrational numbers. */
3147 return scm_rational_p (x
);
3151 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3153 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3154 "otherwise. Note that the set of integer values forms a subset of\n"
3155 "the set of rational numbers, i. e. the predicate will also be\n"
3156 "fulfilled if @var{x} is an integer number.")
3157 #define FUNC_NAME s_scm_rational_p
3159 if (SCM_I_INUMP (x
))
3161 else if (SCM_IMP (x
))
3163 else if (SCM_BIGP (x
))
3165 else if (SCM_FRACTIONP (x
))
3167 else if (SCM_REALP (x
))
3168 /* due to their limited precision, all floating point numbers are
3169 rational as well. */
3176 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3178 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3180 #define FUNC_NAME s_scm_integer_p
3183 if (SCM_I_INUMP (x
))
3189 if (!SCM_INEXACTP (x
))
3191 if (SCM_COMPLEXP (x
))
3193 r
= SCM_REAL_VALUE (x
);
3194 /* +/-inf passes r==floor(r), making those #t */
3202 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3204 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3206 #define FUNC_NAME s_scm_inexact_p
3208 if (SCM_INEXACTP (x
))
3210 if (SCM_NUMBERP (x
))
3212 SCM_WRONG_TYPE_ARG (1, x
);
3217 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3218 /* "Return @code{#t} if all parameters are numerically equal." */
3220 scm_num_eq_p (SCM x
, SCM y
)
3223 if (SCM_I_INUMP (x
))
3225 long xx
= SCM_I_INUM (x
);
3226 if (SCM_I_INUMP (y
))
3228 long yy
= SCM_I_INUM (y
);
3229 return scm_from_bool (xx
== yy
);
3231 else if (SCM_BIGP (y
))
3233 else if (SCM_REALP (y
))
3235 /* On a 32-bit system an inum fits a double, we can cast the inum
3236 to a double and compare.
3238 But on a 64-bit system an inum is bigger than a double and
3239 casting it to a double (call that dxx) will round. dxx is at
3240 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3241 an integer and fits a long. So we cast yy to a long and
3242 compare with plain xx.
3244 An alternative (for any size system actually) would be to check
3245 yy is an integer (with floor) and is in range of an inum
3246 (compare against appropriate powers of 2) then test
3247 xx==(long)yy. It's just a matter of which casts/comparisons
3248 might be fastest or easiest for the cpu. */
3250 double yy
= SCM_REAL_VALUE (y
);
3251 return scm_from_bool ((double) xx
== yy
3252 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3253 || xx
== (long) yy
));
3255 else if (SCM_COMPLEXP (y
))
3256 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3257 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3258 else if (SCM_FRACTIONP (y
))
3261 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3263 else if (SCM_BIGP (x
))
3265 if (SCM_I_INUMP (y
))
3267 else if (SCM_BIGP (y
))
3269 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3270 scm_remember_upto_here_2 (x
, y
);
3271 return scm_from_bool (0 == cmp
);
3273 else if (SCM_REALP (y
))
3276 if (xisnan (SCM_REAL_VALUE (y
)))
3278 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3279 scm_remember_upto_here_1 (x
);
3280 return scm_from_bool (0 == cmp
);
3282 else if (SCM_COMPLEXP (y
))
3285 if (0.0 != SCM_COMPLEX_IMAG (y
))
3287 if (xisnan (SCM_COMPLEX_REAL (y
)))
3289 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3290 scm_remember_upto_here_1 (x
);
3291 return scm_from_bool (0 == cmp
);
3293 else if (SCM_FRACTIONP (y
))
3296 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3298 else if (SCM_REALP (x
))
3300 double xx
= SCM_REAL_VALUE (x
);
3301 if (SCM_I_INUMP (y
))
3303 /* see comments with inum/real above */
3304 long yy
= SCM_I_INUM (y
);
3305 return scm_from_bool (xx
== (double) yy
3306 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3307 || (long) xx
== yy
));
3309 else if (SCM_BIGP (y
))
3312 if (xisnan (SCM_REAL_VALUE (x
)))
3314 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3315 scm_remember_upto_here_1 (y
);
3316 return scm_from_bool (0 == cmp
);
3318 else if (SCM_REALP (y
))
3319 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3320 else if (SCM_COMPLEXP (y
))
3321 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3322 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3323 else if (SCM_FRACTIONP (y
))
3325 double xx
= SCM_REAL_VALUE (x
);
3329 return scm_from_bool (xx
< 0.0);
3330 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3334 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3336 else if (SCM_COMPLEXP (x
))
3338 if (SCM_I_INUMP (y
))
3339 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3340 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3341 else if (SCM_BIGP (y
))
3344 if (0.0 != SCM_COMPLEX_IMAG (x
))
3346 if (xisnan (SCM_COMPLEX_REAL (x
)))
3348 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3349 scm_remember_upto_here_1 (y
);
3350 return scm_from_bool (0 == cmp
);
3352 else if (SCM_REALP (y
))
3353 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3354 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3355 else if (SCM_COMPLEXP (y
))
3356 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3357 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3358 else if (SCM_FRACTIONP (y
))
3361 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3363 xx
= SCM_COMPLEX_REAL (x
);
3367 return scm_from_bool (xx
< 0.0);
3368 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3372 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3374 else if (SCM_FRACTIONP (x
))
3376 if (SCM_I_INUMP (y
))
3378 else if (SCM_BIGP (y
))
3380 else if (SCM_REALP (y
))
3382 double yy
= SCM_REAL_VALUE (y
);
3386 return scm_from_bool (0.0 < yy
);
3387 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3390 else if (SCM_COMPLEXP (y
))
3393 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3395 yy
= SCM_COMPLEX_REAL (y
);
3399 return scm_from_bool (0.0 < yy
);
3400 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3403 else if (SCM_FRACTIONP (y
))
3404 return scm_i_fraction_equalp (x
, y
);
3406 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3409 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3413 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3414 done are good for inums, but for bignums an answer can almost always be
3415 had by just examining a few high bits of the operands, as done by GMP in
3416 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3417 of the float exponent to take into account. */
3419 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3420 /* "Return @code{#t} if the list of parameters is monotonically\n"
3424 scm_less_p (SCM x
, SCM y
)
3427 if (SCM_I_INUMP (x
))
3429 long xx
= SCM_I_INUM (x
);
3430 if (SCM_I_INUMP (y
))
3432 long yy
= SCM_I_INUM (y
);
3433 return scm_from_bool (xx
< yy
);
3435 else if (SCM_BIGP (y
))
3437 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3438 scm_remember_upto_here_1 (y
);
3439 return scm_from_bool (sgn
> 0);
3441 else if (SCM_REALP (y
))
3442 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3443 else if (SCM_FRACTIONP (y
))
3445 /* "x < a/b" becomes "x*b < a" */
3447 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3448 y
= SCM_FRACTION_NUMERATOR (y
);
3452 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3454 else if (SCM_BIGP (x
))
3456 if (SCM_I_INUMP (y
))
3458 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3459 scm_remember_upto_here_1 (x
);
3460 return scm_from_bool (sgn
< 0);
3462 else if (SCM_BIGP (y
))
3464 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3465 scm_remember_upto_here_2 (x
, y
);
3466 return scm_from_bool (cmp
< 0);
3468 else if (SCM_REALP (y
))
3471 if (xisnan (SCM_REAL_VALUE (y
)))
3473 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3474 scm_remember_upto_here_1 (x
);
3475 return scm_from_bool (cmp
< 0);
3477 else if (SCM_FRACTIONP (y
))
3480 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3482 else if (SCM_REALP (x
))
3484 if (SCM_I_INUMP (y
))
3485 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3486 else if (SCM_BIGP (y
))
3489 if (xisnan (SCM_REAL_VALUE (x
)))
3491 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3492 scm_remember_upto_here_1 (y
);
3493 return scm_from_bool (cmp
> 0);
3495 else if (SCM_REALP (y
))
3496 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3497 else if (SCM_FRACTIONP (y
))
3499 double xx
= SCM_REAL_VALUE (x
);
3503 return scm_from_bool (xx
< 0.0);
3504 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3508 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3510 else if (SCM_FRACTIONP (x
))
3512 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3514 /* "a/b < y" becomes "a < y*b" */
3515 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3516 x
= SCM_FRACTION_NUMERATOR (x
);
3519 else if (SCM_REALP (y
))
3521 double yy
= SCM_REAL_VALUE (y
);
3525 return scm_from_bool (0.0 < yy
);
3526 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3529 else if (SCM_FRACTIONP (y
))
3531 /* "a/b < c/d" becomes "a*d < c*b" */
3532 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3533 SCM_FRACTION_DENOMINATOR (y
));
3534 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3535 SCM_FRACTION_DENOMINATOR (x
));
3541 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3544 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3548 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3549 /* "Return @code{#t} if the list of parameters is monotonically\n"
3552 #define FUNC_NAME s_scm_gr_p
3554 scm_gr_p (SCM x
, SCM y
)
3556 if (!SCM_NUMBERP (x
))
3557 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3558 else if (!SCM_NUMBERP (y
))
3559 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3561 return scm_less_p (y
, x
);
3566 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3567 /* "Return @code{#t} if the list of parameters is monotonically\n"
3570 #define FUNC_NAME s_scm_leq_p
3572 scm_leq_p (SCM x
, SCM y
)
3574 if (!SCM_NUMBERP (x
))
3575 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3576 else if (!SCM_NUMBERP (y
))
3577 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3578 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3581 return scm_not (scm_less_p (y
, x
));
3586 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3587 /* "Return @code{#t} if the list of parameters is monotonically\n"
3590 #define FUNC_NAME s_scm_geq_p
3592 scm_geq_p (SCM x
, SCM y
)
3594 if (!SCM_NUMBERP (x
))
3595 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3596 else if (!SCM_NUMBERP (y
))
3597 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3598 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3601 return scm_not (scm_less_p (x
, y
));
3606 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3607 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3613 if (SCM_I_INUMP (z
))
3614 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3615 else if (SCM_BIGP (z
))
3617 else if (SCM_REALP (z
))
3618 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3619 else if (SCM_COMPLEXP (z
))
3620 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3621 && SCM_COMPLEX_IMAG (z
) == 0.0);
3622 else if (SCM_FRACTIONP (z
))
3625 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3629 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3630 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3634 scm_positive_p (SCM x
)
3636 if (SCM_I_INUMP (x
))
3637 return scm_from_bool (SCM_I_INUM (x
) > 0);
3638 else if (SCM_BIGP (x
))
3640 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3641 scm_remember_upto_here_1 (x
);
3642 return scm_from_bool (sgn
> 0);
3644 else if (SCM_REALP (x
))
3645 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3646 else if (SCM_FRACTIONP (x
))
3647 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3649 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3653 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3654 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3658 scm_negative_p (SCM x
)
3660 if (SCM_I_INUMP (x
))
3661 return scm_from_bool (SCM_I_INUM (x
) < 0);
3662 else if (SCM_BIGP (x
))
3664 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3665 scm_remember_upto_here_1 (x
);
3666 return scm_from_bool (sgn
< 0);
3668 else if (SCM_REALP (x
))
3669 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3670 else if (SCM_FRACTIONP (x
))
3671 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3673 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3677 /* scm_min and scm_max return an inexact when either argument is inexact, as
3678 required by r5rs. On that basis, for exact/inexact combinations the
3679 exact is converted to inexact to compare and possibly return. This is
3680 unlike scm_less_p above which takes some trouble to preserve all bits in
3681 its test, such trouble is not required for min and max. */
3683 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3684 /* "Return the maximum of all parameter values."
3687 scm_max (SCM x
, SCM y
)
3692 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3693 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3696 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3699 if (SCM_I_INUMP (x
))
3701 long xx
= SCM_I_INUM (x
);
3702 if (SCM_I_INUMP (y
))
3704 long yy
= SCM_I_INUM (y
);
3705 return (xx
< yy
) ? y
: x
;
3707 else if (SCM_BIGP (y
))
3709 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3710 scm_remember_upto_here_1 (y
);
3711 return (sgn
< 0) ? x
: y
;
3713 else if (SCM_REALP (y
))
3716 /* if y==NaN then ">" is false and we return NaN */
3717 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3719 else if (SCM_FRACTIONP (y
))
3722 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3725 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3727 else if (SCM_BIGP (x
))
3729 if (SCM_I_INUMP (y
))
3731 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3732 scm_remember_upto_here_1 (x
);
3733 return (sgn
< 0) ? y
: x
;
3735 else if (SCM_BIGP (y
))
3737 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3738 scm_remember_upto_here_2 (x
, y
);
3739 return (cmp
> 0) ? x
: y
;
3741 else if (SCM_REALP (y
))
3743 /* if y==NaN then xx>yy is false, so we return the NaN y */
3746 xx
= scm_i_big2dbl (x
);
3747 yy
= SCM_REAL_VALUE (y
);
3748 return (xx
> yy
? scm_from_double (xx
) : y
);
3750 else if (SCM_FRACTIONP (y
))
3755 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3757 else if (SCM_REALP (x
))
3759 if (SCM_I_INUMP (y
))
3761 double z
= SCM_I_INUM (y
);
3762 /* if x==NaN then "<" is false and we return NaN */
3763 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3765 else if (SCM_BIGP (y
))
3770 else if (SCM_REALP (y
))
3772 /* if x==NaN then our explicit check means we return NaN
3773 if y==NaN then ">" is false and we return NaN
3774 calling isnan is unavoidable, since it's the only way to know
3775 which of x or y causes any compares to be false */
3776 double xx
= SCM_REAL_VALUE (x
);
3777 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3779 else if (SCM_FRACTIONP (y
))
3781 double yy
= scm_i_fraction2double (y
);
3782 double xx
= SCM_REAL_VALUE (x
);
3783 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3786 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3788 else if (SCM_FRACTIONP (x
))
3790 if (SCM_I_INUMP (y
))
3794 else if (SCM_BIGP (y
))
3798 else if (SCM_REALP (y
))
3800 double xx
= scm_i_fraction2double (x
);
3801 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3803 else if (SCM_FRACTIONP (y
))
3808 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3811 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3815 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3816 /* "Return the minium of all parameter values."
3819 scm_min (SCM x
, SCM y
)
3824 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3825 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3828 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3831 if (SCM_I_INUMP (x
))
3833 long xx
= SCM_I_INUM (x
);
3834 if (SCM_I_INUMP (y
))
3836 long yy
= SCM_I_INUM (y
);
3837 return (xx
< yy
) ? x
: y
;
3839 else if (SCM_BIGP (y
))
3841 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3842 scm_remember_upto_here_1 (y
);
3843 return (sgn
< 0) ? y
: x
;
3845 else if (SCM_REALP (y
))
3848 /* if y==NaN then "<" is false and we return NaN */
3849 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3851 else if (SCM_FRACTIONP (y
))
3854 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3857 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3859 else if (SCM_BIGP (x
))
3861 if (SCM_I_INUMP (y
))
3863 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3864 scm_remember_upto_here_1 (x
);
3865 return (sgn
< 0) ? x
: y
;
3867 else if (SCM_BIGP (y
))
3869 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3870 scm_remember_upto_here_2 (x
, y
);
3871 return (cmp
> 0) ? y
: x
;
3873 else if (SCM_REALP (y
))
3875 /* if y==NaN then xx<yy is false, so we return the NaN y */
3878 xx
= scm_i_big2dbl (x
);
3879 yy
= SCM_REAL_VALUE (y
);
3880 return (xx
< yy
? scm_from_double (xx
) : y
);
3882 else if (SCM_FRACTIONP (y
))
3887 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3889 else if (SCM_REALP (x
))
3891 if (SCM_I_INUMP (y
))
3893 double z
= SCM_I_INUM (y
);
3894 /* if x==NaN then "<" is false and we return NaN */
3895 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3897 else if (SCM_BIGP (y
))
3902 else if (SCM_REALP (y
))
3904 /* if x==NaN then our explicit check means we return NaN
3905 if y==NaN then "<" is false and we return NaN
3906 calling isnan is unavoidable, since it's the only way to know
3907 which of x or y causes any compares to be false */
3908 double xx
= SCM_REAL_VALUE (x
);
3909 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3911 else if (SCM_FRACTIONP (y
))
3913 double yy
= scm_i_fraction2double (y
);
3914 double xx
= SCM_REAL_VALUE (x
);
3915 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3918 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3920 else if (SCM_FRACTIONP (x
))
3922 if (SCM_I_INUMP (y
))
3926 else if (SCM_BIGP (y
))
3930 else if (SCM_REALP (y
))
3932 double xx
= scm_i_fraction2double (x
);
3933 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3935 else if (SCM_FRACTIONP (y
))
3940 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3943 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3947 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3948 /* "Return the sum of all parameter values. Return 0 if called without\n"
3952 scm_sum (SCM x
, SCM y
)
3954 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3956 if (SCM_NUMBERP (x
)) return x
;
3957 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3958 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3961 if (SCM_LIKELY (SCM_I_INUMP (x
)))
3963 if (SCM_LIKELY (SCM_I_INUMP (y
)))
3965 long xx
= SCM_I_INUM (x
);
3966 long yy
= SCM_I_INUM (y
);
3967 long int z
= xx
+ yy
;
3968 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3970 else if (SCM_BIGP (y
))
3975 else if (SCM_REALP (y
))
3977 long int xx
= SCM_I_INUM (x
);
3978 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3980 else if (SCM_COMPLEXP (y
))
3982 long int xx
= SCM_I_INUM (x
);
3983 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3984 SCM_COMPLEX_IMAG (y
));
3986 else if (SCM_FRACTIONP (y
))
3987 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3988 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3989 SCM_FRACTION_DENOMINATOR (y
));
3991 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3992 } else if (SCM_BIGP (x
))
3994 if (SCM_I_INUMP (y
))
3999 inum
= SCM_I_INUM (y
);
4002 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4005 SCM result
= scm_i_mkbig ();
4006 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4007 scm_remember_upto_here_1 (x
);
4008 /* we know the result will have to be a bignum */
4011 return scm_i_normbig (result
);
4015 SCM result
= scm_i_mkbig ();
4016 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4017 scm_remember_upto_here_1 (x
);
4018 /* we know the result will have to be a bignum */
4021 return scm_i_normbig (result
);
4024 else if (SCM_BIGP (y
))
4026 SCM result
= scm_i_mkbig ();
4027 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4028 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4029 mpz_add (SCM_I_BIG_MPZ (result
),
4032 scm_remember_upto_here_2 (x
, y
);
4033 /* we know the result will have to be a bignum */
4036 return scm_i_normbig (result
);
4038 else if (SCM_REALP (y
))
4040 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4041 scm_remember_upto_here_1 (x
);
4042 return scm_from_double (result
);
4044 else if (SCM_COMPLEXP (y
))
4046 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4047 + SCM_COMPLEX_REAL (y
));
4048 scm_remember_upto_here_1 (x
);
4049 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4051 else if (SCM_FRACTIONP (y
))
4052 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4053 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4054 SCM_FRACTION_DENOMINATOR (y
));
4056 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4058 else if (SCM_REALP (x
))
4060 if (SCM_I_INUMP (y
))
4061 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4062 else if (SCM_BIGP (y
))
4064 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4065 scm_remember_upto_here_1 (y
);
4066 return scm_from_double (result
);
4068 else if (SCM_REALP (y
))
4069 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4070 else if (SCM_COMPLEXP (y
))
4071 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4072 SCM_COMPLEX_IMAG (y
));
4073 else if (SCM_FRACTIONP (y
))
4074 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4076 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4078 else if (SCM_COMPLEXP (x
))
4080 if (SCM_I_INUMP (y
))
4081 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4082 SCM_COMPLEX_IMAG (x
));
4083 else if (SCM_BIGP (y
))
4085 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4086 + SCM_COMPLEX_REAL (x
));
4087 scm_remember_upto_here_1 (y
);
4088 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4090 else if (SCM_REALP (y
))
4091 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4092 SCM_COMPLEX_IMAG (x
));
4093 else if (SCM_COMPLEXP (y
))
4094 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4095 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4096 else if (SCM_FRACTIONP (y
))
4097 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4098 SCM_COMPLEX_IMAG (x
));
4100 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4102 else if (SCM_FRACTIONP (x
))
4104 if (SCM_I_INUMP (y
))
4105 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4106 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4107 SCM_FRACTION_DENOMINATOR (x
));
4108 else if (SCM_BIGP (y
))
4109 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4110 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4111 SCM_FRACTION_DENOMINATOR (x
));
4112 else if (SCM_REALP (y
))
4113 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4114 else if (SCM_COMPLEXP (y
))
4115 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4116 SCM_COMPLEX_IMAG (y
));
4117 else if (SCM_FRACTIONP (y
))
4118 /* a/b + c/d = (ad + bc) / bd */
4119 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4120 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4121 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4123 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4126 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4130 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4132 "Return @math{@var{x}+1}.")
4133 #define FUNC_NAME s_scm_oneplus
4135 return scm_sum (x
, SCM_I_MAKINUM (1));
4140 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4141 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4142 * the sum of all but the first argument are subtracted from the first
4144 #define FUNC_NAME s_difference
4146 scm_difference (SCM x
, SCM y
)
4148 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4151 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4153 if (SCM_I_INUMP (x
))
4155 long xx
= -SCM_I_INUM (x
);
4156 if (SCM_FIXABLE (xx
))
4157 return SCM_I_MAKINUM (xx
);
4159 return scm_i_long2big (xx
);
4161 else if (SCM_BIGP (x
))
4162 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4163 bignum, but negating that gives a fixnum. */
4164 return scm_i_normbig (scm_i_clonebig (x
, 0));
4165 else if (SCM_REALP (x
))
4166 return scm_from_double (-SCM_REAL_VALUE (x
));
4167 else if (SCM_COMPLEXP (x
))
4168 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4169 -SCM_COMPLEX_IMAG (x
));
4170 else if (SCM_FRACTIONP (x
))
4171 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4172 SCM_FRACTION_DENOMINATOR (x
));
4174 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4177 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4179 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4181 long int xx
= SCM_I_INUM (x
);
4182 long int yy
= SCM_I_INUM (y
);
4183 long int z
= xx
- yy
;
4184 if (SCM_FIXABLE (z
))
4185 return SCM_I_MAKINUM (z
);
4187 return scm_i_long2big (z
);
4189 else if (SCM_BIGP (y
))
4191 /* inum-x - big-y */
4192 long xx
= SCM_I_INUM (x
);
4195 return scm_i_clonebig (y
, 0);
4198 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4199 SCM result
= scm_i_mkbig ();
4202 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4205 /* x - y == -(y + -x) */
4206 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4207 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4209 scm_remember_upto_here_1 (y
);
4211 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4212 /* we know the result will have to be a bignum */
4215 return scm_i_normbig (result
);
4218 else if (SCM_REALP (y
))
4220 long int xx
= SCM_I_INUM (x
);
4221 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4223 else if (SCM_COMPLEXP (y
))
4225 long int xx
= SCM_I_INUM (x
);
4226 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4227 - SCM_COMPLEX_IMAG (y
));
4229 else if (SCM_FRACTIONP (y
))
4230 /* a - b/c = (ac - b) / c */
4231 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4232 SCM_FRACTION_NUMERATOR (y
)),
4233 SCM_FRACTION_DENOMINATOR (y
));
4235 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4237 else if (SCM_BIGP (x
))
4239 if (SCM_I_INUMP (y
))
4241 /* big-x - inum-y */
4242 long yy
= SCM_I_INUM (y
);
4243 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4245 scm_remember_upto_here_1 (x
);
4247 return (SCM_FIXABLE (-yy
) ?
4248 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4251 SCM result
= scm_i_mkbig ();
4254 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4256 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4257 scm_remember_upto_here_1 (x
);
4259 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4260 /* we know the result will have to be a bignum */
4263 return scm_i_normbig (result
);
4266 else if (SCM_BIGP (y
))
4268 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4269 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4270 SCM result
= scm_i_mkbig ();
4271 mpz_sub (SCM_I_BIG_MPZ (result
),
4274 scm_remember_upto_here_2 (x
, y
);
4275 /* we know the result will have to be a bignum */
4276 if ((sgn_x
== 1) && (sgn_y
== -1))
4278 if ((sgn_x
== -1) && (sgn_y
== 1))
4280 return scm_i_normbig (result
);
4282 else if (SCM_REALP (y
))
4284 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4285 scm_remember_upto_here_1 (x
);
4286 return scm_from_double (result
);
4288 else if (SCM_COMPLEXP (y
))
4290 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4291 - SCM_COMPLEX_REAL (y
));
4292 scm_remember_upto_here_1 (x
);
4293 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4295 else if (SCM_FRACTIONP (y
))
4296 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4297 SCM_FRACTION_NUMERATOR (y
)),
4298 SCM_FRACTION_DENOMINATOR (y
));
4299 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4301 else if (SCM_REALP (x
))
4303 if (SCM_I_INUMP (y
))
4304 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4305 else if (SCM_BIGP (y
))
4307 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4308 scm_remember_upto_here_1 (x
);
4309 return scm_from_double (result
);
4311 else if (SCM_REALP (y
))
4312 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4313 else if (SCM_COMPLEXP (y
))
4314 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4315 -SCM_COMPLEX_IMAG (y
));
4316 else if (SCM_FRACTIONP (y
))
4317 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4319 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4321 else if (SCM_COMPLEXP (x
))
4323 if (SCM_I_INUMP (y
))
4324 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4325 SCM_COMPLEX_IMAG (x
));
4326 else if (SCM_BIGP (y
))
4328 double real_part
= (SCM_COMPLEX_REAL (x
)
4329 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4330 scm_remember_upto_here_1 (x
);
4331 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4333 else if (SCM_REALP (y
))
4334 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4335 SCM_COMPLEX_IMAG (x
));
4336 else if (SCM_COMPLEXP (y
))
4337 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4338 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4339 else if (SCM_FRACTIONP (y
))
4340 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4341 SCM_COMPLEX_IMAG (x
));
4343 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4345 else if (SCM_FRACTIONP (x
))
4347 if (SCM_I_INUMP (y
))
4348 /* a/b - c = (a - cb) / b */
4349 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4350 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4351 SCM_FRACTION_DENOMINATOR (x
));
4352 else if (SCM_BIGP (y
))
4353 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4354 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4355 SCM_FRACTION_DENOMINATOR (x
));
4356 else if (SCM_REALP (y
))
4357 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4358 else if (SCM_COMPLEXP (y
))
4359 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4360 -SCM_COMPLEX_IMAG (y
));
4361 else if (SCM_FRACTIONP (y
))
4362 /* a/b - c/d = (ad - bc) / bd */
4363 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4364 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4365 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4367 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4370 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4375 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4377 "Return @math{@var{x}-1}.")
4378 #define FUNC_NAME s_scm_oneminus
4380 return scm_difference (x
, SCM_I_MAKINUM (1));
4385 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4386 /* "Return the product of all arguments. If called without arguments,\n"
4390 scm_product (SCM x
, SCM y
)
4392 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4395 return SCM_I_MAKINUM (1L);
4396 else if (SCM_NUMBERP (x
))
4399 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4402 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4407 xx
= SCM_I_INUM (x
);
4411 case 0: return x
; break;
4412 case 1: return y
; break;
4415 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4417 long yy
= SCM_I_INUM (y
);
4419 SCM k
= SCM_I_MAKINUM (kk
);
4420 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4424 SCM result
= scm_i_long2big (xx
);
4425 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4426 return scm_i_normbig (result
);
4429 else if (SCM_BIGP (y
))
4431 SCM result
= scm_i_mkbig ();
4432 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4433 scm_remember_upto_here_1 (y
);
4436 else if (SCM_REALP (y
))
4437 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4438 else if (SCM_COMPLEXP (y
))
4439 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4440 xx
* SCM_COMPLEX_IMAG (y
));
4441 else if (SCM_FRACTIONP (y
))
4442 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4443 SCM_FRACTION_DENOMINATOR (y
));
4445 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4447 else if (SCM_BIGP (x
))
4449 if (SCM_I_INUMP (y
))
4454 else if (SCM_BIGP (y
))
4456 SCM result
= scm_i_mkbig ();
4457 mpz_mul (SCM_I_BIG_MPZ (result
),
4460 scm_remember_upto_here_2 (x
, y
);
4463 else if (SCM_REALP (y
))
4465 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4466 scm_remember_upto_here_1 (x
);
4467 return scm_from_double (result
);
4469 else if (SCM_COMPLEXP (y
))
4471 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4472 scm_remember_upto_here_1 (x
);
4473 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4474 z
* SCM_COMPLEX_IMAG (y
));
4476 else if (SCM_FRACTIONP (y
))
4477 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4478 SCM_FRACTION_DENOMINATOR (y
));
4480 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4482 else if (SCM_REALP (x
))
4484 if (SCM_I_INUMP (y
))
4486 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4487 if (scm_is_eq (y
, SCM_INUM0
))
4489 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4491 else if (SCM_BIGP (y
))
4493 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4494 scm_remember_upto_here_1 (y
);
4495 return scm_from_double (result
);
4497 else if (SCM_REALP (y
))
4498 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4499 else if (SCM_COMPLEXP (y
))
4500 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4501 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4502 else if (SCM_FRACTIONP (y
))
4503 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4505 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4507 else if (SCM_COMPLEXP (x
))
4509 if (SCM_I_INUMP (y
))
4511 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4512 if (scm_is_eq (y
, SCM_INUM0
))
4514 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4515 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4517 else if (SCM_BIGP (y
))
4519 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4520 scm_remember_upto_here_1 (y
);
4521 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4522 z
* SCM_COMPLEX_IMAG (x
));
4524 else if (SCM_REALP (y
))
4525 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4526 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4527 else if (SCM_COMPLEXP (y
))
4529 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4530 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4531 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4532 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4534 else if (SCM_FRACTIONP (y
))
4536 double yy
= scm_i_fraction2double (y
);
4537 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4538 yy
* SCM_COMPLEX_IMAG (x
));
4541 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4543 else if (SCM_FRACTIONP (x
))
4545 if (SCM_I_INUMP (y
))
4546 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4547 SCM_FRACTION_DENOMINATOR (x
));
4548 else if (SCM_BIGP (y
))
4549 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4550 SCM_FRACTION_DENOMINATOR (x
));
4551 else if (SCM_REALP (y
))
4552 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4553 else if (SCM_COMPLEXP (y
))
4555 double xx
= scm_i_fraction2double (x
);
4556 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4557 xx
* SCM_COMPLEX_IMAG (y
));
4559 else if (SCM_FRACTIONP (y
))
4560 /* a/b * c/d = ac / bd */
4561 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4562 SCM_FRACTION_NUMERATOR (y
)),
4563 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4564 SCM_FRACTION_DENOMINATOR (y
)));
4566 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4569 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4572 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4573 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4574 #define ALLOW_DIVIDE_BY_ZERO
4575 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4578 /* The code below for complex division is adapted from the GNU
4579 libstdc++, which adapted it from f2c's libF77, and is subject to
4582 /****************************************************************
4583 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4585 Permission to use, copy, modify, and distribute this software
4586 and its documentation for any purpose and without fee is hereby
4587 granted, provided that the above copyright notice appear in all
4588 copies and that both that the copyright notice and this
4589 permission notice and warranty disclaimer appear in supporting
4590 documentation, and that the names of AT&T Bell Laboratories or
4591 Bellcore or any of their entities not be used in advertising or
4592 publicity pertaining to distribution of the software without
4593 specific, written prior permission.
4595 AT&T and Bellcore disclaim all warranties with regard to this
4596 software, including all implied warranties of merchantability
4597 and fitness. In no event shall AT&T or Bellcore be liable for
4598 any special, indirect or consequential damages or any damages
4599 whatsoever resulting from loss of use, data or profits, whether
4600 in an action of contract, negligence or other tortious action,
4601 arising out of or in connection with the use or performance of
4603 ****************************************************************/
4605 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4606 /* Divide the first argument by the product of the remaining
4607 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4609 #define FUNC_NAME s_divide
4611 scm_i_divide (SCM x
, SCM y
, int inexact
)
4615 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4618 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4619 else if (SCM_I_INUMP (x
))
4621 long xx
= SCM_I_INUM (x
);
4622 if (xx
== 1 || xx
== -1)
4624 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4626 scm_num_overflow (s_divide
);
4631 return scm_from_double (1.0 / (double) xx
);
4632 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4635 else if (SCM_BIGP (x
))
4638 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4639 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4641 else if (SCM_REALP (x
))
4643 double xx
= SCM_REAL_VALUE (x
);
4644 #ifndef ALLOW_DIVIDE_BY_ZERO
4646 scm_num_overflow (s_divide
);
4649 return scm_from_double (1.0 / xx
);
4651 else if (SCM_COMPLEXP (x
))
4653 double r
= SCM_COMPLEX_REAL (x
);
4654 double i
= SCM_COMPLEX_IMAG (x
);
4655 if (fabs(r
) <= fabs(i
))
4658 double d
= i
* (1.0 + t
* t
);
4659 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4664 double d
= r
* (1.0 + t
* t
);
4665 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4668 else if (SCM_FRACTIONP (x
))
4669 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4670 SCM_FRACTION_NUMERATOR (x
));
4672 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4675 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4677 long xx
= SCM_I_INUM (x
);
4678 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4680 long yy
= SCM_I_INUM (y
);
4683 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4684 scm_num_overflow (s_divide
);
4686 return scm_from_double ((double) xx
/ (double) yy
);
4689 else if (xx
% yy
!= 0)
4692 return scm_from_double ((double) xx
/ (double) yy
);
4693 else return scm_i_make_ratio (x
, y
);
4698 if (SCM_FIXABLE (z
))
4699 return SCM_I_MAKINUM (z
);
4701 return scm_i_long2big (z
);
4704 else if (SCM_BIGP (y
))
4707 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4708 else return scm_i_make_ratio (x
, y
);
4710 else if (SCM_REALP (y
))
4712 double yy
= SCM_REAL_VALUE (y
);
4713 #ifndef ALLOW_DIVIDE_BY_ZERO
4715 scm_num_overflow (s_divide
);
4718 return scm_from_double ((double) xx
/ yy
);
4720 else if (SCM_COMPLEXP (y
))
4723 complex_div
: /* y _must_ be a complex number */
4725 double r
= SCM_COMPLEX_REAL (y
);
4726 double i
= SCM_COMPLEX_IMAG (y
);
4727 if (fabs(r
) <= fabs(i
))
4730 double d
= i
* (1.0 + t
* t
);
4731 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4736 double d
= r
* (1.0 + t
* t
);
4737 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4741 else if (SCM_FRACTIONP (y
))
4742 /* a / b/c = ac / b */
4743 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4744 SCM_FRACTION_NUMERATOR (y
));
4746 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4748 else if (SCM_BIGP (x
))
4750 if (SCM_I_INUMP (y
))
4752 long int yy
= SCM_I_INUM (y
);
4755 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4756 scm_num_overflow (s_divide
);
4758 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4759 scm_remember_upto_here_1 (x
);
4760 return (sgn
== 0) ? scm_nan () : scm_inf ();
4767 /* FIXME: HMM, what are the relative performance issues here?
4768 We need to test. Is it faster on average to test
4769 divisible_p, then perform whichever operation, or is it
4770 faster to perform the integer div opportunistically and
4771 switch to real if there's a remainder? For now we take the
4772 middle ground: test, then if divisible, use the faster div
4775 long abs_yy
= yy
< 0 ? -yy
: yy
;
4776 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4780 SCM result
= scm_i_mkbig ();
4781 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4782 scm_remember_upto_here_1 (x
);
4784 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4785 return scm_i_normbig (result
);
4790 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4791 else return scm_i_make_ratio (x
, y
);
4795 else if (SCM_BIGP (y
))
4797 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4800 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4801 scm_num_overflow (s_divide
);
4803 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4804 scm_remember_upto_here_1 (x
);
4805 return (sgn
== 0) ? scm_nan () : scm_inf ();
4813 /* It's easily possible for the ratio x/y to fit a double
4814 but one or both x and y be too big to fit a double,
4815 hence the use of mpq_get_d rather than converting and
4818 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4819 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4820 return scm_from_double (mpq_get_d (q
));
4824 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4828 SCM result
= scm_i_mkbig ();
4829 mpz_divexact (SCM_I_BIG_MPZ (result
),
4832 scm_remember_upto_here_2 (x
, y
);
4833 return scm_i_normbig (result
);
4836 return scm_i_make_ratio (x
, y
);
4840 else if (SCM_REALP (y
))
4842 double yy
= SCM_REAL_VALUE (y
);
4843 #ifndef ALLOW_DIVIDE_BY_ZERO
4845 scm_num_overflow (s_divide
);
4848 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4850 else if (SCM_COMPLEXP (y
))
4852 a
= scm_i_big2dbl (x
);
4855 else if (SCM_FRACTIONP (y
))
4856 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4857 SCM_FRACTION_NUMERATOR (y
));
4859 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4861 else if (SCM_REALP (x
))
4863 double rx
= SCM_REAL_VALUE (x
);
4864 if (SCM_I_INUMP (y
))
4866 long int yy
= SCM_I_INUM (y
);
4867 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4869 scm_num_overflow (s_divide
);
4872 return scm_from_double (rx
/ (double) yy
);
4874 else if (SCM_BIGP (y
))
4876 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4877 scm_remember_upto_here_1 (y
);
4878 return scm_from_double (rx
/ dby
);
4880 else if (SCM_REALP (y
))
4882 double yy
= SCM_REAL_VALUE (y
);
4883 #ifndef ALLOW_DIVIDE_BY_ZERO
4885 scm_num_overflow (s_divide
);
4888 return scm_from_double (rx
/ yy
);
4890 else if (SCM_COMPLEXP (y
))
4895 else if (SCM_FRACTIONP (y
))
4896 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4898 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4900 else if (SCM_COMPLEXP (x
))
4902 double rx
= SCM_COMPLEX_REAL (x
);
4903 double ix
= SCM_COMPLEX_IMAG (x
);
4904 if (SCM_I_INUMP (y
))
4906 long int yy
= SCM_I_INUM (y
);
4907 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4909 scm_num_overflow (s_divide
);
4914 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4917 else if (SCM_BIGP (y
))
4919 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4920 scm_remember_upto_here_1 (y
);
4921 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4923 else if (SCM_REALP (y
))
4925 double yy
= SCM_REAL_VALUE (y
);
4926 #ifndef ALLOW_DIVIDE_BY_ZERO
4928 scm_num_overflow (s_divide
);
4931 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4933 else if (SCM_COMPLEXP (y
))
4935 double ry
= SCM_COMPLEX_REAL (y
);
4936 double iy
= SCM_COMPLEX_IMAG (y
);
4937 if (fabs(ry
) <= fabs(iy
))
4940 double d
= iy
* (1.0 + t
* t
);
4941 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4946 double d
= ry
* (1.0 + t
* t
);
4947 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4950 else if (SCM_FRACTIONP (y
))
4952 double yy
= scm_i_fraction2double (y
);
4953 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4956 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4958 else if (SCM_FRACTIONP (x
))
4960 if (SCM_I_INUMP (y
))
4962 long int yy
= SCM_I_INUM (y
);
4963 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4965 scm_num_overflow (s_divide
);
4968 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4969 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4971 else if (SCM_BIGP (y
))
4973 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4974 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4976 else if (SCM_REALP (y
))
4978 double yy
= SCM_REAL_VALUE (y
);
4979 #ifndef ALLOW_DIVIDE_BY_ZERO
4981 scm_num_overflow (s_divide
);
4984 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4986 else if (SCM_COMPLEXP (y
))
4988 a
= scm_i_fraction2double (x
);
4991 else if (SCM_FRACTIONP (y
))
4992 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4993 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4995 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4998 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5002 scm_divide (SCM x
, SCM y
)
5004 return scm_i_divide (x
, y
, 0);
5007 static SCM
scm_divide2real (SCM x
, SCM y
)
5009 return scm_i_divide (x
, y
, 1);
5015 scm_asinh (double x
)
5020 #define asinh scm_asinh
5021 return log (x
+ sqrt (x
* x
+ 1));
5024 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5025 /* "Return the inverse hyperbolic sine of @var{x}."
5030 scm_acosh (double x
)
5035 #define acosh scm_acosh
5036 return log (x
+ sqrt (x
* x
- 1));
5039 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5040 /* "Return the inverse hyperbolic cosine of @var{x}."
5045 scm_atanh (double x
)
5050 #define atanh scm_atanh
5051 return 0.5 * log ((1 + x
) / (1 - x
));
5054 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5055 /* "Return the inverse hyperbolic tangent of @var{x}."
5060 scm_c_truncate (double x
)
5071 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5072 half-way case (ie. when x is an integer plus 0.5) going upwards.
5073 Then half-way cases are identified and adjusted down if the
5074 round-upwards didn't give the desired even integer.
5076 "plus_half == result" identifies a half-way case. If plus_half, which is
5077 x + 0.5, is an integer then x must be an integer plus 0.5.
5079 An odd "result" value is identified with result/2 != floor(result/2).
5080 This is done with plus_half, since that value is ready for use sooner in
5081 a pipelined cpu, and we're already requiring plus_half == result.
5083 Note however that we need to be careful when x is big and already an
5084 integer. In that case "x+0.5" may round to an adjacent integer, causing
5085 us to return such a value, incorrectly. For instance if the hardware is
5086 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5087 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5088 returned. Or if the hardware is in round-upwards mode, then other bigger
5089 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5090 representable value, 2^128+2^76 (or whatever), again incorrect.
5092 These bad roundings of x+0.5 are avoided by testing at the start whether
5093 x is already an integer. If it is then clearly that's the desired result
5094 already. And if it's not then the exponent must be small enough to allow
5095 an 0.5 to be represented, and hence added without a bad rounding. */
5098 scm_c_round (double x
)
5100 double plus_half
, result
;
5105 plus_half
= x
+ 0.5;
5106 result
= floor (plus_half
);
5107 /* Adjust so that the rounding is towards even. */
5108 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5113 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5115 "Round the number @var{x} towards zero.")
5116 #define FUNC_NAME s_scm_truncate_number
5118 if (scm_is_false (scm_negative_p (x
)))
5119 return scm_floor (x
);
5121 return scm_ceiling (x
);
5125 static SCM exactly_one_half
;
5127 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5129 "Round the number @var{x} towards the nearest integer. "
5130 "When it is exactly halfway between two integers, "
5131 "round towards the even one.")
5132 #define FUNC_NAME s_scm_round_number
5134 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5136 else if (SCM_REALP (x
))
5137 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5140 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5141 single quotient+remainder division then examining to see which way
5142 the rounding should go. */
5143 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5144 SCM result
= scm_floor (plus_half
);
5145 /* Adjust so that the rounding is towards even. */
5146 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5147 && scm_is_true (scm_odd_p (result
)))
5148 return scm_difference (result
, SCM_I_MAKINUM (1));
5155 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5157 "Round the number @var{x} towards minus infinity.")
5158 #define FUNC_NAME s_scm_floor
5160 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5162 else if (SCM_REALP (x
))
5163 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5164 else if (SCM_FRACTIONP (x
))
5166 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5167 SCM_FRACTION_DENOMINATOR (x
));
5168 if (scm_is_false (scm_negative_p (x
)))
5170 /* For positive x, rounding towards zero is correct. */
5175 /* For negative x, we need to return q-1 unless x is an
5176 integer. But fractions are never integer, per our
5178 return scm_difference (q
, SCM_I_MAKINUM (1));
5182 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5186 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5188 "Round the number @var{x} towards infinity.")
5189 #define FUNC_NAME s_scm_ceiling
5191 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5193 else if (SCM_REALP (x
))
5194 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5195 else if (SCM_FRACTIONP (x
))
5197 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5198 SCM_FRACTION_DENOMINATOR (x
));
5199 if (scm_is_false (scm_positive_p (x
)))
5201 /* For negative x, rounding towards zero is correct. */
5206 /* For positive x, we need to return q+1 unless x is an
5207 integer. But fractions are never integer, per our
5209 return scm_sum (q
, SCM_I_MAKINUM (1));
5213 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5217 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5218 /* "Return the square root of the real number @var{x}."
5220 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5221 /* "Return the absolute value of the real number @var{x}."
5223 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5224 /* "Return the @var{x}th power of e."
5226 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5227 /* "Return the natural logarithm of the real number @var{x}."
5229 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5230 /* "Return the sine of the real number @var{x}."
5232 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5233 /* "Return the cosine of the real number @var{x}."
5235 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5236 /* "Return the tangent of the real number @var{x}."
5238 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5239 /* "Return the arc sine of the real number @var{x}."
5241 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5242 /* "Return the arc cosine of the real number @var{x}."
5244 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5245 /* "Return the arc tangent of the real number @var{x}."
5247 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5248 /* "Return the hyperbolic sine of the real number @var{x}."
5250 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5251 /* "Return the hyperbolic cosine of the real number @var{x}."
5253 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5254 /* "Return the hyperbolic tangent of the real number @var{x}."
5262 static void scm_two_doubles (SCM x
,
5264 const char *sstring
,
5268 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5270 if (SCM_I_INUMP (x
))
5271 xy
->x
= SCM_I_INUM (x
);
5272 else if (SCM_BIGP (x
))
5273 xy
->x
= scm_i_big2dbl (x
);
5274 else if (SCM_REALP (x
))
5275 xy
->x
= SCM_REAL_VALUE (x
);
5276 else if (SCM_FRACTIONP (x
))
5277 xy
->x
= scm_i_fraction2double (x
);
5279 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5281 if (SCM_I_INUMP (y
))
5282 xy
->y
= SCM_I_INUM (y
);
5283 else if (SCM_BIGP (y
))
5284 xy
->y
= scm_i_big2dbl (y
);
5285 else if (SCM_REALP (y
))
5286 xy
->y
= SCM_REAL_VALUE (y
);
5287 else if (SCM_FRACTIONP (y
))
5288 xy
->y
= scm_i_fraction2double (y
);
5290 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5294 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5296 "Return @var{x} raised to the power of @var{y}. This\n"
5297 "procedure does not accept complex arguments.")
5298 #define FUNC_NAME s_scm_sys_expt
5301 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5302 return scm_from_double (pow (xy
.x
, xy
.y
));
5307 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5309 "Return the arc tangent of the two arguments @var{x} and\n"
5310 "@var{y}. This is similar to calculating the arc tangent of\n"
5311 "@var{x} / @var{y}, except that the signs of both arguments\n"
5312 "are used to determine the quadrant of the result. This\n"
5313 "procedure does not accept complex arguments.")
5314 #define FUNC_NAME s_scm_sys_atan2
5317 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5318 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5323 scm_c_make_rectangular (double re
, double im
)
5326 return scm_from_double (re
);
5330 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5332 SCM_COMPLEX_REAL (z
) = re
;
5333 SCM_COMPLEX_IMAG (z
) = im
;
5338 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5339 (SCM real_part
, SCM imaginary_part
),
5340 "Return a complex number constructed of the given @var{real-part} "
5341 "and @var{imaginary-part} parts.")
5342 #define FUNC_NAME s_scm_make_rectangular
5345 scm_two_doubles (real_part
, imaginary_part
, FUNC_NAME
, &xy
);
5346 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5351 scm_c_make_polar (double mag
, double ang
)
5355 sincos (ang
, &s
, &c
);
5360 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5363 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5365 "Return the complex number @var{x} * e^(i * @var{y}).")
5366 #define FUNC_NAME s_scm_make_polar
5369 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5370 return scm_c_make_polar (xy
.x
, xy
.y
);
5375 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5376 /* "Return the real part of the number @var{z}."
5379 scm_real_part (SCM z
)
5381 if (SCM_I_INUMP (z
))
5383 else if (SCM_BIGP (z
))
5385 else if (SCM_REALP (z
))
5387 else if (SCM_COMPLEXP (z
))
5388 return scm_from_double (SCM_COMPLEX_REAL (z
));
5389 else if (SCM_FRACTIONP (z
))
5392 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5396 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5397 /* "Return the imaginary part of the number @var{z}."
5400 scm_imag_part (SCM z
)
5402 if (SCM_I_INUMP (z
))
5404 else if (SCM_BIGP (z
))
5406 else if (SCM_REALP (z
))
5408 else if (SCM_COMPLEXP (z
))
5409 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5410 else if (SCM_FRACTIONP (z
))
5413 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5416 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5417 /* "Return the numerator of the number @var{z}."
5420 scm_numerator (SCM z
)
5422 if (SCM_I_INUMP (z
))
5424 else if (SCM_BIGP (z
))
5426 else if (SCM_FRACTIONP (z
))
5427 return SCM_FRACTION_NUMERATOR (z
);
5428 else if (SCM_REALP (z
))
5429 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5431 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5435 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5436 /* "Return the denominator of the number @var{z}."
5439 scm_denominator (SCM z
)
5441 if (SCM_I_INUMP (z
))
5442 return SCM_I_MAKINUM (1);
5443 else if (SCM_BIGP (z
))
5444 return SCM_I_MAKINUM (1);
5445 else if (SCM_FRACTIONP (z
))
5446 return SCM_FRACTION_DENOMINATOR (z
);
5447 else if (SCM_REALP (z
))
5448 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5450 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5453 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5454 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5455 * "@code{abs} for real arguments, but also allows complex numbers."
5458 scm_magnitude (SCM z
)
5460 if (SCM_I_INUMP (z
))
5462 long int zz
= SCM_I_INUM (z
);
5465 else if (SCM_POSFIXABLE (-zz
))
5466 return SCM_I_MAKINUM (-zz
);
5468 return scm_i_long2big (-zz
);
5470 else if (SCM_BIGP (z
))
5472 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5473 scm_remember_upto_here_1 (z
);
5475 return scm_i_clonebig (z
, 0);
5479 else if (SCM_REALP (z
))
5480 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5481 else if (SCM_COMPLEXP (z
))
5482 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5483 else if (SCM_FRACTIONP (z
))
5485 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5487 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5488 SCM_FRACTION_DENOMINATOR (z
));
5491 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5495 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5496 /* "Return the angle of the complex number @var{z}."
5501 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5502 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5503 But if atan2 follows the floating point rounding mode, then the value
5504 is not a constant. Maybe it'd be close enough though. */
5505 if (SCM_I_INUMP (z
))
5507 if (SCM_I_INUM (z
) >= 0)
5510 return scm_from_double (atan2 (0.0, -1.0));
5512 else if (SCM_BIGP (z
))
5514 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5515 scm_remember_upto_here_1 (z
);
5517 return scm_from_double (atan2 (0.0, -1.0));
5521 else if (SCM_REALP (z
))
5523 if (SCM_REAL_VALUE (z
) >= 0)
5526 return scm_from_double (atan2 (0.0, -1.0));
5528 else if (SCM_COMPLEXP (z
))
5529 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5530 else if (SCM_FRACTIONP (z
))
5532 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5534 else return scm_from_double (atan2 (0.0, -1.0));
5537 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5541 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5542 /* Convert the number @var{x} to its inexact representation.\n"
5545 scm_exact_to_inexact (SCM z
)
5547 if (SCM_I_INUMP (z
))
5548 return scm_from_double ((double) SCM_I_INUM (z
));
5549 else if (SCM_BIGP (z
))
5550 return scm_from_double (scm_i_big2dbl (z
));
5551 else if (SCM_FRACTIONP (z
))
5552 return scm_from_double (scm_i_fraction2double (z
));
5553 else if (SCM_INEXACTP (z
))
5556 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5560 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5562 "Return an exact number that is numerically closest to @var{z}.")
5563 #define FUNC_NAME s_scm_inexact_to_exact
5565 if (SCM_I_INUMP (z
))
5567 else if (SCM_BIGP (z
))
5569 else if (SCM_REALP (z
))
5571 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5572 SCM_OUT_OF_RANGE (1, z
);
5579 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5580 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5581 scm_i_mpz2num (mpq_denref (frac
)));
5583 /* When scm_i_make_ratio throws, we leak the memory allocated
5590 else if (SCM_FRACTIONP (z
))
5593 SCM_WRONG_TYPE_ARG (1, z
);
5597 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5599 "Return an exact number that is within @var{err} of @var{x}.")
5600 #define FUNC_NAME s_scm_rationalize
5602 if (SCM_I_INUMP (x
))
5604 else if (SCM_BIGP (x
))
5606 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5608 /* Use continued fractions to find closest ratio. All
5609 arithmetic is done with exact numbers.
5612 SCM ex
= scm_inexact_to_exact (x
);
5613 SCM int_part
= scm_floor (ex
);
5614 SCM tt
= SCM_I_MAKINUM (1);
5615 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5616 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5620 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5623 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5624 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5626 /* We stop after a million iterations just to be absolutely sure
5627 that we don't go into an infinite loop. The process normally
5628 converges after less than a dozen iterations.
5631 err
= scm_abs (err
);
5632 while (++i
< 1000000)
5634 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5635 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5636 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5638 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5639 err
))) /* abs(x-a/b) <= err */
5641 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5642 if (scm_is_false (scm_exact_p (x
))
5643 || scm_is_false (scm_exact_p (err
)))
5644 return scm_exact_to_inexact (res
);
5648 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5650 tt
= scm_floor (rx
); /* tt = floor (rx) */
5656 scm_num_overflow (s_scm_rationalize
);
5659 SCM_WRONG_TYPE_ARG (1, x
);
5663 /* conversion functions */
5666 scm_is_integer (SCM val
)
5668 return scm_is_true (scm_integer_p (val
));
5672 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5674 if (SCM_I_INUMP (val
))
5676 scm_t_signed_bits n
= SCM_I_INUM (val
);
5677 return n
>= min
&& n
<= max
;
5679 else if (SCM_BIGP (val
))
5681 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5683 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5685 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5687 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5688 return n
>= min
&& n
<= max
;
5698 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5699 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5702 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5703 SCM_I_BIG_MPZ (val
));
5705 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5717 return n
>= min
&& n
<= max
;
5725 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5727 if (SCM_I_INUMP (val
))
5729 scm_t_signed_bits n
= SCM_I_INUM (val
);
5730 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5732 else if (SCM_BIGP (val
))
5734 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5736 else if (max
<= ULONG_MAX
)
5738 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5740 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5741 return n
>= min
&& n
<= max
;
5751 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5754 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5755 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5758 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5759 SCM_I_BIG_MPZ (val
));
5761 return n
>= min
&& n
<= max
;
5769 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5771 scm_error (scm_out_of_range_key
,
5773 "Value out of range ~S to ~S: ~S",
5774 scm_list_3 (min
, max
, bad_val
),
5775 scm_list_1 (bad_val
));
5778 #define TYPE scm_t_intmax
5779 #define TYPE_MIN min
5780 #define TYPE_MAX max
5781 #define SIZEOF_TYPE 0
5782 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5783 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5784 #include "libguile/conv-integer.i.c"
5786 #define TYPE scm_t_uintmax
5787 #define TYPE_MIN min
5788 #define TYPE_MAX max
5789 #define SIZEOF_TYPE 0
5790 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5791 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5792 #include "libguile/conv-uinteger.i.c"
5794 #define TYPE scm_t_int8
5795 #define TYPE_MIN SCM_T_INT8_MIN
5796 #define TYPE_MAX SCM_T_INT8_MAX
5797 #define SIZEOF_TYPE 1
5798 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5799 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5800 #include "libguile/conv-integer.i.c"
5802 #define TYPE scm_t_uint8
5804 #define TYPE_MAX SCM_T_UINT8_MAX
5805 #define SIZEOF_TYPE 1
5806 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5807 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5808 #include "libguile/conv-uinteger.i.c"
5810 #define TYPE scm_t_int16
5811 #define TYPE_MIN SCM_T_INT16_MIN
5812 #define TYPE_MAX SCM_T_INT16_MAX
5813 #define SIZEOF_TYPE 2
5814 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5815 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5816 #include "libguile/conv-integer.i.c"
5818 #define TYPE scm_t_uint16
5820 #define TYPE_MAX SCM_T_UINT16_MAX
5821 #define SIZEOF_TYPE 2
5822 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5823 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5824 #include "libguile/conv-uinteger.i.c"
5826 #define TYPE scm_t_int32
5827 #define TYPE_MIN SCM_T_INT32_MIN
5828 #define TYPE_MAX SCM_T_INT32_MAX
5829 #define SIZEOF_TYPE 4
5830 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5831 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5832 #include "libguile/conv-integer.i.c"
5834 #define TYPE scm_t_uint32
5836 #define TYPE_MAX SCM_T_UINT32_MAX
5837 #define SIZEOF_TYPE 4
5838 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5839 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5840 #include "libguile/conv-uinteger.i.c"
5842 #if SCM_HAVE_T_INT64
5844 #define TYPE scm_t_int64
5845 #define TYPE_MIN SCM_T_INT64_MIN
5846 #define TYPE_MAX SCM_T_INT64_MAX
5847 #define SIZEOF_TYPE 8
5848 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5849 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5850 #include "libguile/conv-integer.i.c"
5852 #define TYPE scm_t_uint64
5854 #define TYPE_MAX SCM_T_UINT64_MAX
5855 #define SIZEOF_TYPE 8
5856 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5857 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5858 #include "libguile/conv-uinteger.i.c"
5863 scm_to_mpz (SCM val
, mpz_t rop
)
5865 if (SCM_I_INUMP (val
))
5866 mpz_set_si (rop
, SCM_I_INUM (val
));
5867 else if (SCM_BIGP (val
))
5868 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5870 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5874 scm_from_mpz (mpz_t val
)
5876 return scm_i_mpz2num (val
);
5880 scm_is_real (SCM val
)
5882 return scm_is_true (scm_real_p (val
));
5886 scm_is_rational (SCM val
)
5888 return scm_is_true (scm_rational_p (val
));
5892 scm_to_double (SCM val
)
5894 if (SCM_I_INUMP (val
))
5895 return SCM_I_INUM (val
);
5896 else if (SCM_BIGP (val
))
5897 return scm_i_big2dbl (val
);
5898 else if (SCM_FRACTIONP (val
))
5899 return scm_i_fraction2double (val
);
5900 else if (SCM_REALP (val
))
5901 return SCM_REAL_VALUE (val
);
5903 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5907 scm_from_double (double val
)
5909 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5910 SCM_REAL_VALUE (z
) = val
;
5914 #if SCM_ENABLE_DISCOURAGED == 1
5917 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5921 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5925 scm_out_of_range (NULL
, num
);
5928 return scm_to_double (num
);
5932 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5936 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5940 scm_out_of_range (NULL
, num
);
5943 return scm_to_double (num
);
5949 scm_is_complex (SCM val
)
5951 return scm_is_true (scm_complex_p (val
));
5955 scm_c_real_part (SCM z
)
5957 if (SCM_COMPLEXP (z
))
5958 return SCM_COMPLEX_REAL (z
);
5961 /* Use the scm_real_part to get proper error checking and
5964 return scm_to_double (scm_real_part (z
));
5969 scm_c_imag_part (SCM z
)
5971 if (SCM_COMPLEXP (z
))
5972 return SCM_COMPLEX_IMAG (z
);
5975 /* Use the scm_imag_part to get proper error checking and
5976 dispatching. The result will almost always be 0.0, but not
5979 return scm_to_double (scm_imag_part (z
));
5984 scm_c_magnitude (SCM z
)
5986 return scm_to_double (scm_magnitude (z
));
5992 return scm_to_double (scm_angle (z
));
5996 scm_is_number (SCM z
)
5998 return scm_is_true (scm_number_p (z
));
6002 /* In the following functions we dispatch to the real-arg funcs like log()
6003 when we know the arg is real, instead of just handing everything to
6004 clog() for instance. This is in case clog() doesn't optimize for a
6005 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6006 well use it to go straight to the applicable C func. */
6008 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6010 "Return the natural logarithm of @var{z}.")
6011 #define FUNC_NAME s_scm_log
6013 if (SCM_COMPLEXP (z
))
6015 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6016 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6018 double re
= SCM_COMPLEX_REAL (z
);
6019 double im
= SCM_COMPLEX_IMAG (z
);
6020 return scm_c_make_rectangular (log (hypot (re
, im
)),
6026 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6027 although the value itself overflows. */
6028 double re
= scm_to_double (z
);
6029 double l
= log (fabs (re
));
6031 return scm_from_double (l
);
6033 return scm_c_make_rectangular (l
, M_PI
);
6039 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6041 "Return the base 10 logarithm of @var{z}.")
6042 #define FUNC_NAME s_scm_log10
6044 if (SCM_COMPLEXP (z
))
6046 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6047 clog() and a multiply by M_LOG10E, rather than the fallback
6048 log10+hypot+atan2.) */
6049 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6050 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6052 double re
= SCM_COMPLEX_REAL (z
);
6053 double im
= SCM_COMPLEX_IMAG (z
);
6054 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6055 M_LOG10E
* atan2 (im
, re
));
6060 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6061 although the value itself overflows. */
6062 double re
= scm_to_double (z
);
6063 double l
= log10 (fabs (re
));
6065 return scm_from_double (l
);
6067 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6073 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6075 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6076 "base of natural logarithms (2.71828@dots{}).")
6077 #define FUNC_NAME s_scm_exp
6079 if (SCM_COMPLEXP (z
))
6081 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6082 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6084 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6085 SCM_COMPLEX_IMAG (z
));
6090 /* When z is a negative bignum the conversion to double overflows,
6091 giving -infinity, but that's ok, the exp is still 0.0. */
6092 return scm_from_double (exp (scm_to_double (z
)));
6098 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6100 "Return the square root of @var{z}. Of the two possible roots\n"
6101 "(positive and negative), the one with the a positive real part\n"
6102 "is returned, or if that's zero then a positive imaginary part.\n"
6106 "(sqrt 9.0) @result{} 3.0\n"
6107 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6108 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6109 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6111 #define FUNC_NAME s_scm_sqrt
6113 if (SCM_COMPLEXP (x
))
6115 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6116 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6118 double re
= SCM_COMPLEX_REAL (x
);
6119 double im
= SCM_COMPLEX_IMAG (x
);
6120 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6121 0.5 * atan2 (im
, re
));
6126 double xx
= scm_to_double (x
);
6128 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6130 return scm_from_double (sqrt (xx
));
6142 mpz_init_set_si (z_negative_one
, -1);
6144 /* It may be possible to tune the performance of some algorithms by using
6145 * the following constants to avoid the creation of bignums. Please, before
6146 * using these values, remember the two rules of program optimization:
6147 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6148 scm_c_define ("most-positive-fixnum",
6149 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6150 scm_c_define ("most-negative-fixnum",
6151 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6153 scm_add_feature ("complex");
6154 scm_add_feature ("inexact");
6155 scm_flo0
= scm_from_double (0.0);
6157 /* determine floating point precision */
6158 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6160 init_dblprec(&scm_dblprec
[i
-2],i
);
6161 init_fx_radix(fx_per_radix
[i
-2],i
);
6164 /* hard code precision for base 10 if the preprocessor tells us to... */
6165 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6168 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6169 SCM_I_MAKINUM (2)));
6170 #include "libguile/numbers.x"