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.
55 #include "libguile/_scm.h"
56 #include "libguile/feature.h"
57 #include "libguile/ports.h"
58 #include "libguile/root.h"
59 #include "libguile/smob.h"
60 #include "libguile/strings.h"
62 #include "libguile/validate.h"
63 #include "libguile/numbers.h"
64 #include "libguile/deprecation.h"
66 #include "libguile/eq.h"
68 #include "libguile/discouraged.h"
70 /* values per glibc, if not already defined */
72 #define M_LOG10E 0.43429448190325182765
75 #define M_PI 3.14159265358979323846
81 Wonder if this might be faster for some of our code? A switch on
82 the numtag would jump directly to the right case, and the
83 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
85 #define SCM_I_NUMTAG_NOTNUM 0
86 #define SCM_I_NUMTAG_INUM 1
87 #define SCM_I_NUMTAG_BIG scm_tc16_big
88 #define SCM_I_NUMTAG_REAL scm_tc16_real
89 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
90 #define SCM_I_NUMTAG(x) \
91 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
92 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
93 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
94 : SCM_I_NUMTAG_NOTNUM)))
96 /* the macro above will not work as is with fractions */
99 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
101 /* FLOBUFLEN is the maximum number of characters neccessary for the
102 * printed or scm_string representation of an inexact number.
104 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
107 #if ! defined (HAVE_ISNAN)
112 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
115 #if ! defined (HAVE_ISINF)
120 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
127 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
128 an explicit check. In some future gmp (don't know what version number),
129 mpz_cmp_d is supposed to do this itself. */
131 #define xmpz_cmp_d(z, d) \
132 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
134 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
137 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
138 isinf. It does have finite and isnan though, hence the use of those.
139 fpclass would be a possibility on that system too. */
143 #if defined (HAVE_ISINF)
145 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
146 return (! (finite (x
) || isnan (x
)));
155 #if defined (HAVE_ISNAN)
162 #if defined (GUILE_I)
163 #if HAVE_COMPLEX_DOUBLE
165 /* For an SCM object Z which is a complex number (ie. satisfies
166 SCM_COMPLEXP), return its value as a C level "complex double". */
167 #define SCM_COMPLEX_VALUE(z) \
168 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
170 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
172 /* Convert a C "complex double" to an SCM value. */
174 scm_from_complex_double (complex double z
)
176 return scm_c_make_rectangular (creal (z
), cimag (z
));
179 #endif /* HAVE_COMPLEX_DOUBLE */
184 static mpz_t z_negative_one
;
191 /* Return a newly created bignum. */
192 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
193 mpz_init (SCM_I_BIG_MPZ (z
));
198 scm_i_long2big (long x
)
200 /* Return a newly created bignum initialized to X. */
201 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
202 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
207 scm_i_ulong2big (unsigned long x
)
209 /* Return a newly created bignum initialized to X. */
210 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
211 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
216 scm_i_clonebig (SCM src_big
, int same_sign_p
)
218 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
219 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
220 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
222 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
227 scm_i_bigcmp (SCM x
, SCM y
)
229 /* Return neg if x < y, pos if x > y, and 0 if x == y */
230 /* presume we already know x and y are bignums */
231 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
232 scm_remember_upto_here_2 (x
, y
);
237 scm_i_dbl2big (double d
)
239 /* results are only defined if d is an integer */
240 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
241 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
245 /* Convert a integer in double representation to a SCM number. */
248 scm_i_dbl2num (double u
)
250 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
251 powers of 2, so there's no rounding when making "double" values
252 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
253 get rounded on a 64-bit machine, hence the "+1".
255 The use of floor() to force to an integer value ensures we get a
256 "numerically closest" value without depending on how a
257 double->long cast or how mpz_set_d will round. For reference,
258 double->long probably follows the hardware rounding mode,
259 mpz_set_d truncates towards zero. */
261 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
262 representable as a double? */
264 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
265 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
266 return SCM_I_MAKINUM ((long) u
);
268 return scm_i_dbl2big (u
);
271 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
272 with R5RS exact->inexact.
274 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
275 (ie. truncate towards zero), then adjust to get the closest double by
276 examining the next lower bit and adding 1 (to the absolute value) if
279 Bignums exactly half way between representable doubles are rounded to the
280 next higher absolute value (ie. away from zero). This seems like an
281 adequate interpretation of R5RS "numerically closest", and it's easier
282 and faster than a full "nearest-even" style.
284 The bit test must be done on the absolute value of the mpz_t, which means
285 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
286 negatives as twos complement.
288 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
289 following the hardware rounding mode, but applied to the absolute value
290 of the mpz_t operand. This is not what we want so we put the high
291 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
292 mpz_get_d is supposed to always truncate towards zero.
294 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
295 is a slowdown. It'd be faster to pick out the relevant high bits with
296 mpz_getlimbn if we could be bothered coding that, and if the new
297 truncating gmp doesn't come out. */
300 scm_i_big2dbl (SCM b
)
305 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
309 /* Current GMP, eg. 4.1.3, force truncation towards zero */
311 if (bits
> DBL_MANT_DIG
)
313 size_t shift
= bits
- DBL_MANT_DIG
;
314 mpz_init2 (tmp
, DBL_MANT_DIG
);
315 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
316 result
= ldexp (mpz_get_d (tmp
), shift
);
321 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
326 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
329 if (bits
> DBL_MANT_DIG
)
331 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
332 /* test bit number "pos" in absolute value */
333 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
334 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
336 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
340 scm_remember_upto_here_1 (b
);
345 scm_i_normbig (SCM b
)
347 /* convert a big back to a fixnum if it'll fit */
348 /* presume b is a bignum */
349 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
351 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
352 if (SCM_FIXABLE (val
))
353 b
= SCM_I_MAKINUM (val
);
358 static SCM_C_INLINE_KEYWORD SCM
359 scm_i_mpz2num (mpz_t b
)
361 /* convert a mpz number to a SCM number. */
362 if (mpz_fits_slong_p (b
))
364 long val
= mpz_get_si (b
);
365 if (SCM_FIXABLE (val
))
366 return SCM_I_MAKINUM (val
);
370 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
371 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
376 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
377 static SCM
scm_divide2real (SCM x
, SCM y
);
380 scm_i_make_ratio (SCM numerator
, SCM denominator
)
381 #define FUNC_NAME "make-ratio"
383 /* First make sure the arguments are proper.
385 if (SCM_I_INUMP (denominator
))
387 if (scm_is_eq (denominator
, SCM_INUM0
))
388 scm_num_overflow ("make-ratio");
389 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
394 if (!(SCM_BIGP(denominator
)))
395 SCM_WRONG_TYPE_ARG (2, denominator
);
397 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
398 SCM_WRONG_TYPE_ARG (1, numerator
);
400 /* Then flip signs so that the denominator is positive.
402 if (scm_is_true (scm_negative_p (denominator
)))
404 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
405 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
408 /* Now consider for each of the four fixnum/bignum combinations
409 whether the rational number is really an integer.
411 if (SCM_I_INUMP (numerator
))
413 long x
= SCM_I_INUM (numerator
);
414 if (scm_is_eq (numerator
, SCM_INUM0
))
416 if (SCM_I_INUMP (denominator
))
419 y
= SCM_I_INUM (denominator
);
421 return SCM_I_MAKINUM(1);
423 return SCM_I_MAKINUM (x
/ y
);
427 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
428 of that value for the denominator, as a bignum. Apart from
429 that case, abs(bignum) > abs(inum) so inum/bignum is not an
431 if (x
== SCM_MOST_NEGATIVE_FIXNUM
432 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
433 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
434 return SCM_I_MAKINUM(-1);
437 else if (SCM_BIGP (numerator
))
439 if (SCM_I_INUMP (denominator
))
441 long yy
= SCM_I_INUM (denominator
);
442 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
443 return scm_divide (numerator
, denominator
);
447 if (scm_is_eq (numerator
, denominator
))
448 return SCM_I_MAKINUM(1);
449 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
450 SCM_I_BIG_MPZ (denominator
)))
451 return scm_divide(numerator
, denominator
);
455 /* No, it's a proper fraction.
458 SCM divisor
= scm_gcd (numerator
, denominator
);
459 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
461 numerator
= scm_divide (numerator
, divisor
);
462 denominator
= scm_divide (denominator
, divisor
);
465 return scm_double_cell (scm_tc16_fraction
,
466 SCM_UNPACK (numerator
),
467 SCM_UNPACK (denominator
), 0);
473 scm_i_fraction2double (SCM z
)
475 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
476 SCM_FRACTION_DENOMINATOR (z
)));
479 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
481 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
483 #define FUNC_NAME s_scm_exact_p
489 if (SCM_FRACTIONP (x
))
493 SCM_WRONG_TYPE_ARG (1, x
);
498 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
500 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
502 #define FUNC_NAME s_scm_odd_p
506 long val
= SCM_I_INUM (n
);
507 return scm_from_bool ((val
& 1L) != 0);
509 else if (SCM_BIGP (n
))
511 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
512 scm_remember_upto_here_1 (n
);
513 return scm_from_bool (odd_p
);
515 else if (scm_is_true (scm_inf_p (n
)))
517 else if (SCM_REALP (n
))
519 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
525 SCM_WRONG_TYPE_ARG (1, n
);
528 SCM_WRONG_TYPE_ARG (1, n
);
533 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
535 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
537 #define FUNC_NAME s_scm_even_p
541 long val
= SCM_I_INUM (n
);
542 return scm_from_bool ((val
& 1L) == 0);
544 else if (SCM_BIGP (n
))
546 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
547 scm_remember_upto_here_1 (n
);
548 return scm_from_bool (even_p
);
550 else if (scm_is_true (scm_inf_p (n
)))
552 else if (SCM_REALP (n
))
554 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
560 SCM_WRONG_TYPE_ARG (1, n
);
563 SCM_WRONG_TYPE_ARG (1, n
);
567 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
569 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
570 "or @samp{-inf.0}, @code{#f} otherwise.")
571 #define FUNC_NAME s_scm_inf_p
574 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
575 else if (SCM_COMPLEXP (x
))
576 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
577 || xisinf (SCM_COMPLEX_IMAG (x
)));
583 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
585 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
587 #define FUNC_NAME s_scm_nan_p
590 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
591 else if (SCM_COMPLEXP (n
))
592 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
593 || xisnan (SCM_COMPLEX_IMAG (n
)));
599 /* Guile's idea of infinity. */
600 static double guile_Inf
;
602 /* Guile's idea of not a number. */
603 static double guile_NaN
;
606 guile_ieee_init (void)
608 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
610 /* Some version of gcc on some old version of Linux used to crash when
611 trying to make Inf and NaN. */
614 /* C99 INFINITY, when available.
615 FIXME: The standard allows for INFINITY to be something that overflows
616 at compile time. We ought to have a configure test to check for that
617 before trying to use it. (But in practice we believe this is not a
618 problem on any system guile is likely to target.) */
619 guile_Inf
= INFINITY
;
622 extern unsigned int DINFINITY
[2];
623 guile_Inf
= (*((double *) (DINFINITY
)));
630 if (guile_Inf
== tmp
)
638 #if defined (HAVE_ISNAN)
641 /* C99 NAN, when available */
646 extern unsigned int DQNAN
[2];
647 guile_NaN
= (*((double *)(DQNAN
)));
650 guile_NaN
= guile_Inf
/ guile_Inf
;
656 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
659 #define FUNC_NAME s_scm_inf
661 static int initialized
= 0;
667 return scm_from_double (guile_Inf
);
671 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
674 #define FUNC_NAME s_scm_nan
676 static int initialized
= 0;
682 return scm_from_double (guile_NaN
);
687 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
689 "Return the absolute value of @var{x}.")
694 long int xx
= SCM_I_INUM (x
);
697 else if (SCM_POSFIXABLE (-xx
))
698 return SCM_I_MAKINUM (-xx
);
700 return scm_i_long2big (-xx
);
702 else if (SCM_BIGP (x
))
704 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
706 return scm_i_clonebig (x
, 0);
710 else if (SCM_REALP (x
))
712 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
713 double xx
= SCM_REAL_VALUE (x
);
715 return scm_from_double (-xx
);
719 else if (SCM_FRACTIONP (x
))
721 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
723 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
724 SCM_FRACTION_DENOMINATOR (x
));
727 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
732 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
733 /* "Return the quotient of the numbers @var{x} and @var{y}."
736 scm_quotient (SCM x
, SCM y
)
740 long xx
= SCM_I_INUM (x
);
743 long yy
= SCM_I_INUM (y
);
745 scm_num_overflow (s_quotient
);
750 return SCM_I_MAKINUM (z
);
752 return scm_i_long2big (z
);
755 else if (SCM_BIGP (y
))
757 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
758 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
759 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
761 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
762 scm_remember_upto_here_1 (y
);
763 return SCM_I_MAKINUM (-1);
766 return SCM_I_MAKINUM (0);
769 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
771 else if (SCM_BIGP (x
))
775 long yy
= SCM_I_INUM (y
);
777 scm_num_overflow (s_quotient
);
782 SCM result
= scm_i_mkbig ();
785 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
788 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
791 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
792 scm_remember_upto_here_1 (x
);
793 return scm_i_normbig (result
);
796 else if (SCM_BIGP (y
))
798 SCM result
= scm_i_mkbig ();
799 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
802 scm_remember_upto_here_2 (x
, y
);
803 return scm_i_normbig (result
);
806 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
809 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
812 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
813 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
815 * "(remainder 13 4) @result{} 1\n"
816 * "(remainder -13 4) @result{} -1\n"
820 scm_remainder (SCM x
, SCM y
)
826 long yy
= SCM_I_INUM (y
);
828 scm_num_overflow (s_remainder
);
831 long z
= SCM_I_INUM (x
) % yy
;
832 return SCM_I_MAKINUM (z
);
835 else if (SCM_BIGP (y
))
837 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
838 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
839 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
841 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
842 scm_remember_upto_here_1 (y
);
843 return SCM_I_MAKINUM (0);
849 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
851 else if (SCM_BIGP (x
))
855 long yy
= SCM_I_INUM (y
);
857 scm_num_overflow (s_remainder
);
860 SCM result
= scm_i_mkbig ();
863 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
864 scm_remember_upto_here_1 (x
);
865 return scm_i_normbig (result
);
868 else if (SCM_BIGP (y
))
870 SCM result
= scm_i_mkbig ();
871 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
874 scm_remember_upto_here_2 (x
, y
);
875 return scm_i_normbig (result
);
878 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
881 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
885 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
886 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
888 * "(modulo 13 4) @result{} 1\n"
889 * "(modulo -13 4) @result{} 3\n"
893 scm_modulo (SCM x
, SCM y
)
897 long xx
= SCM_I_INUM (x
);
900 long yy
= SCM_I_INUM (y
);
902 scm_num_overflow (s_modulo
);
905 /* C99 specifies that "%" is the remainder corresponding to a
906 quotient rounded towards zero, and that's also traditional
907 for machine division, so z here should be well defined. */
925 return SCM_I_MAKINUM (result
);
928 else if (SCM_BIGP (y
))
930 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
937 SCM pos_y
= scm_i_clonebig (y
, 0);
938 /* do this after the last scm_op */
939 mpz_init_set_si (z_x
, xx
);
940 result
= pos_y
; /* re-use this bignum */
941 mpz_mod (SCM_I_BIG_MPZ (result
),
943 SCM_I_BIG_MPZ (pos_y
));
944 scm_remember_upto_here_1 (pos_y
);
948 result
= scm_i_mkbig ();
949 /* do this after the last scm_op */
950 mpz_init_set_si (z_x
, xx
);
951 mpz_mod (SCM_I_BIG_MPZ (result
),
954 scm_remember_upto_here_1 (y
);
957 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
958 mpz_add (SCM_I_BIG_MPZ (result
),
960 SCM_I_BIG_MPZ (result
));
961 scm_remember_upto_here_1 (y
);
962 /* and do this before the next one */
964 return scm_i_normbig (result
);
968 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
970 else if (SCM_BIGP (x
))
974 long yy
= SCM_I_INUM (y
);
976 scm_num_overflow (s_modulo
);
979 SCM result
= scm_i_mkbig ();
980 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
982 (yy
< 0) ? - yy
: yy
);
983 scm_remember_upto_here_1 (x
);
984 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
985 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
986 SCM_I_BIG_MPZ (result
),
988 return scm_i_normbig (result
);
991 else if (SCM_BIGP (y
))
994 SCM result
= scm_i_mkbig ();
995 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
996 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
997 mpz_mod (SCM_I_BIG_MPZ (result
),
999 SCM_I_BIG_MPZ (pos_y
));
1001 scm_remember_upto_here_1 (x
);
1002 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1003 mpz_add (SCM_I_BIG_MPZ (result
),
1005 SCM_I_BIG_MPZ (result
));
1006 scm_remember_upto_here_2 (y
, pos_y
);
1007 return scm_i_normbig (result
);
1011 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1014 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1017 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
1018 /* "Return the greatest common divisor of all arguments.\n"
1019 * "If called without arguments, 0 is returned."
1022 scm_gcd (SCM x
, SCM y
)
1025 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1027 if (SCM_I_INUMP (x
))
1029 if (SCM_I_INUMP (y
))
1031 long xx
= SCM_I_INUM (x
);
1032 long yy
= SCM_I_INUM (y
);
1033 long u
= xx
< 0 ? -xx
: xx
;
1034 long v
= yy
< 0 ? -yy
: yy
;
1044 /* Determine a common factor 2^k */
1045 while (!(1 & (u
| v
)))
1051 /* Now, any factor 2^n can be eliminated */
1071 return (SCM_POSFIXABLE (result
)
1072 ? SCM_I_MAKINUM (result
)
1073 : scm_i_long2big (result
));
1075 else if (SCM_BIGP (y
))
1081 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1083 else if (SCM_BIGP (x
))
1085 if (SCM_I_INUMP (y
))
1087 unsigned long result
;
1090 yy
= SCM_I_INUM (y
);
1095 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1096 scm_remember_upto_here_1 (x
);
1097 return (SCM_POSFIXABLE (result
)
1098 ? SCM_I_MAKINUM (result
)
1099 : scm_from_ulong (result
));
1101 else if (SCM_BIGP (y
))
1103 SCM result
= scm_i_mkbig ();
1104 mpz_gcd (SCM_I_BIG_MPZ (result
),
1107 scm_remember_upto_here_2 (x
, y
);
1108 return scm_i_normbig (result
);
1111 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1114 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1117 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1118 /* "Return the least common multiple of the arguments.\n"
1119 * "If called without arguments, 1 is returned."
1122 scm_lcm (SCM n1
, SCM n2
)
1124 if (SCM_UNBNDP (n2
))
1126 if (SCM_UNBNDP (n1
))
1127 return SCM_I_MAKINUM (1L);
1128 n2
= SCM_I_MAKINUM (1L);
1131 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1132 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1133 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1134 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1136 if (SCM_I_INUMP (n1
))
1138 if (SCM_I_INUMP (n2
))
1140 SCM d
= scm_gcd (n1
, n2
);
1141 if (scm_is_eq (d
, SCM_INUM0
))
1144 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1148 /* inum n1, big n2 */
1151 SCM result
= scm_i_mkbig ();
1152 long nn1
= SCM_I_INUM (n1
);
1153 if (nn1
== 0) return SCM_INUM0
;
1154 if (nn1
< 0) nn1
= - nn1
;
1155 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1156 scm_remember_upto_here_1 (n2
);
1164 if (SCM_I_INUMP (n2
))
1171 SCM result
= scm_i_mkbig ();
1172 mpz_lcm(SCM_I_BIG_MPZ (result
),
1174 SCM_I_BIG_MPZ (n2
));
1175 scm_remember_upto_here_2(n1
, n2
);
1176 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1182 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1187 + + + x (map digit:logand X Y)
1188 + - + x (map digit:logand X (lognot (+ -1 Y)))
1189 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1190 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1195 + + + (map digit:logior X Y)
1196 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1197 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1198 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1203 + + + (map digit:logxor X Y)
1204 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1205 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1206 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1211 + + (any digit:logand X Y)
1212 + - (any digit:logand X (lognot (+ -1 Y)))
1213 - + (any digit:logand (lognot (+ -1 X)) Y)
1218 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1220 "Return the bitwise AND of the integer arguments.\n\n"
1222 "(logand) @result{} -1\n"
1223 "(logand 7) @result{} 7\n"
1224 "(logand #b111 #b011 #b001) @result{} 1\n"
1226 #define FUNC_NAME s_scm_logand
1230 if (SCM_UNBNDP (n2
))
1232 if (SCM_UNBNDP (n1
))
1233 return SCM_I_MAKINUM (-1);
1234 else if (!SCM_NUMBERP (n1
))
1235 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1236 else if (SCM_NUMBERP (n1
))
1239 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1242 if (SCM_I_INUMP (n1
))
1244 nn1
= SCM_I_INUM (n1
);
1245 if (SCM_I_INUMP (n2
))
1247 long nn2
= SCM_I_INUM (n2
);
1248 return SCM_I_MAKINUM (nn1
& nn2
);
1250 else if SCM_BIGP (n2
)
1256 SCM result_z
= scm_i_mkbig ();
1258 mpz_init_set_si (nn1_z
, nn1
);
1259 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1260 scm_remember_upto_here_1 (n2
);
1262 return scm_i_normbig (result_z
);
1266 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1268 else if (SCM_BIGP (n1
))
1270 if (SCM_I_INUMP (n2
))
1273 nn1
= SCM_I_INUM (n1
);
1276 else if (SCM_BIGP (n2
))
1278 SCM result_z
= scm_i_mkbig ();
1279 mpz_and (SCM_I_BIG_MPZ (result_z
),
1281 SCM_I_BIG_MPZ (n2
));
1282 scm_remember_upto_here_2 (n1
, n2
);
1283 return scm_i_normbig (result_z
);
1286 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1294 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1296 "Return the bitwise OR of the integer arguments.\n\n"
1298 "(logior) @result{} 0\n"
1299 "(logior 7) @result{} 7\n"
1300 "(logior #b000 #b001 #b011) @result{} 3\n"
1302 #define FUNC_NAME s_scm_logior
1306 if (SCM_UNBNDP (n2
))
1308 if (SCM_UNBNDP (n1
))
1310 else if (SCM_NUMBERP (n1
))
1313 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1316 if (SCM_I_INUMP (n1
))
1318 nn1
= SCM_I_INUM (n1
);
1319 if (SCM_I_INUMP (n2
))
1321 long nn2
= SCM_I_INUM (n2
);
1322 return SCM_I_MAKINUM (nn1
| nn2
);
1324 else if (SCM_BIGP (n2
))
1330 SCM result_z
= scm_i_mkbig ();
1332 mpz_init_set_si (nn1_z
, nn1
);
1333 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1334 scm_remember_upto_here_1 (n2
);
1336 return scm_i_normbig (result_z
);
1340 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1342 else if (SCM_BIGP (n1
))
1344 if (SCM_I_INUMP (n2
))
1347 nn1
= SCM_I_INUM (n1
);
1350 else if (SCM_BIGP (n2
))
1352 SCM result_z
= scm_i_mkbig ();
1353 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1355 SCM_I_BIG_MPZ (n2
));
1356 scm_remember_upto_here_2 (n1
, n2
);
1357 return scm_i_normbig (result_z
);
1360 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1363 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1368 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1370 "Return the bitwise XOR of the integer arguments. A bit is\n"
1371 "set in the result if it is set in an odd number of arguments.\n"
1373 "(logxor) @result{} 0\n"
1374 "(logxor 7) @result{} 7\n"
1375 "(logxor #b000 #b001 #b011) @result{} 2\n"
1376 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1378 #define FUNC_NAME s_scm_logxor
1382 if (SCM_UNBNDP (n2
))
1384 if (SCM_UNBNDP (n1
))
1386 else if (SCM_NUMBERP (n1
))
1389 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1392 if (SCM_I_INUMP (n1
))
1394 nn1
= SCM_I_INUM (n1
);
1395 if (SCM_I_INUMP (n2
))
1397 long nn2
= SCM_I_INUM (n2
);
1398 return SCM_I_MAKINUM (nn1
^ nn2
);
1400 else if (SCM_BIGP (n2
))
1404 SCM result_z
= scm_i_mkbig ();
1406 mpz_init_set_si (nn1_z
, nn1
);
1407 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1408 scm_remember_upto_here_1 (n2
);
1410 return scm_i_normbig (result_z
);
1414 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1416 else if (SCM_BIGP (n1
))
1418 if (SCM_I_INUMP (n2
))
1421 nn1
= SCM_I_INUM (n1
);
1424 else if (SCM_BIGP (n2
))
1426 SCM result_z
= scm_i_mkbig ();
1427 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1429 SCM_I_BIG_MPZ (n2
));
1430 scm_remember_upto_here_2 (n1
, n2
);
1431 return scm_i_normbig (result_z
);
1434 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1437 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1442 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1444 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1445 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1446 "without actually calculating the @code{logand}, just testing\n"
1450 "(logtest #b0100 #b1011) @result{} #f\n"
1451 "(logtest #b0100 #b0111) @result{} #t\n"
1453 #define FUNC_NAME s_scm_logtest
1457 if (SCM_I_INUMP (j
))
1459 nj
= SCM_I_INUM (j
);
1460 if (SCM_I_INUMP (k
))
1462 long nk
= SCM_I_INUM (k
);
1463 return scm_from_bool (nj
& nk
);
1465 else if (SCM_BIGP (k
))
1473 mpz_init_set_si (nj_z
, nj
);
1474 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1475 scm_remember_upto_here_1 (k
);
1476 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1482 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1484 else if (SCM_BIGP (j
))
1486 if (SCM_I_INUMP (k
))
1489 nj
= SCM_I_INUM (j
);
1492 else if (SCM_BIGP (k
))
1496 mpz_init (result_z
);
1500 scm_remember_upto_here_2 (j
, k
);
1501 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1502 mpz_clear (result_z
);
1506 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1509 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1514 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1516 "Test whether bit number @var{index} in @var{j} is set.\n"
1517 "@var{index} starts from 0 for the least significant bit.\n"
1520 "(logbit? 0 #b1101) @result{} #t\n"
1521 "(logbit? 1 #b1101) @result{} #f\n"
1522 "(logbit? 2 #b1101) @result{} #t\n"
1523 "(logbit? 3 #b1101) @result{} #t\n"
1524 "(logbit? 4 #b1101) @result{} #f\n"
1526 #define FUNC_NAME s_scm_logbit_p
1528 unsigned long int iindex
;
1529 iindex
= scm_to_ulong (index
);
1531 if (SCM_I_INUMP (j
))
1533 /* bits above what's in an inum follow the sign bit */
1534 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1535 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1537 else if (SCM_BIGP (j
))
1539 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1540 scm_remember_upto_here_1 (j
);
1541 return scm_from_bool (val
);
1544 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1549 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1551 "Return the integer which is the ones-complement of the integer\n"
1555 "(number->string (lognot #b10000000) 2)\n"
1556 " @result{} \"-10000001\"\n"
1557 "(number->string (lognot #b0) 2)\n"
1558 " @result{} \"-1\"\n"
1560 #define FUNC_NAME s_scm_lognot
1562 if (SCM_I_INUMP (n
)) {
1563 /* No overflow here, just need to toggle all the bits making up the inum.
1564 Enhancement: No need to strip the tag and add it back, could just xor
1565 a block of 1 bits, if that worked with the various debug versions of
1567 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1569 } else if (SCM_BIGP (n
)) {
1570 SCM result
= scm_i_mkbig ();
1571 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1572 scm_remember_upto_here_1 (n
);
1576 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1581 /* returns 0 if IN is not an integer. OUT must already be
1584 coerce_to_big (SCM in
, mpz_t out
)
1587 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1588 else if (SCM_I_INUMP (in
))
1589 mpz_set_si (out
, SCM_I_INUM (in
));
1596 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1597 (SCM n
, SCM k
, SCM m
),
1598 "Return @var{n} raised to the integer exponent\n"
1599 "@var{k}, modulo @var{m}.\n"
1602 "(modulo-expt 2 3 5)\n"
1605 #define FUNC_NAME s_scm_modulo_expt
1611 /* There are two classes of error we might encounter --
1612 1) Math errors, which we'll report by calling scm_num_overflow,
1614 2) wrong-type errors, which of course we'll report by calling
1616 We don't report those errors immediately, however; instead we do
1617 some cleanup first. These variables tell us which error (if
1618 any) we should report after cleaning up.
1620 int report_overflow
= 0;
1622 int position_of_wrong_type
= 0;
1623 SCM value_of_wrong_type
= SCM_INUM0
;
1625 SCM result
= SCM_UNDEFINED
;
1631 if (scm_is_eq (m
, SCM_INUM0
))
1633 report_overflow
= 1;
1637 if (!coerce_to_big (n
, n_tmp
))
1639 value_of_wrong_type
= n
;
1640 position_of_wrong_type
= 1;
1644 if (!coerce_to_big (k
, k_tmp
))
1646 value_of_wrong_type
= k
;
1647 position_of_wrong_type
= 2;
1651 if (!coerce_to_big (m
, m_tmp
))
1653 value_of_wrong_type
= m
;
1654 position_of_wrong_type
= 3;
1658 /* if the exponent K is negative, and we simply call mpz_powm, we
1659 will get a divide-by-zero exception when an inverse 1/n mod m
1660 doesn't exist (or is not unique). Since exceptions are hard to
1661 handle, we'll attempt the inversion "by hand" -- that way, we get
1662 a simple failure code, which is easy to handle. */
1664 if (-1 == mpz_sgn (k_tmp
))
1666 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1668 report_overflow
= 1;
1671 mpz_neg (k_tmp
, k_tmp
);
1674 result
= scm_i_mkbig ();
1675 mpz_powm (SCM_I_BIG_MPZ (result
),
1680 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1681 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1688 if (report_overflow
)
1689 scm_num_overflow (FUNC_NAME
);
1691 if (position_of_wrong_type
)
1692 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1693 value_of_wrong_type
);
1695 return scm_i_normbig (result
);
1699 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1701 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1702 "exact integer, @var{n} can be any number.\n"
1704 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1705 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1706 "includes @math{0^0} is 1.\n"
1709 "(integer-expt 2 5) @result{} 32\n"
1710 "(integer-expt -3 3) @result{} -27\n"
1711 "(integer-expt 5 -3) @result{} 1/125\n"
1712 "(integer-expt 0 0) @result{} 1\n"
1714 #define FUNC_NAME s_scm_integer_expt
1717 SCM z_i2
= SCM_BOOL_F
;
1719 SCM acc
= SCM_I_MAKINUM (1L);
1721 /* 0^0 == 1 according to R5RS */
1722 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1723 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1724 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1725 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1727 if (SCM_I_INUMP (k
))
1728 i2
= SCM_I_INUM (k
);
1729 else if (SCM_BIGP (k
))
1731 z_i2
= scm_i_clonebig (k
, 1);
1732 scm_remember_upto_here_1 (k
);
1736 SCM_WRONG_TYPE_ARG (2, k
);
1740 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1742 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1743 n
= scm_divide (n
, SCM_UNDEFINED
);
1747 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1751 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1753 return scm_product (acc
, n
);
1755 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1756 acc
= scm_product (acc
, n
);
1757 n
= scm_product (n
, n
);
1758 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1766 n
= scm_divide (n
, SCM_UNDEFINED
);
1773 return scm_product (acc
, n
);
1775 acc
= scm_product (acc
, n
);
1776 n
= scm_product (n
, n
);
1783 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1785 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1786 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1788 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1789 "@var{cnt} is negative it's a division, rounded towards negative\n"
1790 "infinity. (Note that this is not the same rounding as\n"
1791 "@code{quotient} does.)\n"
1793 "With @var{n} viewed as an infinite precision twos complement,\n"
1794 "@code{ash} means a left shift introducing zero bits, or a right\n"
1795 "shift dropping bits.\n"
1798 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1799 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1801 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1802 "(ash -23 -2) @result{} -6\n"
1804 #define FUNC_NAME s_scm_ash
1807 bits_to_shift
= scm_to_long (cnt
);
1809 if (SCM_I_INUMP (n
))
1811 long nn
= SCM_I_INUM (n
);
1813 if (bits_to_shift
> 0)
1815 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1816 overflow a non-zero fixnum. For smaller shifts we check the
1817 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1818 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1819 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1825 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1827 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1830 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1834 SCM result
= scm_i_long2big (nn
);
1835 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1842 bits_to_shift
= -bits_to_shift
;
1843 if (bits_to_shift
>= SCM_LONG_BIT
)
1844 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1846 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1850 else if (SCM_BIGP (n
))
1854 if (bits_to_shift
== 0)
1857 result
= scm_i_mkbig ();
1858 if (bits_to_shift
>= 0)
1860 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1866 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1867 we have to allocate a bignum even if the result is going to be a
1869 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1871 return scm_i_normbig (result
);
1877 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1883 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1884 (SCM n
, SCM start
, SCM end
),
1885 "Return the integer composed of the @var{start} (inclusive)\n"
1886 "through @var{end} (exclusive) bits of @var{n}. The\n"
1887 "@var{start}th bit becomes the 0-th bit in the result.\n"
1890 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1891 " @result{} \"1010\"\n"
1892 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1893 " @result{} \"10110\"\n"
1895 #define FUNC_NAME s_scm_bit_extract
1897 unsigned long int istart
, iend
, bits
;
1898 istart
= scm_to_ulong (start
);
1899 iend
= scm_to_ulong (end
);
1900 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1902 /* how many bits to keep */
1903 bits
= iend
- istart
;
1905 if (SCM_I_INUMP (n
))
1907 long int in
= SCM_I_INUM (n
);
1909 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1910 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1911 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1913 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1915 /* Since we emulate two's complement encoded numbers, this
1916 * special case requires us to produce a result that has
1917 * more bits than can be stored in a fixnum.
1919 SCM result
= scm_i_long2big (in
);
1920 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1925 /* mask down to requisite bits */
1926 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1927 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1929 else if (SCM_BIGP (n
))
1934 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1938 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1939 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1940 such bits into a ulong. */
1941 result
= scm_i_mkbig ();
1942 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1943 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1944 result
= scm_i_normbig (result
);
1946 scm_remember_upto_here_1 (n
);
1950 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1955 static const char scm_logtab
[] = {
1956 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1959 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1961 "Return the number of bits in integer @var{n}. If integer is\n"
1962 "positive, the 1-bits in its binary representation are counted.\n"
1963 "If negative, the 0-bits in its two's-complement binary\n"
1964 "representation are counted. If 0, 0 is returned.\n"
1967 "(logcount #b10101010)\n"
1974 #define FUNC_NAME s_scm_logcount
1976 if (SCM_I_INUMP (n
))
1978 unsigned long int c
= 0;
1979 long int nn
= SCM_I_INUM (n
);
1984 c
+= scm_logtab
[15 & nn
];
1987 return SCM_I_MAKINUM (c
);
1989 else if (SCM_BIGP (n
))
1991 unsigned long count
;
1992 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1993 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1995 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1996 scm_remember_upto_here_1 (n
);
1997 return SCM_I_MAKINUM (count
);
2000 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2005 static const char scm_ilentab
[] = {
2006 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2010 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2012 "Return the number of bits necessary to represent @var{n}.\n"
2015 "(integer-length #b10101010)\n"
2017 "(integer-length 0)\n"
2019 "(integer-length #b1111)\n"
2022 #define FUNC_NAME s_scm_integer_length
2024 if (SCM_I_INUMP (n
))
2026 unsigned long int c
= 0;
2028 long int nn
= SCM_I_INUM (n
);
2034 l
= scm_ilentab
[15 & nn
];
2037 return SCM_I_MAKINUM (c
- 4 + l
);
2039 else if (SCM_BIGP (n
))
2041 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2042 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2043 1 too big, so check for that and adjust. */
2044 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2045 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2046 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2047 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2049 scm_remember_upto_here_1 (n
);
2050 return SCM_I_MAKINUM (size
);
2053 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2057 /*** NUMBERS -> STRINGS ***/
2058 #define SCM_MAX_DBL_PREC 60
2059 #define SCM_MAX_DBL_RADIX 36
2061 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2062 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2063 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2066 void init_dblprec(int *prec
, int radix
) {
2067 /* determine floating point precision by adding successively
2068 smaller increments to 1.0 until it is considered == 1.0 */
2069 double f
= ((double)1.0)/radix
;
2070 double fsum
= 1.0 + f
;
2075 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2087 void init_fx_radix(double *fx_list
, int radix
)
2089 /* initialize a per-radix list of tolerances. When added
2090 to a number < 1.0, we can determine if we should raund
2091 up and quit converting a number to a string. */
2095 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2096 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2099 /* use this array as a way to generate a single digit */
2100 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2103 idbl2str (double f
, char *a
, int radix
)
2105 int efmt
, dpt
, d
, i
, wp
;
2107 #ifdef DBL_MIN_10_EXP
2110 #endif /* DBL_MIN_10_EXP */
2115 radix
> SCM_MAX_DBL_RADIX
)
2117 /* revert to existing behavior */
2121 wp
= scm_dblprec
[radix
-2];
2122 fx
= fx_per_radix
[radix
-2];
2126 #ifdef HAVE_COPYSIGN
2127 double sgn
= copysign (1.0, f
);
2132 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2138 strcpy (a
, "-inf.0");
2140 strcpy (a
, "+inf.0");
2143 else if (xisnan (f
))
2145 strcpy (a
, "+nan.0");
2155 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2156 make-uniform-vector, from causing infinite loops. */
2157 /* just do the checking...if it passes, we do the conversion for our
2158 radix again below */
2165 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2173 while (f_cpy
> 10.0)
2176 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2197 if (f
+ fx
[wp
] >= radix
)
2204 /* adding 9999 makes this equivalent to abs(x) % 3 */
2205 dpt
= (exp
+ 9999) % 3;
2209 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2231 a
[ch
++] = number_chars
[d
];
2234 if (f
+ fx
[wp
] >= 1.0)
2236 a
[ch
- 1] = number_chars
[d
+1];
2248 if ((dpt
> 4) && (exp
> 6))
2250 d
= (a
[0] == '-' ? 2 : 1);
2251 for (i
= ch
++; i
> d
; i
--)
2264 if (a
[ch
- 1] == '.')
2265 a
[ch
++] = '0'; /* trailing zero */
2274 for (i
= radix
; i
<= exp
; i
*= radix
);
2275 for (i
/= radix
; i
; i
/= radix
)
2277 a
[ch
++] = number_chars
[exp
/ i
];
2286 icmplx2str (double real
, double imag
, char *str
, int radix
)
2290 i
= idbl2str (real
, str
, radix
);
2293 /* Don't output a '+' for negative numbers or for Inf and
2294 NaN. They will provide their own sign. */
2295 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2297 i
+= idbl2str (imag
, &str
[i
], radix
);
2304 iflo2str (SCM flt
, char *str
, int radix
)
2307 if (SCM_REALP (flt
))
2308 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2310 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2315 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2316 characters in the result.
2318 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2320 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2325 return scm_iuint2str (-num
, rad
, p
) + 1;
2328 return scm_iuint2str (num
, rad
, p
);
2331 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2332 characters in the result.
2334 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2336 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2340 scm_t_uintmax n
= num
;
2342 for (n
/= rad
; n
> 0; n
/= rad
)
2352 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2357 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2359 "Return a string holding the external representation of the\n"
2360 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2361 "inexact, a radix of 10 will be used.")
2362 #define FUNC_NAME s_scm_number_to_string
2366 if (SCM_UNBNDP (radix
))
2369 base
= scm_to_signed_integer (radix
, 2, 36);
2371 if (SCM_I_INUMP (n
))
2373 char num_buf
[SCM_INTBUFLEN
];
2374 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2375 return scm_from_locale_stringn (num_buf
, length
);
2377 else if (SCM_BIGP (n
))
2379 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2380 scm_remember_upto_here_1 (n
);
2381 return scm_take_locale_string (str
);
2383 else if (SCM_FRACTIONP (n
))
2385 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2386 scm_from_locale_string ("/"),
2387 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2389 else if (SCM_INEXACTP (n
))
2391 char num_buf
[FLOBUFLEN
];
2392 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2395 SCM_WRONG_TYPE_ARG (1, n
);
2400 /* These print routines used to be stubbed here so that scm_repl.c
2401 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2404 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2406 char num_buf
[FLOBUFLEN
];
2407 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2412 scm_i_print_double (double val
, SCM port
)
2414 char num_buf
[FLOBUFLEN
];
2415 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2419 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2422 char num_buf
[FLOBUFLEN
];
2423 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2428 scm_i_print_complex (double real
, double imag
, SCM port
)
2430 char num_buf
[FLOBUFLEN
];
2431 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2435 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2438 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2439 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2440 scm_remember_upto_here_1 (str
);
2445 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2447 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2448 scm_remember_upto_here_1 (exp
);
2449 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2453 /*** END nums->strs ***/
2456 /*** STRINGS -> NUMBERS ***/
2458 /* The following functions implement the conversion from strings to numbers.
2459 * The implementation somehow follows the grammar for numbers as it is given
2460 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2461 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2462 * points should be noted about the implementation:
2463 * * Each function keeps a local index variable 'idx' that points at the
2464 * current position within the parsed string. The global index is only
2465 * updated if the function could parse the corresponding syntactic unit
2467 * * Similarly, the functions keep track of indicators of inexactness ('#',
2468 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2469 * global exactness information is only updated after each part has been
2470 * successfully parsed.
2471 * * Sequences of digits are parsed into temporary variables holding fixnums.
2472 * Only if these fixnums would overflow, the result variables are updated
2473 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2474 * the temporary variables holding the fixnums are cleared, and the process
2475 * starts over again. If for example fixnums were able to store five decimal
2476 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2477 * and the result was computed as 12345 * 100000 + 67890. In other words,
2478 * only every five digits two bignum operations were performed.
2481 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2483 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2485 /* In non ASCII-style encodings the following macro might not work. */
2486 #define XDIGIT2UINT(d) \
2487 (isdigit ((int) (unsigned char) d) \
2489 : tolower ((int) (unsigned char) d) - 'a' + 10)
2492 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2493 unsigned int radix
, enum t_exactness
*p_exactness
)
2495 unsigned int idx
= *p_idx
;
2496 unsigned int hash_seen
= 0;
2497 scm_t_bits shift
= 1;
2499 unsigned int digit_value
;
2507 if (!isxdigit ((int) (unsigned char) c
))
2509 digit_value
= XDIGIT2UINT (c
);
2510 if (digit_value
>= radix
)
2514 result
= SCM_I_MAKINUM (digit_value
);
2518 if (isxdigit ((int) (unsigned char) c
))
2522 digit_value
= XDIGIT2UINT (c
);
2523 if (digit_value
>= radix
)
2535 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2537 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2539 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2546 shift
= shift
* radix
;
2547 add
= add
* radix
+ digit_value
;
2552 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2554 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2558 *p_exactness
= INEXACT
;
2564 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2565 * covers the parts of the rules that start at a potential point. The value
2566 * of the digits up to the point have been parsed by the caller and are given
2567 * in variable result. The content of *p_exactness indicates, whether a hash
2568 * has already been seen in the digits before the point.
2571 /* In non ASCII-style encodings the following macro might not work. */
2572 #define DIGIT2UINT(d) ((d) - '0')
2575 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2576 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2578 unsigned int idx
= *p_idx
;
2579 enum t_exactness x
= *p_exactness
;
2584 if (mem
[idx
] == '.')
2586 scm_t_bits shift
= 1;
2588 unsigned int digit_value
;
2589 SCM big_shift
= SCM_I_MAKINUM (1);
2595 if (isdigit ((int) (unsigned char) c
))
2600 digit_value
= DIGIT2UINT (c
);
2611 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2613 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2614 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2616 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2624 add
= add
* 10 + digit_value
;
2630 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2631 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2632 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2635 result
= scm_divide (result
, big_shift
);
2637 /* We've seen a decimal point, thus the value is implicitly inexact. */
2649 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2676 if (!isdigit ((int) (unsigned char) c
))
2680 exponent
= DIGIT2UINT (c
);
2684 if (isdigit ((int) (unsigned char) c
))
2687 if (exponent
<= SCM_MAXEXP
)
2688 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2694 if (exponent
> SCM_MAXEXP
)
2696 size_t exp_len
= idx
- start
;
2697 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2698 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2699 scm_out_of_range ("string->number", exp_num
);
2702 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2704 result
= scm_product (result
, e
);
2706 result
= scm_divide2real (result
, e
);
2708 /* We've seen an exponent, thus the value is implicitly inexact. */
2726 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2729 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2730 unsigned int radix
, enum t_exactness
*p_exactness
)
2732 unsigned int idx
= *p_idx
;
2738 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2744 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2746 enum t_exactness x
= EXACT
;
2748 /* Cobble up the fractional part. We might want to set the
2749 NaN's mantissa from it. */
2751 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2756 if (mem
[idx
] == '.')
2760 else if (idx
+ 1 == len
)
2762 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2765 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2766 p_idx
, p_exactness
);
2770 enum t_exactness x
= EXACT
;
2773 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2774 if (scm_is_false (uinteger
))
2779 else if (mem
[idx
] == '/')
2785 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2786 if (scm_is_false (divisor
))
2789 /* both are int/big here, I assume */
2790 result
= scm_i_make_ratio (uinteger
, divisor
);
2792 else if (radix
== 10)
2794 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2795 if (scm_is_false (result
))
2806 /* When returning an inexact zero, make sure it is represented as a
2807 floating point value so that we can change its sign.
2809 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2810 result
= scm_from_double (0.0);
2816 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2819 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2820 unsigned int radix
, enum t_exactness
*p_exactness
)
2844 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2845 if (scm_is_false (ureal
))
2847 /* input must be either +i or -i */
2852 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2858 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2865 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2866 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2875 /* either +<ureal>i or -<ureal>i */
2882 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2885 /* polar input: <real>@<real>. */
2910 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2911 if (scm_is_false (angle
))
2916 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2917 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2919 result
= scm_make_polar (ureal
, angle
);
2924 /* expecting input matching <real>[+-]<ureal>?i */
2931 int sign
= (c
== '+') ? 1 : -1;
2932 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2934 if (scm_is_false (imag
))
2935 imag
= SCM_I_MAKINUM (sign
);
2936 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2937 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2941 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2948 return scm_make_rectangular (ureal
, imag
);
2957 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2959 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2962 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2963 unsigned int default_radix
)
2965 unsigned int idx
= 0;
2966 unsigned int radix
= NO_RADIX
;
2967 enum t_exactness forced_x
= NO_EXACTNESS
;
2968 enum t_exactness implicit_x
= EXACT
;
2971 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2972 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2974 switch (mem
[idx
+ 1])
2977 if (radix
!= NO_RADIX
)
2982 if (radix
!= NO_RADIX
)
2987 if (forced_x
!= NO_EXACTNESS
)
2992 if (forced_x
!= NO_EXACTNESS
)
2997 if (radix
!= NO_RADIX
)
3002 if (radix
!= NO_RADIX
)
3012 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3013 if (radix
== NO_RADIX
)
3014 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
3016 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
3018 if (scm_is_false (result
))
3024 if (SCM_INEXACTP (result
))
3025 return scm_inexact_to_exact (result
);
3029 if (SCM_INEXACTP (result
))
3032 return scm_exact_to_inexact (result
);
3035 if (implicit_x
== INEXACT
)
3037 if (SCM_INEXACTP (result
))
3040 return scm_exact_to_inexact (result
);
3048 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3049 (SCM string
, SCM radix
),
3050 "Return a number of the maximally precise representation\n"
3051 "expressed by the given @var{string}. @var{radix} must be an\n"
3052 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3053 "is a default radix that may be overridden by an explicit radix\n"
3054 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3055 "supplied, then the default radix is 10. If string is not a\n"
3056 "syntactically valid notation for a number, then\n"
3057 "@code{string->number} returns @code{#f}.")
3058 #define FUNC_NAME s_scm_string_to_number
3062 SCM_VALIDATE_STRING (1, string
);
3064 if (SCM_UNBNDP (radix
))
3067 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3069 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3070 scm_i_string_length (string
),
3072 scm_remember_upto_here_1 (string
);
3078 /*** END strs->nums ***/
3082 scm_bigequal (SCM x
, SCM y
)
3084 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3085 scm_remember_upto_here_2 (x
, y
);
3086 return scm_from_bool (0 == result
);
3090 scm_real_equalp (SCM x
, SCM y
)
3092 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3096 scm_complex_equalp (SCM x
, SCM y
)
3098 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3099 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3103 scm_i_fraction_equalp (SCM x
, SCM y
)
3105 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3106 SCM_FRACTION_NUMERATOR (y
)))
3107 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3108 SCM_FRACTION_DENOMINATOR (y
))))
3115 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3117 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3119 #define FUNC_NAME s_scm_number_p
3121 return scm_from_bool (SCM_NUMBERP (x
));
3125 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3127 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3128 "otherwise. Note that the sets of real, rational and integer\n"
3129 "values form subsets of the set of complex numbers, i. e. the\n"
3130 "predicate will also be fulfilled if @var{x} is a real,\n"
3131 "rational or integer number.")
3132 #define FUNC_NAME s_scm_complex_p
3134 /* all numbers are complex. */
3135 return scm_number_p (x
);
3139 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3141 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3142 "otherwise. Note that the set of integer values forms a subset of\n"
3143 "the set of real numbers, i. e. the predicate will also be\n"
3144 "fulfilled if @var{x} is an integer number.")
3145 #define FUNC_NAME s_scm_real_p
3147 /* we can't represent irrational numbers. */
3148 return scm_rational_p (x
);
3152 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3154 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3155 "otherwise. Note that the set of integer values forms a subset of\n"
3156 "the set of rational numbers, i. e. the predicate will also be\n"
3157 "fulfilled if @var{x} is an integer number.")
3158 #define FUNC_NAME s_scm_rational_p
3160 if (SCM_I_INUMP (x
))
3162 else if (SCM_IMP (x
))
3164 else if (SCM_BIGP (x
))
3166 else if (SCM_FRACTIONP (x
))
3168 else if (SCM_REALP (x
))
3169 /* due to their limited precision, all floating point numbers are
3170 rational as well. */
3177 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3179 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3181 #define FUNC_NAME s_scm_integer_p
3184 if (SCM_I_INUMP (x
))
3190 if (!SCM_INEXACTP (x
))
3192 if (SCM_COMPLEXP (x
))
3194 r
= SCM_REAL_VALUE (x
);
3195 /* +/-inf passes r==floor(r), making those #t */
3203 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3205 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3207 #define FUNC_NAME s_scm_inexact_p
3209 if (SCM_INEXACTP (x
))
3211 if (SCM_NUMBERP (x
))
3213 SCM_WRONG_TYPE_ARG (1, x
);
3218 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3219 /* "Return @code{#t} if all parameters are numerically equal." */
3221 scm_num_eq_p (SCM x
, SCM y
)
3224 if (SCM_I_INUMP (x
))
3226 long xx
= SCM_I_INUM (x
);
3227 if (SCM_I_INUMP (y
))
3229 long yy
= SCM_I_INUM (y
);
3230 return scm_from_bool (xx
== yy
);
3232 else if (SCM_BIGP (y
))
3234 else if (SCM_REALP (y
))
3236 /* On a 32-bit system an inum fits a double, we can cast the inum
3237 to a double and compare.
3239 But on a 64-bit system an inum is bigger than a double and
3240 casting it to a double (call that dxx) will round. dxx is at
3241 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3242 an integer and fits a long. So we cast yy to a long and
3243 compare with plain xx.
3245 An alternative (for any size system actually) would be to check
3246 yy is an integer (with floor) and is in range of an inum
3247 (compare against appropriate powers of 2) then test
3248 xx==(long)yy. It's just a matter of which casts/comparisons
3249 might be fastest or easiest for the cpu. */
3251 double yy
= SCM_REAL_VALUE (y
);
3252 return scm_from_bool ((double) xx
== yy
3253 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3254 || xx
== (long) yy
));
3256 else if (SCM_COMPLEXP (y
))
3257 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3258 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3259 else if (SCM_FRACTIONP (y
))
3262 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3264 else if (SCM_BIGP (x
))
3266 if (SCM_I_INUMP (y
))
3268 else if (SCM_BIGP (y
))
3270 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3271 scm_remember_upto_here_2 (x
, y
);
3272 return scm_from_bool (0 == cmp
);
3274 else if (SCM_REALP (y
))
3277 if (xisnan (SCM_REAL_VALUE (y
)))
3279 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3280 scm_remember_upto_here_1 (x
);
3281 return scm_from_bool (0 == cmp
);
3283 else if (SCM_COMPLEXP (y
))
3286 if (0.0 != SCM_COMPLEX_IMAG (y
))
3288 if (xisnan (SCM_COMPLEX_REAL (y
)))
3290 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3291 scm_remember_upto_here_1 (x
);
3292 return scm_from_bool (0 == cmp
);
3294 else if (SCM_FRACTIONP (y
))
3297 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3299 else if (SCM_REALP (x
))
3301 double xx
= SCM_REAL_VALUE (x
);
3302 if (SCM_I_INUMP (y
))
3304 /* see comments with inum/real above */
3305 long yy
= SCM_I_INUM (y
);
3306 return scm_from_bool (xx
== (double) yy
3307 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3308 || (long) xx
== yy
));
3310 else if (SCM_BIGP (y
))
3313 if (xisnan (SCM_REAL_VALUE (x
)))
3315 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3316 scm_remember_upto_here_1 (y
);
3317 return scm_from_bool (0 == cmp
);
3319 else if (SCM_REALP (y
))
3320 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3321 else if (SCM_COMPLEXP (y
))
3322 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3323 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3324 else if (SCM_FRACTIONP (y
))
3326 double xx
= SCM_REAL_VALUE (x
);
3330 return scm_from_bool (xx
< 0.0);
3331 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3335 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3337 else if (SCM_COMPLEXP (x
))
3339 if (SCM_I_INUMP (y
))
3340 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3341 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3342 else if (SCM_BIGP (y
))
3345 if (0.0 != SCM_COMPLEX_IMAG (x
))
3347 if (xisnan (SCM_COMPLEX_REAL (x
)))
3349 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3350 scm_remember_upto_here_1 (y
);
3351 return scm_from_bool (0 == cmp
);
3353 else if (SCM_REALP (y
))
3354 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3355 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3356 else if (SCM_COMPLEXP (y
))
3357 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3358 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3359 else if (SCM_FRACTIONP (y
))
3362 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3364 xx
= SCM_COMPLEX_REAL (x
);
3368 return scm_from_bool (xx
< 0.0);
3369 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3373 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3375 else if (SCM_FRACTIONP (x
))
3377 if (SCM_I_INUMP (y
))
3379 else if (SCM_BIGP (y
))
3381 else if (SCM_REALP (y
))
3383 double yy
= SCM_REAL_VALUE (y
);
3387 return scm_from_bool (0.0 < yy
);
3388 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3391 else if (SCM_COMPLEXP (y
))
3394 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3396 yy
= SCM_COMPLEX_REAL (y
);
3400 return scm_from_bool (0.0 < yy
);
3401 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3404 else if (SCM_FRACTIONP (y
))
3405 return scm_i_fraction_equalp (x
, y
);
3407 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3410 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3414 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3415 done are good for inums, but for bignums an answer can almost always be
3416 had by just examining a few high bits of the operands, as done by GMP in
3417 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3418 of the float exponent to take into account. */
3420 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3421 /* "Return @code{#t} if the list of parameters is monotonically\n"
3425 scm_less_p (SCM x
, SCM y
)
3428 if (SCM_I_INUMP (x
))
3430 long xx
= SCM_I_INUM (x
);
3431 if (SCM_I_INUMP (y
))
3433 long yy
= SCM_I_INUM (y
);
3434 return scm_from_bool (xx
< yy
);
3436 else if (SCM_BIGP (y
))
3438 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3439 scm_remember_upto_here_1 (y
);
3440 return scm_from_bool (sgn
> 0);
3442 else if (SCM_REALP (y
))
3443 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3444 else if (SCM_FRACTIONP (y
))
3446 /* "x < a/b" becomes "x*b < a" */
3448 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3449 y
= SCM_FRACTION_NUMERATOR (y
);
3453 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3455 else if (SCM_BIGP (x
))
3457 if (SCM_I_INUMP (y
))
3459 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3460 scm_remember_upto_here_1 (x
);
3461 return scm_from_bool (sgn
< 0);
3463 else if (SCM_BIGP (y
))
3465 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3466 scm_remember_upto_here_2 (x
, y
);
3467 return scm_from_bool (cmp
< 0);
3469 else if (SCM_REALP (y
))
3472 if (xisnan (SCM_REAL_VALUE (y
)))
3474 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3475 scm_remember_upto_here_1 (x
);
3476 return scm_from_bool (cmp
< 0);
3478 else if (SCM_FRACTIONP (y
))
3481 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3483 else if (SCM_REALP (x
))
3485 if (SCM_I_INUMP (y
))
3486 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3487 else if (SCM_BIGP (y
))
3490 if (xisnan (SCM_REAL_VALUE (x
)))
3492 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3493 scm_remember_upto_here_1 (y
);
3494 return scm_from_bool (cmp
> 0);
3496 else if (SCM_REALP (y
))
3497 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3498 else if (SCM_FRACTIONP (y
))
3500 double xx
= SCM_REAL_VALUE (x
);
3504 return scm_from_bool (xx
< 0.0);
3505 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3509 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3511 else if (SCM_FRACTIONP (x
))
3513 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3515 /* "a/b < y" becomes "a < y*b" */
3516 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3517 x
= SCM_FRACTION_NUMERATOR (x
);
3520 else if (SCM_REALP (y
))
3522 double yy
= SCM_REAL_VALUE (y
);
3526 return scm_from_bool (0.0 < yy
);
3527 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3530 else if (SCM_FRACTIONP (y
))
3532 /* "a/b < c/d" becomes "a*d < c*b" */
3533 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3534 SCM_FRACTION_DENOMINATOR (y
));
3535 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3536 SCM_FRACTION_DENOMINATOR (x
));
3542 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3545 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3549 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3550 /* "Return @code{#t} if the list of parameters is monotonically\n"
3553 #define FUNC_NAME s_scm_gr_p
3555 scm_gr_p (SCM x
, SCM y
)
3557 if (!SCM_NUMBERP (x
))
3558 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3559 else if (!SCM_NUMBERP (y
))
3560 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3562 return scm_less_p (y
, x
);
3567 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3568 /* "Return @code{#t} if the list of parameters is monotonically\n"
3571 #define FUNC_NAME s_scm_leq_p
3573 scm_leq_p (SCM x
, SCM y
)
3575 if (!SCM_NUMBERP (x
))
3576 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3577 else if (!SCM_NUMBERP (y
))
3578 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3579 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3582 return scm_not (scm_less_p (y
, x
));
3587 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3588 /* "Return @code{#t} if the list of parameters is monotonically\n"
3591 #define FUNC_NAME s_scm_geq_p
3593 scm_geq_p (SCM x
, SCM y
)
3595 if (!SCM_NUMBERP (x
))
3596 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3597 else if (!SCM_NUMBERP (y
))
3598 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3599 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3602 return scm_not (scm_less_p (x
, y
));
3607 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3608 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3614 if (SCM_I_INUMP (z
))
3615 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3616 else if (SCM_BIGP (z
))
3618 else if (SCM_REALP (z
))
3619 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3620 else if (SCM_COMPLEXP (z
))
3621 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3622 && SCM_COMPLEX_IMAG (z
) == 0.0);
3623 else if (SCM_FRACTIONP (z
))
3626 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3630 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3631 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3635 scm_positive_p (SCM x
)
3637 if (SCM_I_INUMP (x
))
3638 return scm_from_bool (SCM_I_INUM (x
) > 0);
3639 else if (SCM_BIGP (x
))
3641 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3642 scm_remember_upto_here_1 (x
);
3643 return scm_from_bool (sgn
> 0);
3645 else if (SCM_REALP (x
))
3646 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3647 else if (SCM_FRACTIONP (x
))
3648 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3650 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3654 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3655 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3659 scm_negative_p (SCM x
)
3661 if (SCM_I_INUMP (x
))
3662 return scm_from_bool (SCM_I_INUM (x
) < 0);
3663 else if (SCM_BIGP (x
))
3665 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3666 scm_remember_upto_here_1 (x
);
3667 return scm_from_bool (sgn
< 0);
3669 else if (SCM_REALP (x
))
3670 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3671 else if (SCM_FRACTIONP (x
))
3672 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3674 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3678 /* scm_min and scm_max return an inexact when either argument is inexact, as
3679 required by r5rs. On that basis, for exact/inexact combinations the
3680 exact is converted to inexact to compare and possibly return. This is
3681 unlike scm_less_p above which takes some trouble to preserve all bits in
3682 its test, such trouble is not required for min and max. */
3684 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3685 /* "Return the maximum of all parameter values."
3688 scm_max (SCM x
, SCM y
)
3693 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3694 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3697 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3700 if (SCM_I_INUMP (x
))
3702 long xx
= SCM_I_INUM (x
);
3703 if (SCM_I_INUMP (y
))
3705 long yy
= SCM_I_INUM (y
);
3706 return (xx
< yy
) ? y
: x
;
3708 else if (SCM_BIGP (y
))
3710 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3711 scm_remember_upto_here_1 (y
);
3712 return (sgn
< 0) ? x
: y
;
3714 else if (SCM_REALP (y
))
3717 /* if y==NaN then ">" is false and we return NaN */
3718 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3720 else if (SCM_FRACTIONP (y
))
3723 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3726 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3728 else if (SCM_BIGP (x
))
3730 if (SCM_I_INUMP (y
))
3732 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3733 scm_remember_upto_here_1 (x
);
3734 return (sgn
< 0) ? y
: x
;
3736 else if (SCM_BIGP (y
))
3738 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3739 scm_remember_upto_here_2 (x
, y
);
3740 return (cmp
> 0) ? x
: y
;
3742 else if (SCM_REALP (y
))
3744 /* if y==NaN then xx>yy is false, so we return the NaN y */
3747 xx
= scm_i_big2dbl (x
);
3748 yy
= SCM_REAL_VALUE (y
);
3749 return (xx
> yy
? scm_from_double (xx
) : y
);
3751 else if (SCM_FRACTIONP (y
))
3756 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3758 else if (SCM_REALP (x
))
3760 if (SCM_I_INUMP (y
))
3762 double z
= SCM_I_INUM (y
);
3763 /* if x==NaN then "<" is false and we return NaN */
3764 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3766 else if (SCM_BIGP (y
))
3771 else if (SCM_REALP (y
))
3773 /* if x==NaN then our explicit check means we return NaN
3774 if y==NaN then ">" is false and we return NaN
3775 calling isnan is unavoidable, since it's the only way to know
3776 which of x or y causes any compares to be false */
3777 double xx
= SCM_REAL_VALUE (x
);
3778 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3780 else if (SCM_FRACTIONP (y
))
3782 double yy
= scm_i_fraction2double (y
);
3783 double xx
= SCM_REAL_VALUE (x
);
3784 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3787 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3789 else if (SCM_FRACTIONP (x
))
3791 if (SCM_I_INUMP (y
))
3795 else if (SCM_BIGP (y
))
3799 else if (SCM_REALP (y
))
3801 double xx
= scm_i_fraction2double (x
);
3802 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3804 else if (SCM_FRACTIONP (y
))
3809 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3812 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3816 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3817 /* "Return the minium of all parameter values."
3820 scm_min (SCM x
, SCM y
)
3825 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3826 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3829 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3832 if (SCM_I_INUMP (x
))
3834 long xx
= SCM_I_INUM (x
);
3835 if (SCM_I_INUMP (y
))
3837 long yy
= SCM_I_INUM (y
);
3838 return (xx
< yy
) ? x
: y
;
3840 else if (SCM_BIGP (y
))
3842 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3843 scm_remember_upto_here_1 (y
);
3844 return (sgn
< 0) ? y
: x
;
3846 else if (SCM_REALP (y
))
3849 /* if y==NaN then "<" is false and we return NaN */
3850 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3852 else if (SCM_FRACTIONP (y
))
3855 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3858 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3860 else if (SCM_BIGP (x
))
3862 if (SCM_I_INUMP (y
))
3864 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3865 scm_remember_upto_here_1 (x
);
3866 return (sgn
< 0) ? x
: y
;
3868 else if (SCM_BIGP (y
))
3870 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3871 scm_remember_upto_here_2 (x
, y
);
3872 return (cmp
> 0) ? y
: x
;
3874 else if (SCM_REALP (y
))
3876 /* if y==NaN then xx<yy is false, so we return the NaN y */
3879 xx
= scm_i_big2dbl (x
);
3880 yy
= SCM_REAL_VALUE (y
);
3881 return (xx
< yy
? scm_from_double (xx
) : y
);
3883 else if (SCM_FRACTIONP (y
))
3888 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3890 else if (SCM_REALP (x
))
3892 if (SCM_I_INUMP (y
))
3894 double z
= SCM_I_INUM (y
);
3895 /* if x==NaN then "<" is false and we return NaN */
3896 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3898 else if (SCM_BIGP (y
))
3903 else if (SCM_REALP (y
))
3905 /* if x==NaN then our explicit check means we return NaN
3906 if y==NaN then "<" is false and we return NaN
3907 calling isnan is unavoidable, since it's the only way to know
3908 which of x or y causes any compares to be false */
3909 double xx
= SCM_REAL_VALUE (x
);
3910 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3912 else if (SCM_FRACTIONP (y
))
3914 double yy
= scm_i_fraction2double (y
);
3915 double xx
= SCM_REAL_VALUE (x
);
3916 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3919 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3921 else if (SCM_FRACTIONP (x
))
3923 if (SCM_I_INUMP (y
))
3927 else if (SCM_BIGP (y
))
3931 else if (SCM_REALP (y
))
3933 double xx
= scm_i_fraction2double (x
);
3934 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3936 else if (SCM_FRACTIONP (y
))
3941 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3944 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3948 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3949 /* "Return the sum of all parameter values. Return 0 if called without\n"
3953 scm_sum (SCM x
, SCM y
)
3955 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3957 if (SCM_NUMBERP (x
)) return x
;
3958 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3959 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3962 if (SCM_LIKELY (SCM_I_INUMP (x
)))
3964 if (SCM_LIKELY (SCM_I_INUMP (y
)))
3966 long xx
= SCM_I_INUM (x
);
3967 long yy
= SCM_I_INUM (y
);
3968 long int z
= xx
+ yy
;
3969 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3971 else if (SCM_BIGP (y
))
3976 else if (SCM_REALP (y
))
3978 long int xx
= SCM_I_INUM (x
);
3979 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3981 else if (SCM_COMPLEXP (y
))
3983 long int xx
= SCM_I_INUM (x
);
3984 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3985 SCM_COMPLEX_IMAG (y
));
3987 else if (SCM_FRACTIONP (y
))
3988 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3989 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3990 SCM_FRACTION_DENOMINATOR (y
));
3992 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3993 } else if (SCM_BIGP (x
))
3995 if (SCM_I_INUMP (y
))
4000 inum
= SCM_I_INUM (y
);
4003 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4006 SCM result
= scm_i_mkbig ();
4007 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4008 scm_remember_upto_here_1 (x
);
4009 /* we know the result will have to be a bignum */
4012 return scm_i_normbig (result
);
4016 SCM result
= scm_i_mkbig ();
4017 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4018 scm_remember_upto_here_1 (x
);
4019 /* we know the result will have to be a bignum */
4022 return scm_i_normbig (result
);
4025 else if (SCM_BIGP (y
))
4027 SCM result
= scm_i_mkbig ();
4028 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4029 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4030 mpz_add (SCM_I_BIG_MPZ (result
),
4033 scm_remember_upto_here_2 (x
, y
);
4034 /* we know the result will have to be a bignum */
4037 return scm_i_normbig (result
);
4039 else if (SCM_REALP (y
))
4041 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4042 scm_remember_upto_here_1 (x
);
4043 return scm_from_double (result
);
4045 else if (SCM_COMPLEXP (y
))
4047 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4048 + SCM_COMPLEX_REAL (y
));
4049 scm_remember_upto_here_1 (x
);
4050 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4052 else if (SCM_FRACTIONP (y
))
4053 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4054 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4055 SCM_FRACTION_DENOMINATOR (y
));
4057 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4059 else if (SCM_REALP (x
))
4061 if (SCM_I_INUMP (y
))
4062 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4063 else if (SCM_BIGP (y
))
4065 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4066 scm_remember_upto_here_1 (y
);
4067 return scm_from_double (result
);
4069 else if (SCM_REALP (y
))
4070 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4071 else if (SCM_COMPLEXP (y
))
4072 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4073 SCM_COMPLEX_IMAG (y
));
4074 else if (SCM_FRACTIONP (y
))
4075 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4077 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4079 else if (SCM_COMPLEXP (x
))
4081 if (SCM_I_INUMP (y
))
4082 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4083 SCM_COMPLEX_IMAG (x
));
4084 else if (SCM_BIGP (y
))
4086 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4087 + SCM_COMPLEX_REAL (x
));
4088 scm_remember_upto_here_1 (y
);
4089 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4091 else if (SCM_REALP (y
))
4092 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4093 SCM_COMPLEX_IMAG (x
));
4094 else if (SCM_COMPLEXP (y
))
4095 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4096 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4097 else if (SCM_FRACTIONP (y
))
4098 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4099 SCM_COMPLEX_IMAG (x
));
4101 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4103 else if (SCM_FRACTIONP (x
))
4105 if (SCM_I_INUMP (y
))
4106 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4107 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4108 SCM_FRACTION_DENOMINATOR (x
));
4109 else if (SCM_BIGP (y
))
4110 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4111 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4112 SCM_FRACTION_DENOMINATOR (x
));
4113 else if (SCM_REALP (y
))
4114 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4115 else if (SCM_COMPLEXP (y
))
4116 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4117 SCM_COMPLEX_IMAG (y
));
4118 else if (SCM_FRACTIONP (y
))
4119 /* a/b + c/d = (ad + bc) / bd */
4120 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4121 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4122 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4124 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4127 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4131 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4133 "Return @math{@var{x}+1}.")
4134 #define FUNC_NAME s_scm_oneplus
4136 return scm_sum (x
, SCM_I_MAKINUM (1));
4141 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4142 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4143 * the sum of all but the first argument are subtracted from the first
4145 #define FUNC_NAME s_difference
4147 scm_difference (SCM x
, SCM y
)
4149 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4152 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4154 if (SCM_I_INUMP (x
))
4156 long xx
= -SCM_I_INUM (x
);
4157 if (SCM_FIXABLE (xx
))
4158 return SCM_I_MAKINUM (xx
);
4160 return scm_i_long2big (xx
);
4162 else if (SCM_BIGP (x
))
4163 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4164 bignum, but negating that gives a fixnum. */
4165 return scm_i_normbig (scm_i_clonebig (x
, 0));
4166 else if (SCM_REALP (x
))
4167 return scm_from_double (-SCM_REAL_VALUE (x
));
4168 else if (SCM_COMPLEXP (x
))
4169 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4170 -SCM_COMPLEX_IMAG (x
));
4171 else if (SCM_FRACTIONP (x
))
4172 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4173 SCM_FRACTION_DENOMINATOR (x
));
4175 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4178 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4180 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4182 long int xx
= SCM_I_INUM (x
);
4183 long int yy
= SCM_I_INUM (y
);
4184 long int z
= xx
- yy
;
4185 if (SCM_FIXABLE (z
))
4186 return SCM_I_MAKINUM (z
);
4188 return scm_i_long2big (z
);
4190 else if (SCM_BIGP (y
))
4192 /* inum-x - big-y */
4193 long xx
= SCM_I_INUM (x
);
4196 return scm_i_clonebig (y
, 0);
4199 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4200 SCM result
= scm_i_mkbig ();
4203 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4206 /* x - y == -(y + -x) */
4207 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4208 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4210 scm_remember_upto_here_1 (y
);
4212 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4213 /* we know the result will have to be a bignum */
4216 return scm_i_normbig (result
);
4219 else if (SCM_REALP (y
))
4221 long int xx
= SCM_I_INUM (x
);
4222 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4224 else if (SCM_COMPLEXP (y
))
4226 long int xx
= SCM_I_INUM (x
);
4227 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4228 - SCM_COMPLEX_IMAG (y
));
4230 else if (SCM_FRACTIONP (y
))
4231 /* a - b/c = (ac - b) / c */
4232 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4233 SCM_FRACTION_NUMERATOR (y
)),
4234 SCM_FRACTION_DENOMINATOR (y
));
4236 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4238 else if (SCM_BIGP (x
))
4240 if (SCM_I_INUMP (y
))
4242 /* big-x - inum-y */
4243 long yy
= SCM_I_INUM (y
);
4244 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4246 scm_remember_upto_here_1 (x
);
4248 return (SCM_FIXABLE (-yy
) ?
4249 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4252 SCM result
= scm_i_mkbig ();
4255 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4257 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4258 scm_remember_upto_here_1 (x
);
4260 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4261 /* we know the result will have to be a bignum */
4264 return scm_i_normbig (result
);
4267 else if (SCM_BIGP (y
))
4269 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4270 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4271 SCM result
= scm_i_mkbig ();
4272 mpz_sub (SCM_I_BIG_MPZ (result
),
4275 scm_remember_upto_here_2 (x
, y
);
4276 /* we know the result will have to be a bignum */
4277 if ((sgn_x
== 1) && (sgn_y
== -1))
4279 if ((sgn_x
== -1) && (sgn_y
== 1))
4281 return scm_i_normbig (result
);
4283 else if (SCM_REALP (y
))
4285 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4286 scm_remember_upto_here_1 (x
);
4287 return scm_from_double (result
);
4289 else if (SCM_COMPLEXP (y
))
4291 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4292 - SCM_COMPLEX_REAL (y
));
4293 scm_remember_upto_here_1 (x
);
4294 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4296 else if (SCM_FRACTIONP (y
))
4297 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4298 SCM_FRACTION_NUMERATOR (y
)),
4299 SCM_FRACTION_DENOMINATOR (y
));
4300 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4302 else if (SCM_REALP (x
))
4304 if (SCM_I_INUMP (y
))
4305 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4306 else if (SCM_BIGP (y
))
4308 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4309 scm_remember_upto_here_1 (x
);
4310 return scm_from_double (result
);
4312 else if (SCM_REALP (y
))
4313 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4314 else if (SCM_COMPLEXP (y
))
4315 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4316 -SCM_COMPLEX_IMAG (y
));
4317 else if (SCM_FRACTIONP (y
))
4318 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4320 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4322 else if (SCM_COMPLEXP (x
))
4324 if (SCM_I_INUMP (y
))
4325 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4326 SCM_COMPLEX_IMAG (x
));
4327 else if (SCM_BIGP (y
))
4329 double real_part
= (SCM_COMPLEX_REAL (x
)
4330 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4331 scm_remember_upto_here_1 (x
);
4332 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4334 else if (SCM_REALP (y
))
4335 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4336 SCM_COMPLEX_IMAG (x
));
4337 else if (SCM_COMPLEXP (y
))
4338 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4339 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4340 else if (SCM_FRACTIONP (y
))
4341 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4342 SCM_COMPLEX_IMAG (x
));
4344 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4346 else if (SCM_FRACTIONP (x
))
4348 if (SCM_I_INUMP (y
))
4349 /* a/b - c = (a - cb) / b */
4350 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4351 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4352 SCM_FRACTION_DENOMINATOR (x
));
4353 else if (SCM_BIGP (y
))
4354 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4355 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4356 SCM_FRACTION_DENOMINATOR (x
));
4357 else if (SCM_REALP (y
))
4358 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4359 else if (SCM_COMPLEXP (y
))
4360 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4361 -SCM_COMPLEX_IMAG (y
));
4362 else if (SCM_FRACTIONP (y
))
4363 /* a/b - c/d = (ad - bc) / bd */
4364 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4365 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4366 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4368 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4371 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4376 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4378 "Return @math{@var{x}-1}.")
4379 #define FUNC_NAME s_scm_oneminus
4381 return scm_difference (x
, SCM_I_MAKINUM (1));
4386 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4387 /* "Return the product of all arguments. If called without arguments,\n"
4391 scm_product (SCM x
, SCM y
)
4393 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4396 return SCM_I_MAKINUM (1L);
4397 else if (SCM_NUMBERP (x
))
4400 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4403 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4408 xx
= SCM_I_INUM (x
);
4412 case 0: return x
; break;
4413 case 1: return y
; break;
4416 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4418 long yy
= SCM_I_INUM (y
);
4420 SCM k
= SCM_I_MAKINUM (kk
);
4421 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4425 SCM result
= scm_i_long2big (xx
);
4426 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4427 return scm_i_normbig (result
);
4430 else if (SCM_BIGP (y
))
4432 SCM result
= scm_i_mkbig ();
4433 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4434 scm_remember_upto_here_1 (y
);
4437 else if (SCM_REALP (y
))
4438 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4439 else if (SCM_COMPLEXP (y
))
4440 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4441 xx
* SCM_COMPLEX_IMAG (y
));
4442 else if (SCM_FRACTIONP (y
))
4443 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4444 SCM_FRACTION_DENOMINATOR (y
));
4446 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4448 else if (SCM_BIGP (x
))
4450 if (SCM_I_INUMP (y
))
4455 else if (SCM_BIGP (y
))
4457 SCM result
= scm_i_mkbig ();
4458 mpz_mul (SCM_I_BIG_MPZ (result
),
4461 scm_remember_upto_here_2 (x
, y
);
4464 else if (SCM_REALP (y
))
4466 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4467 scm_remember_upto_here_1 (x
);
4468 return scm_from_double (result
);
4470 else if (SCM_COMPLEXP (y
))
4472 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4473 scm_remember_upto_here_1 (x
);
4474 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4475 z
* SCM_COMPLEX_IMAG (y
));
4477 else if (SCM_FRACTIONP (y
))
4478 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4479 SCM_FRACTION_DENOMINATOR (y
));
4481 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4483 else if (SCM_REALP (x
))
4485 if (SCM_I_INUMP (y
))
4487 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4488 if (scm_is_eq (y
, SCM_INUM0
))
4490 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4492 else if (SCM_BIGP (y
))
4494 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4495 scm_remember_upto_here_1 (y
);
4496 return scm_from_double (result
);
4498 else if (SCM_REALP (y
))
4499 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4500 else if (SCM_COMPLEXP (y
))
4501 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4502 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4503 else if (SCM_FRACTIONP (y
))
4504 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4506 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4508 else if (SCM_COMPLEXP (x
))
4510 if (SCM_I_INUMP (y
))
4512 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4513 if (scm_is_eq (y
, SCM_INUM0
))
4515 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4516 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4518 else if (SCM_BIGP (y
))
4520 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4521 scm_remember_upto_here_1 (y
);
4522 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4523 z
* SCM_COMPLEX_IMAG (x
));
4525 else if (SCM_REALP (y
))
4526 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4527 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4528 else if (SCM_COMPLEXP (y
))
4530 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4531 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4532 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4533 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4535 else if (SCM_FRACTIONP (y
))
4537 double yy
= scm_i_fraction2double (y
);
4538 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4539 yy
* SCM_COMPLEX_IMAG (x
));
4542 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4544 else if (SCM_FRACTIONP (x
))
4546 if (SCM_I_INUMP (y
))
4547 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4548 SCM_FRACTION_DENOMINATOR (x
));
4549 else if (SCM_BIGP (y
))
4550 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4551 SCM_FRACTION_DENOMINATOR (x
));
4552 else if (SCM_REALP (y
))
4553 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4554 else if (SCM_COMPLEXP (y
))
4556 double xx
= scm_i_fraction2double (x
);
4557 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4558 xx
* SCM_COMPLEX_IMAG (y
));
4560 else if (SCM_FRACTIONP (y
))
4561 /* a/b * c/d = ac / bd */
4562 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4563 SCM_FRACTION_NUMERATOR (y
)),
4564 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4565 SCM_FRACTION_DENOMINATOR (y
)));
4567 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4570 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4573 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4574 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4575 #define ALLOW_DIVIDE_BY_ZERO
4576 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4579 /* The code below for complex division is adapted from the GNU
4580 libstdc++, which adapted it from f2c's libF77, and is subject to
4583 /****************************************************************
4584 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4586 Permission to use, copy, modify, and distribute this software
4587 and its documentation for any purpose and without fee is hereby
4588 granted, provided that the above copyright notice appear in all
4589 copies and that both that the copyright notice and this
4590 permission notice and warranty disclaimer appear in supporting
4591 documentation, and that the names of AT&T Bell Laboratories or
4592 Bellcore or any of their entities not be used in advertising or
4593 publicity pertaining to distribution of the software without
4594 specific, written prior permission.
4596 AT&T and Bellcore disclaim all warranties with regard to this
4597 software, including all implied warranties of merchantability
4598 and fitness. In no event shall AT&T or Bellcore be liable for
4599 any special, indirect or consequential damages or any damages
4600 whatsoever resulting from loss of use, data or profits, whether
4601 in an action of contract, negligence or other tortious action,
4602 arising out of or in connection with the use or performance of
4604 ****************************************************************/
4606 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4607 /* Divide the first argument by the product of the remaining
4608 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4610 #define FUNC_NAME s_divide
4612 scm_i_divide (SCM x
, SCM y
, int inexact
)
4616 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4619 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4620 else if (SCM_I_INUMP (x
))
4622 long xx
= SCM_I_INUM (x
);
4623 if (xx
== 1 || xx
== -1)
4625 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4627 scm_num_overflow (s_divide
);
4632 return scm_from_double (1.0 / (double) xx
);
4633 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4636 else if (SCM_BIGP (x
))
4639 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4640 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4642 else if (SCM_REALP (x
))
4644 double xx
= SCM_REAL_VALUE (x
);
4645 #ifndef ALLOW_DIVIDE_BY_ZERO
4647 scm_num_overflow (s_divide
);
4650 return scm_from_double (1.0 / xx
);
4652 else if (SCM_COMPLEXP (x
))
4654 double r
= SCM_COMPLEX_REAL (x
);
4655 double i
= SCM_COMPLEX_IMAG (x
);
4656 if (fabs(r
) <= fabs(i
))
4659 double d
= i
* (1.0 + t
* t
);
4660 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4665 double d
= r
* (1.0 + t
* t
);
4666 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4669 else if (SCM_FRACTIONP (x
))
4670 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4671 SCM_FRACTION_NUMERATOR (x
));
4673 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4676 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4678 long xx
= SCM_I_INUM (x
);
4679 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4681 long yy
= SCM_I_INUM (y
);
4684 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4685 scm_num_overflow (s_divide
);
4687 return scm_from_double ((double) xx
/ (double) yy
);
4690 else if (xx
% yy
!= 0)
4693 return scm_from_double ((double) xx
/ (double) yy
);
4694 else return scm_i_make_ratio (x
, y
);
4699 if (SCM_FIXABLE (z
))
4700 return SCM_I_MAKINUM (z
);
4702 return scm_i_long2big (z
);
4705 else if (SCM_BIGP (y
))
4708 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4709 else return scm_i_make_ratio (x
, y
);
4711 else if (SCM_REALP (y
))
4713 double yy
= SCM_REAL_VALUE (y
);
4714 #ifndef ALLOW_DIVIDE_BY_ZERO
4716 scm_num_overflow (s_divide
);
4719 return scm_from_double ((double) xx
/ yy
);
4721 else if (SCM_COMPLEXP (y
))
4724 complex_div
: /* y _must_ be a complex number */
4726 double r
= SCM_COMPLEX_REAL (y
);
4727 double i
= SCM_COMPLEX_IMAG (y
);
4728 if (fabs(r
) <= fabs(i
))
4731 double d
= i
* (1.0 + t
* t
);
4732 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4737 double d
= r
* (1.0 + t
* t
);
4738 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4742 else if (SCM_FRACTIONP (y
))
4743 /* a / b/c = ac / b */
4744 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4745 SCM_FRACTION_NUMERATOR (y
));
4747 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4749 else if (SCM_BIGP (x
))
4751 if (SCM_I_INUMP (y
))
4753 long int yy
= SCM_I_INUM (y
);
4756 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4757 scm_num_overflow (s_divide
);
4759 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4760 scm_remember_upto_here_1 (x
);
4761 return (sgn
== 0) ? scm_nan () : scm_inf ();
4768 /* FIXME: HMM, what are the relative performance issues here?
4769 We need to test. Is it faster on average to test
4770 divisible_p, then perform whichever operation, or is it
4771 faster to perform the integer div opportunistically and
4772 switch to real if there's a remainder? For now we take the
4773 middle ground: test, then if divisible, use the faster div
4776 long abs_yy
= yy
< 0 ? -yy
: yy
;
4777 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4781 SCM result
= scm_i_mkbig ();
4782 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4783 scm_remember_upto_here_1 (x
);
4785 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4786 return scm_i_normbig (result
);
4791 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4792 else return scm_i_make_ratio (x
, y
);
4796 else if (SCM_BIGP (y
))
4798 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4801 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4802 scm_num_overflow (s_divide
);
4804 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4805 scm_remember_upto_here_1 (x
);
4806 return (sgn
== 0) ? scm_nan () : scm_inf ();
4814 /* It's easily possible for the ratio x/y to fit a double
4815 but one or both x and y be too big to fit a double,
4816 hence the use of mpq_get_d rather than converting and
4819 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4820 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4821 return scm_from_double (mpq_get_d (q
));
4825 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4829 SCM result
= scm_i_mkbig ();
4830 mpz_divexact (SCM_I_BIG_MPZ (result
),
4833 scm_remember_upto_here_2 (x
, y
);
4834 return scm_i_normbig (result
);
4837 return scm_i_make_ratio (x
, y
);
4841 else if (SCM_REALP (y
))
4843 double yy
= SCM_REAL_VALUE (y
);
4844 #ifndef ALLOW_DIVIDE_BY_ZERO
4846 scm_num_overflow (s_divide
);
4849 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4851 else if (SCM_COMPLEXP (y
))
4853 a
= scm_i_big2dbl (x
);
4856 else if (SCM_FRACTIONP (y
))
4857 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4858 SCM_FRACTION_NUMERATOR (y
));
4860 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4862 else if (SCM_REALP (x
))
4864 double rx
= SCM_REAL_VALUE (x
);
4865 if (SCM_I_INUMP (y
))
4867 long int yy
= SCM_I_INUM (y
);
4868 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4870 scm_num_overflow (s_divide
);
4873 return scm_from_double (rx
/ (double) yy
);
4875 else if (SCM_BIGP (y
))
4877 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4878 scm_remember_upto_here_1 (y
);
4879 return scm_from_double (rx
/ dby
);
4881 else if (SCM_REALP (y
))
4883 double yy
= SCM_REAL_VALUE (y
);
4884 #ifndef ALLOW_DIVIDE_BY_ZERO
4886 scm_num_overflow (s_divide
);
4889 return scm_from_double (rx
/ yy
);
4891 else if (SCM_COMPLEXP (y
))
4896 else if (SCM_FRACTIONP (y
))
4897 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4899 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4901 else if (SCM_COMPLEXP (x
))
4903 double rx
= SCM_COMPLEX_REAL (x
);
4904 double ix
= SCM_COMPLEX_IMAG (x
);
4905 if (SCM_I_INUMP (y
))
4907 long int yy
= SCM_I_INUM (y
);
4908 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4910 scm_num_overflow (s_divide
);
4915 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4918 else if (SCM_BIGP (y
))
4920 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4921 scm_remember_upto_here_1 (y
);
4922 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4924 else if (SCM_REALP (y
))
4926 double yy
= SCM_REAL_VALUE (y
);
4927 #ifndef ALLOW_DIVIDE_BY_ZERO
4929 scm_num_overflow (s_divide
);
4932 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4934 else if (SCM_COMPLEXP (y
))
4936 double ry
= SCM_COMPLEX_REAL (y
);
4937 double iy
= SCM_COMPLEX_IMAG (y
);
4938 if (fabs(ry
) <= fabs(iy
))
4941 double d
= iy
* (1.0 + t
* t
);
4942 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4947 double d
= ry
* (1.0 + t
* t
);
4948 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4951 else if (SCM_FRACTIONP (y
))
4953 double yy
= scm_i_fraction2double (y
);
4954 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4957 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4959 else if (SCM_FRACTIONP (x
))
4961 if (SCM_I_INUMP (y
))
4963 long int yy
= SCM_I_INUM (y
);
4964 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4966 scm_num_overflow (s_divide
);
4969 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4970 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4972 else if (SCM_BIGP (y
))
4974 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4975 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4977 else if (SCM_REALP (y
))
4979 double yy
= SCM_REAL_VALUE (y
);
4980 #ifndef ALLOW_DIVIDE_BY_ZERO
4982 scm_num_overflow (s_divide
);
4985 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4987 else if (SCM_COMPLEXP (y
))
4989 a
= scm_i_fraction2double (x
);
4992 else if (SCM_FRACTIONP (y
))
4993 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4994 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4996 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4999 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5003 scm_divide (SCM x
, SCM y
)
5005 return scm_i_divide (x
, y
, 0);
5008 static SCM
scm_divide2real (SCM x
, SCM y
)
5010 return scm_i_divide (x
, y
, 1);
5016 scm_asinh (double x
)
5021 #define asinh scm_asinh
5022 return log (x
+ sqrt (x
* x
+ 1));
5025 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5026 /* "Return the inverse hyperbolic sine of @var{x}."
5031 scm_acosh (double x
)
5036 #define acosh scm_acosh
5037 return log (x
+ sqrt (x
* x
- 1));
5040 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5041 /* "Return the inverse hyperbolic cosine of @var{x}."
5046 scm_atanh (double x
)
5051 #define atanh scm_atanh
5052 return 0.5 * log ((1 + x
) / (1 - x
));
5055 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5056 /* "Return the inverse hyperbolic tangent of @var{x}."
5061 scm_c_truncate (double x
)
5072 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5073 half-way case (ie. when x is an integer plus 0.5) going upwards.
5074 Then half-way cases are identified and adjusted down if the
5075 round-upwards didn't give the desired even integer.
5077 "plus_half == result" identifies a half-way case. If plus_half, which is
5078 x + 0.5, is an integer then x must be an integer plus 0.5.
5080 An odd "result" value is identified with result/2 != floor(result/2).
5081 This is done with plus_half, since that value is ready for use sooner in
5082 a pipelined cpu, and we're already requiring plus_half == result.
5084 Note however that we need to be careful when x is big and already an
5085 integer. In that case "x+0.5" may round to an adjacent integer, causing
5086 us to return such a value, incorrectly. For instance if the hardware is
5087 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5088 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5089 returned. Or if the hardware is in round-upwards mode, then other bigger
5090 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5091 representable value, 2^128+2^76 (or whatever), again incorrect.
5093 These bad roundings of x+0.5 are avoided by testing at the start whether
5094 x is already an integer. If it is then clearly that's the desired result
5095 already. And if it's not then the exponent must be small enough to allow
5096 an 0.5 to be represented, and hence added without a bad rounding. */
5099 scm_c_round (double x
)
5101 double plus_half
, result
;
5106 plus_half
= x
+ 0.5;
5107 result
= floor (plus_half
);
5108 /* Adjust so that the rounding is towards even. */
5109 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5114 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5116 "Round the number @var{x} towards zero.")
5117 #define FUNC_NAME s_scm_truncate_number
5119 if (scm_is_false (scm_negative_p (x
)))
5120 return scm_floor (x
);
5122 return scm_ceiling (x
);
5126 static SCM exactly_one_half
;
5128 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5130 "Round the number @var{x} towards the nearest integer. "
5131 "When it is exactly halfway between two integers, "
5132 "round towards the even one.")
5133 #define FUNC_NAME s_scm_round_number
5135 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5137 else if (SCM_REALP (x
))
5138 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5141 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5142 single quotient+remainder division then examining to see which way
5143 the rounding should go. */
5144 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5145 SCM result
= scm_floor (plus_half
);
5146 /* Adjust so that the rounding is towards even. */
5147 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5148 && scm_is_true (scm_odd_p (result
)))
5149 return scm_difference (result
, SCM_I_MAKINUM (1));
5156 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5158 "Round the number @var{x} towards minus infinity.")
5159 #define FUNC_NAME s_scm_floor
5161 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5163 else if (SCM_REALP (x
))
5164 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5165 else if (SCM_FRACTIONP (x
))
5167 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5168 SCM_FRACTION_DENOMINATOR (x
));
5169 if (scm_is_false (scm_negative_p (x
)))
5171 /* For positive x, rounding towards zero is correct. */
5176 /* For negative x, we need to return q-1 unless x is an
5177 integer. But fractions are never integer, per our
5179 return scm_difference (q
, SCM_I_MAKINUM (1));
5183 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5187 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5189 "Round the number @var{x} towards infinity.")
5190 #define FUNC_NAME s_scm_ceiling
5192 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5194 else if (SCM_REALP (x
))
5195 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5196 else if (SCM_FRACTIONP (x
))
5198 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5199 SCM_FRACTION_DENOMINATOR (x
));
5200 if (scm_is_false (scm_positive_p (x
)))
5202 /* For negative x, rounding towards zero is correct. */
5207 /* For positive x, we need to return q+1 unless x is an
5208 integer. But fractions are never integer, per our
5210 return scm_sum (q
, SCM_I_MAKINUM (1));
5214 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5218 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5219 /* "Return the square root of the real number @var{x}."
5221 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5222 /* "Return the absolute value of the real number @var{x}."
5224 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5225 /* "Return the @var{x}th power of e."
5227 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5228 /* "Return the natural logarithm of the real number @var{x}."
5230 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5231 /* "Return the sine of the real number @var{x}."
5233 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5234 /* "Return the cosine of the real number @var{x}."
5236 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5237 /* "Return the tangent of the real number @var{x}."
5239 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5240 /* "Return the arc sine of the real number @var{x}."
5242 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5243 /* "Return the arc cosine of the real number @var{x}."
5245 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5246 /* "Return the arc tangent of the real number @var{x}."
5248 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5249 /* "Return the hyperbolic sine of the real number @var{x}."
5251 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5252 /* "Return the hyperbolic cosine of the real number @var{x}."
5254 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5255 /* "Return the hyperbolic tangent of the real number @var{x}."
5263 static void scm_two_doubles (SCM x
,
5265 const char *sstring
,
5269 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5271 if (SCM_I_INUMP (x
))
5272 xy
->x
= SCM_I_INUM (x
);
5273 else if (SCM_BIGP (x
))
5274 xy
->x
= scm_i_big2dbl (x
);
5275 else if (SCM_REALP (x
))
5276 xy
->x
= SCM_REAL_VALUE (x
);
5277 else if (SCM_FRACTIONP (x
))
5278 xy
->x
= scm_i_fraction2double (x
);
5280 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5282 if (SCM_I_INUMP (y
))
5283 xy
->y
= SCM_I_INUM (y
);
5284 else if (SCM_BIGP (y
))
5285 xy
->y
= scm_i_big2dbl (y
);
5286 else if (SCM_REALP (y
))
5287 xy
->y
= SCM_REAL_VALUE (y
);
5288 else if (SCM_FRACTIONP (y
))
5289 xy
->y
= scm_i_fraction2double (y
);
5291 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5295 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5297 "Return @var{x} raised to the power of @var{y}. This\n"
5298 "procedure does not accept complex arguments.")
5299 #define FUNC_NAME s_scm_sys_expt
5302 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5303 return scm_from_double (pow (xy
.x
, xy
.y
));
5308 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5310 "Return the arc tangent of the two arguments @var{x} and\n"
5311 "@var{y}. This is similar to calculating the arc tangent of\n"
5312 "@var{x} / @var{y}, except that the signs of both arguments\n"
5313 "are used to determine the quadrant of the result. This\n"
5314 "procedure does not accept complex arguments.")
5315 #define FUNC_NAME s_scm_sys_atan2
5318 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5319 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5324 scm_c_make_rectangular (double re
, double im
)
5327 return scm_from_double (re
);
5331 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5333 SCM_COMPLEX_REAL (z
) = re
;
5334 SCM_COMPLEX_IMAG (z
) = im
;
5339 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5340 (SCM real_part
, SCM imaginary_part
),
5341 "Return a complex number constructed of the given @var{real-part} "
5342 "and @var{imaginary-part} parts.")
5343 #define FUNC_NAME s_scm_make_rectangular
5346 scm_two_doubles (real_part
, imaginary_part
, FUNC_NAME
, &xy
);
5347 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5352 scm_c_make_polar (double mag
, double ang
)
5356 sincos (ang
, &s
, &c
);
5361 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5364 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5366 "Return the complex number @var{x} * e^(i * @var{y}).")
5367 #define FUNC_NAME s_scm_make_polar
5370 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5371 return scm_c_make_polar (xy
.x
, xy
.y
);
5376 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5377 /* "Return the real part of the number @var{z}."
5380 scm_real_part (SCM z
)
5382 if (SCM_I_INUMP (z
))
5384 else if (SCM_BIGP (z
))
5386 else if (SCM_REALP (z
))
5388 else if (SCM_COMPLEXP (z
))
5389 return scm_from_double (SCM_COMPLEX_REAL (z
));
5390 else if (SCM_FRACTIONP (z
))
5393 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5397 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5398 /* "Return the imaginary part of the number @var{z}."
5401 scm_imag_part (SCM z
)
5403 if (SCM_I_INUMP (z
))
5405 else if (SCM_BIGP (z
))
5407 else if (SCM_REALP (z
))
5409 else if (SCM_COMPLEXP (z
))
5410 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5411 else if (SCM_FRACTIONP (z
))
5414 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5417 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5418 /* "Return the numerator of the number @var{z}."
5421 scm_numerator (SCM z
)
5423 if (SCM_I_INUMP (z
))
5425 else if (SCM_BIGP (z
))
5427 else if (SCM_FRACTIONP (z
))
5428 return SCM_FRACTION_NUMERATOR (z
);
5429 else if (SCM_REALP (z
))
5430 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5432 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5436 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5437 /* "Return the denominator of the number @var{z}."
5440 scm_denominator (SCM z
)
5442 if (SCM_I_INUMP (z
))
5443 return SCM_I_MAKINUM (1);
5444 else if (SCM_BIGP (z
))
5445 return SCM_I_MAKINUM (1);
5446 else if (SCM_FRACTIONP (z
))
5447 return SCM_FRACTION_DENOMINATOR (z
);
5448 else if (SCM_REALP (z
))
5449 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5451 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5454 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5455 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5456 * "@code{abs} for real arguments, but also allows complex numbers."
5459 scm_magnitude (SCM z
)
5461 if (SCM_I_INUMP (z
))
5463 long int zz
= SCM_I_INUM (z
);
5466 else if (SCM_POSFIXABLE (-zz
))
5467 return SCM_I_MAKINUM (-zz
);
5469 return scm_i_long2big (-zz
);
5471 else if (SCM_BIGP (z
))
5473 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5474 scm_remember_upto_here_1 (z
);
5476 return scm_i_clonebig (z
, 0);
5480 else if (SCM_REALP (z
))
5481 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5482 else if (SCM_COMPLEXP (z
))
5483 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5484 else if (SCM_FRACTIONP (z
))
5486 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5488 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5489 SCM_FRACTION_DENOMINATOR (z
));
5492 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5496 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5497 /* "Return the angle of the complex number @var{z}."
5502 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5503 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5504 But if atan2 follows the floating point rounding mode, then the value
5505 is not a constant. Maybe it'd be close enough though. */
5506 if (SCM_I_INUMP (z
))
5508 if (SCM_I_INUM (z
) >= 0)
5511 return scm_from_double (atan2 (0.0, -1.0));
5513 else if (SCM_BIGP (z
))
5515 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5516 scm_remember_upto_here_1 (z
);
5518 return scm_from_double (atan2 (0.0, -1.0));
5522 else if (SCM_REALP (z
))
5524 if (SCM_REAL_VALUE (z
) >= 0)
5527 return scm_from_double (atan2 (0.0, -1.0));
5529 else if (SCM_COMPLEXP (z
))
5530 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5531 else if (SCM_FRACTIONP (z
))
5533 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5535 else return scm_from_double (atan2 (0.0, -1.0));
5538 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5542 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5543 /* Convert the number @var{x} to its inexact representation.\n"
5546 scm_exact_to_inexact (SCM z
)
5548 if (SCM_I_INUMP (z
))
5549 return scm_from_double ((double) SCM_I_INUM (z
));
5550 else if (SCM_BIGP (z
))
5551 return scm_from_double (scm_i_big2dbl (z
));
5552 else if (SCM_FRACTIONP (z
))
5553 return scm_from_double (scm_i_fraction2double (z
));
5554 else if (SCM_INEXACTP (z
))
5557 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5561 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5563 "Return an exact number that is numerically closest to @var{z}.")
5564 #define FUNC_NAME s_scm_inexact_to_exact
5566 if (SCM_I_INUMP (z
))
5568 else if (SCM_BIGP (z
))
5570 else if (SCM_REALP (z
))
5572 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5573 SCM_OUT_OF_RANGE (1, z
);
5580 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5581 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5582 scm_i_mpz2num (mpq_denref (frac
)));
5584 /* When scm_i_make_ratio throws, we leak the memory allocated
5591 else if (SCM_FRACTIONP (z
))
5594 SCM_WRONG_TYPE_ARG (1, z
);
5598 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5600 "Returns the @emph{simplest} rational number differing\n"
5601 "from @var{x} by no more than @var{eps}.\n"
5603 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5604 "exact result when both its arguments are exact. Thus, you might need\n"
5605 "to use @code{inexact->exact} on the arguments.\n"
5608 "(rationalize (inexact->exact 1.2) 1/100)\n"
5611 #define FUNC_NAME s_scm_rationalize
5613 if (SCM_I_INUMP (x
))
5615 else if (SCM_BIGP (x
))
5617 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5619 /* Use continued fractions to find closest ratio. All
5620 arithmetic is done with exact numbers.
5623 SCM ex
= scm_inexact_to_exact (x
);
5624 SCM int_part
= scm_floor (ex
);
5625 SCM tt
= SCM_I_MAKINUM (1);
5626 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5627 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5631 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5634 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5635 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5637 /* We stop after a million iterations just to be absolutely sure
5638 that we don't go into an infinite loop. The process normally
5639 converges after less than a dozen iterations.
5642 eps
= scm_abs (eps
);
5643 while (++i
< 1000000)
5645 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5646 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5647 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5649 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5650 eps
))) /* abs(x-a/b) <= eps */
5652 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5653 if (scm_is_false (scm_exact_p (x
))
5654 || scm_is_false (scm_exact_p (eps
)))
5655 return scm_exact_to_inexact (res
);
5659 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5661 tt
= scm_floor (rx
); /* tt = floor (rx) */
5667 scm_num_overflow (s_scm_rationalize
);
5670 SCM_WRONG_TYPE_ARG (1, x
);
5674 /* conversion functions */
5677 scm_is_integer (SCM val
)
5679 return scm_is_true (scm_integer_p (val
));
5683 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5685 if (SCM_I_INUMP (val
))
5687 scm_t_signed_bits n
= SCM_I_INUM (val
);
5688 return n
>= min
&& n
<= max
;
5690 else if (SCM_BIGP (val
))
5692 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5694 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5696 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5698 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5699 return n
>= min
&& n
<= max
;
5709 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5710 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5713 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5714 SCM_I_BIG_MPZ (val
));
5716 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5728 return n
>= min
&& n
<= max
;
5736 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5738 if (SCM_I_INUMP (val
))
5740 scm_t_signed_bits n
= SCM_I_INUM (val
);
5741 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5743 else if (SCM_BIGP (val
))
5745 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5747 else if (max
<= ULONG_MAX
)
5749 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5751 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5752 return n
>= min
&& n
<= max
;
5762 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5765 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5766 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5769 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5770 SCM_I_BIG_MPZ (val
));
5772 return n
>= min
&& n
<= max
;
5780 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5782 scm_error (scm_out_of_range_key
,
5784 "Value out of range ~S to ~S: ~S",
5785 scm_list_3 (min
, max
, bad_val
),
5786 scm_list_1 (bad_val
));
5789 #define TYPE scm_t_intmax
5790 #define TYPE_MIN min
5791 #define TYPE_MAX max
5792 #define SIZEOF_TYPE 0
5793 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5794 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5795 #include "libguile/conv-integer.i.c"
5797 #define TYPE scm_t_uintmax
5798 #define TYPE_MIN min
5799 #define TYPE_MAX max
5800 #define SIZEOF_TYPE 0
5801 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5802 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5803 #include "libguile/conv-uinteger.i.c"
5805 #define TYPE scm_t_int8
5806 #define TYPE_MIN SCM_T_INT8_MIN
5807 #define TYPE_MAX SCM_T_INT8_MAX
5808 #define SIZEOF_TYPE 1
5809 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5810 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5811 #include "libguile/conv-integer.i.c"
5813 #define TYPE scm_t_uint8
5815 #define TYPE_MAX SCM_T_UINT8_MAX
5816 #define SIZEOF_TYPE 1
5817 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5818 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5819 #include "libguile/conv-uinteger.i.c"
5821 #define TYPE scm_t_int16
5822 #define TYPE_MIN SCM_T_INT16_MIN
5823 #define TYPE_MAX SCM_T_INT16_MAX
5824 #define SIZEOF_TYPE 2
5825 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5826 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5827 #include "libguile/conv-integer.i.c"
5829 #define TYPE scm_t_uint16
5831 #define TYPE_MAX SCM_T_UINT16_MAX
5832 #define SIZEOF_TYPE 2
5833 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5834 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5835 #include "libguile/conv-uinteger.i.c"
5837 #define TYPE scm_t_int32
5838 #define TYPE_MIN SCM_T_INT32_MIN
5839 #define TYPE_MAX SCM_T_INT32_MAX
5840 #define SIZEOF_TYPE 4
5841 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5842 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5843 #include "libguile/conv-integer.i.c"
5845 #define TYPE scm_t_uint32
5847 #define TYPE_MAX SCM_T_UINT32_MAX
5848 #define SIZEOF_TYPE 4
5849 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5850 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5851 #include "libguile/conv-uinteger.i.c"
5853 #if SCM_HAVE_T_INT64
5855 #define TYPE scm_t_int64
5856 #define TYPE_MIN SCM_T_INT64_MIN
5857 #define TYPE_MAX SCM_T_INT64_MAX
5858 #define SIZEOF_TYPE 8
5859 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5860 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5861 #include "libguile/conv-integer.i.c"
5863 #define TYPE scm_t_uint64
5865 #define TYPE_MAX SCM_T_UINT64_MAX
5866 #define SIZEOF_TYPE 8
5867 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5868 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5869 #include "libguile/conv-uinteger.i.c"
5874 scm_to_mpz (SCM val
, mpz_t rop
)
5876 if (SCM_I_INUMP (val
))
5877 mpz_set_si (rop
, SCM_I_INUM (val
));
5878 else if (SCM_BIGP (val
))
5879 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5881 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5885 scm_from_mpz (mpz_t val
)
5887 return scm_i_mpz2num (val
);
5891 scm_is_real (SCM val
)
5893 return scm_is_true (scm_real_p (val
));
5897 scm_is_rational (SCM val
)
5899 return scm_is_true (scm_rational_p (val
));
5903 scm_to_double (SCM val
)
5905 if (SCM_I_INUMP (val
))
5906 return SCM_I_INUM (val
);
5907 else if (SCM_BIGP (val
))
5908 return scm_i_big2dbl (val
);
5909 else if (SCM_FRACTIONP (val
))
5910 return scm_i_fraction2double (val
);
5911 else if (SCM_REALP (val
))
5912 return SCM_REAL_VALUE (val
);
5914 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5918 scm_from_double (double val
)
5920 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5921 SCM_REAL_VALUE (z
) = val
;
5925 #if SCM_ENABLE_DISCOURAGED == 1
5928 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5932 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5936 scm_out_of_range (NULL
, num
);
5939 return scm_to_double (num
);
5943 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5947 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5951 scm_out_of_range (NULL
, num
);
5954 return scm_to_double (num
);
5960 scm_is_complex (SCM val
)
5962 return scm_is_true (scm_complex_p (val
));
5966 scm_c_real_part (SCM z
)
5968 if (SCM_COMPLEXP (z
))
5969 return SCM_COMPLEX_REAL (z
);
5972 /* Use the scm_real_part to get proper error checking and
5975 return scm_to_double (scm_real_part (z
));
5980 scm_c_imag_part (SCM z
)
5982 if (SCM_COMPLEXP (z
))
5983 return SCM_COMPLEX_IMAG (z
);
5986 /* Use the scm_imag_part to get proper error checking and
5987 dispatching. The result will almost always be 0.0, but not
5990 return scm_to_double (scm_imag_part (z
));
5995 scm_c_magnitude (SCM z
)
5997 return scm_to_double (scm_magnitude (z
));
6003 return scm_to_double (scm_angle (z
));
6007 scm_is_number (SCM z
)
6009 return scm_is_true (scm_number_p (z
));
6013 /* In the following functions we dispatch to the real-arg funcs like log()
6014 when we know the arg is real, instead of just handing everything to
6015 clog() for instance. This is in case clog() doesn't optimize for a
6016 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6017 well use it to go straight to the applicable C func. */
6019 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6021 "Return the natural logarithm of @var{z}.")
6022 #define FUNC_NAME s_scm_log
6024 if (SCM_COMPLEXP (z
))
6026 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6027 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6029 double re
= SCM_COMPLEX_REAL (z
);
6030 double im
= SCM_COMPLEX_IMAG (z
);
6031 return scm_c_make_rectangular (log (hypot (re
, im
)),
6037 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6038 although the value itself overflows. */
6039 double re
= scm_to_double (z
);
6040 double l
= log (fabs (re
));
6042 return scm_from_double (l
);
6044 return scm_c_make_rectangular (l
, M_PI
);
6050 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6052 "Return the base 10 logarithm of @var{z}.")
6053 #define FUNC_NAME s_scm_log10
6055 if (SCM_COMPLEXP (z
))
6057 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6058 clog() and a multiply by M_LOG10E, rather than the fallback
6059 log10+hypot+atan2.) */
6060 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6061 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6063 double re
= SCM_COMPLEX_REAL (z
);
6064 double im
= SCM_COMPLEX_IMAG (z
);
6065 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6066 M_LOG10E
* atan2 (im
, re
));
6071 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6072 although the value itself overflows. */
6073 double re
= scm_to_double (z
);
6074 double l
= log10 (fabs (re
));
6076 return scm_from_double (l
);
6078 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6084 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6086 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6087 "base of natural logarithms (2.71828@dots{}).")
6088 #define FUNC_NAME s_scm_exp
6090 if (SCM_COMPLEXP (z
))
6092 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6093 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6095 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6096 SCM_COMPLEX_IMAG (z
));
6101 /* When z is a negative bignum the conversion to double overflows,
6102 giving -infinity, but that's ok, the exp is still 0.0. */
6103 return scm_from_double (exp (scm_to_double (z
)));
6109 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6111 "Return the square root of @var{z}. Of the two possible roots\n"
6112 "(positive and negative), the one with the a positive real part\n"
6113 "is returned, or if that's zero then a positive imaginary part.\n"
6117 "(sqrt 9.0) @result{} 3.0\n"
6118 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6119 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6120 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6122 #define FUNC_NAME s_scm_sqrt
6124 if (SCM_COMPLEXP (x
))
6126 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6127 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6129 double re
= SCM_COMPLEX_REAL (x
);
6130 double im
= SCM_COMPLEX_IMAG (x
);
6131 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6132 0.5 * atan2 (im
, re
));
6137 double xx
= scm_to_double (x
);
6139 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6141 return scm_from_double (sqrt (xx
));
6153 mpz_init_set_si (z_negative_one
, -1);
6155 /* It may be possible to tune the performance of some algorithms by using
6156 * the following constants to avoid the creation of bignums. Please, before
6157 * using these values, remember the two rules of program optimization:
6158 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6159 scm_c_define ("most-positive-fixnum",
6160 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6161 scm_c_define ("most-negative-fixnum",
6162 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6164 scm_add_feature ("complex");
6165 scm_add_feature ("inexact");
6166 scm_flo0
= scm_from_double (0.0);
6168 /* determine floating point precision */
6169 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6171 init_dblprec(&scm_dblprec
[i
-2],i
);
6172 init_fx_radix(fx_per_radix
[i
-2],i
);
6175 /* hard code precision for base 10 if the preprocessor tells us to... */
6176 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6179 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6180 SCM_I_MAKINUM (2)));
6181 #include "libguile/numbers.x"