1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public License
9 * as published by the Free Software Foundation; either version 3 of
10 * the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
24 /* General assumptions:
25 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
26 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
27 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
28 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
29 * All objects satisfying SCM_FRACTIONP are never an integer.
34 - see if special casing bignums and reals in integer-exponent when
35 possible (to use mpz_pow and mpf_pow_ui) is faster.
37 - look in to better short-circuiting of common cases in
38 integer-expt and elsewhere.
40 - see if direct mpz operations can help in ash and elsewhere.
57 #include "libguile/_scm.h"
58 #include "libguile/feature.h"
59 #include "libguile/ports.h"
60 #include "libguile/root.h"
61 #include "libguile/smob.h"
62 #include "libguile/strings.h"
64 #include "libguile/validate.h"
65 #include "libguile/numbers.h"
66 #include "libguile/deprecation.h"
68 #include "libguile/eq.h"
70 #include "libguile/discouraged.h"
72 /* values per glibc, if not already defined */
74 #define M_LOG10E 0.43429448190325182765
77 #define M_PI 3.14159265358979323846
83 Wonder if this might be faster for some of our code? A switch on
84 the numtag would jump directly to the right case, and the
85 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
87 #define SCM_I_NUMTAG_NOTNUM 0
88 #define SCM_I_NUMTAG_INUM 1
89 #define SCM_I_NUMTAG_BIG scm_tc16_big
90 #define SCM_I_NUMTAG_REAL scm_tc16_real
91 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
92 #define SCM_I_NUMTAG(x) \
93 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
94 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
95 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
96 : SCM_I_NUMTAG_NOTNUM)))
98 /* the macro above will not work as is with fractions */
101 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
103 /* FLOBUFLEN is the maximum number of characters neccessary for the
104 * printed or scm_string representation of an inexact number.
106 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
109 #if ! defined (HAVE_ISNAN)
114 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
117 #if ! defined (HAVE_ISINF)
122 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
129 #if !defined (HAVE_ASINH)
130 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
132 #if !defined (HAVE_ACOSH)
133 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
135 #if !defined (HAVE_ATANH)
136 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
139 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
140 an explicit check. In some future gmp (don't know what version number),
141 mpz_cmp_d is supposed to do this itself. */
143 #define xmpz_cmp_d(z, d) \
144 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
146 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
149 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
150 isinf. It does have finite and isnan though, hence the use of those.
151 fpclass would be a possibility on that system too. */
155 #if defined (HAVE_ISINF)
157 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
158 return (! (finite (x
) || isnan (x
)));
167 #if defined (HAVE_ISNAN)
174 #if defined (GUILE_I)
175 #if HAVE_COMPLEX_DOUBLE
177 /* For an SCM object Z which is a complex number (ie. satisfies
178 SCM_COMPLEXP), return its value as a C level "complex double". */
179 #define SCM_COMPLEX_VALUE(z) \
180 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
182 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
184 /* Convert a C "complex double" to an SCM value. */
186 scm_from_complex_double (complex double z
)
188 return scm_c_make_rectangular (creal (z
), cimag (z
));
191 #endif /* HAVE_COMPLEX_DOUBLE */
196 static mpz_t z_negative_one
;
203 /* Return a newly created bignum. */
204 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
205 mpz_init (SCM_I_BIG_MPZ (z
));
210 scm_i_long2big (long x
)
212 /* Return a newly created bignum initialized to X. */
213 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
214 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
219 scm_i_ulong2big (unsigned long x
)
221 /* Return a newly created bignum initialized to X. */
222 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
223 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
228 scm_i_clonebig (SCM src_big
, int same_sign_p
)
230 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
231 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
232 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
234 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
239 scm_i_bigcmp (SCM x
, SCM y
)
241 /* Return neg if x < y, pos if x > y, and 0 if x == y */
242 /* presume we already know x and y are bignums */
243 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
244 scm_remember_upto_here_2 (x
, y
);
249 scm_i_dbl2big (double d
)
251 /* results are only defined if d is an integer */
252 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
253 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
257 /* Convert a integer in double representation to a SCM number. */
260 scm_i_dbl2num (double u
)
262 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
263 powers of 2, so there's no rounding when making "double" values
264 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
265 get rounded on a 64-bit machine, hence the "+1".
267 The use of floor() to force to an integer value ensures we get a
268 "numerically closest" value without depending on how a
269 double->long cast or how mpz_set_d will round. For reference,
270 double->long probably follows the hardware rounding mode,
271 mpz_set_d truncates towards zero. */
273 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
274 representable as a double? */
276 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
277 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
278 return SCM_I_MAKINUM ((long) u
);
280 return scm_i_dbl2big (u
);
283 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
284 with R5RS exact->inexact.
286 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
287 (ie. truncate towards zero), then adjust to get the closest double by
288 examining the next lower bit and adding 1 (to the absolute value) if
291 Bignums exactly half way between representable doubles are rounded to the
292 next higher absolute value (ie. away from zero). This seems like an
293 adequate interpretation of R5RS "numerically closest", and it's easier
294 and faster than a full "nearest-even" style.
296 The bit test must be done on the absolute value of the mpz_t, which means
297 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
298 negatives as twos complement.
300 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
301 following the hardware rounding mode, but applied to the absolute value
302 of the mpz_t operand. This is not what we want so we put the high
303 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
304 mpz_get_d is supposed to always truncate towards zero.
306 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
307 is a slowdown. It'd be faster to pick out the relevant high bits with
308 mpz_getlimbn if we could be bothered coding that, and if the new
309 truncating gmp doesn't come out. */
312 scm_i_big2dbl (SCM b
)
317 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
321 /* Current GMP, eg. 4.1.3, force truncation towards zero */
323 if (bits
> DBL_MANT_DIG
)
325 size_t shift
= bits
- DBL_MANT_DIG
;
326 mpz_init2 (tmp
, DBL_MANT_DIG
);
327 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
328 result
= ldexp (mpz_get_d (tmp
), shift
);
333 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
338 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
341 if (bits
> DBL_MANT_DIG
)
343 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
344 /* test bit number "pos" in absolute value */
345 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
346 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
348 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
352 scm_remember_upto_here_1 (b
);
357 scm_i_normbig (SCM b
)
359 /* convert a big back to a fixnum if it'll fit */
360 /* presume b is a bignum */
361 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
363 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
364 if (SCM_FIXABLE (val
))
365 b
= SCM_I_MAKINUM (val
);
370 static SCM_C_INLINE_KEYWORD SCM
371 scm_i_mpz2num (mpz_t b
)
373 /* convert a mpz number to a SCM number. */
374 if (mpz_fits_slong_p (b
))
376 long val
= mpz_get_si (b
);
377 if (SCM_FIXABLE (val
))
378 return SCM_I_MAKINUM (val
);
382 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
383 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
388 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
389 static SCM
scm_divide2real (SCM x
, SCM y
);
392 scm_i_make_ratio (SCM numerator
, SCM denominator
)
393 #define FUNC_NAME "make-ratio"
395 /* First make sure the arguments are proper.
397 if (SCM_I_INUMP (denominator
))
399 if (scm_is_eq (denominator
, SCM_INUM0
))
400 scm_num_overflow ("make-ratio");
401 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
406 if (!(SCM_BIGP(denominator
)))
407 SCM_WRONG_TYPE_ARG (2, denominator
);
409 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
410 SCM_WRONG_TYPE_ARG (1, numerator
);
412 /* Then flip signs so that the denominator is positive.
414 if (scm_is_true (scm_negative_p (denominator
)))
416 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
417 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
420 /* Now consider for each of the four fixnum/bignum combinations
421 whether the rational number is really an integer.
423 if (SCM_I_INUMP (numerator
))
425 long x
= SCM_I_INUM (numerator
);
426 if (scm_is_eq (numerator
, SCM_INUM0
))
428 if (SCM_I_INUMP (denominator
))
431 y
= SCM_I_INUM (denominator
);
433 return SCM_I_MAKINUM(1);
435 return SCM_I_MAKINUM (x
/ y
);
439 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
440 of that value for the denominator, as a bignum. Apart from
441 that case, abs(bignum) > abs(inum) so inum/bignum is not an
443 if (x
== SCM_MOST_NEGATIVE_FIXNUM
444 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
445 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
446 return SCM_I_MAKINUM(-1);
449 else if (SCM_BIGP (numerator
))
451 if (SCM_I_INUMP (denominator
))
453 long yy
= SCM_I_INUM (denominator
);
454 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
455 return scm_divide (numerator
, denominator
);
459 if (scm_is_eq (numerator
, denominator
))
460 return SCM_I_MAKINUM(1);
461 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
462 SCM_I_BIG_MPZ (denominator
)))
463 return scm_divide(numerator
, denominator
);
467 /* No, it's a proper fraction.
470 SCM divisor
= scm_gcd (numerator
, denominator
);
471 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
473 numerator
= scm_divide (numerator
, divisor
);
474 denominator
= scm_divide (denominator
, divisor
);
477 return scm_double_cell (scm_tc16_fraction
,
478 SCM_UNPACK (numerator
),
479 SCM_UNPACK (denominator
), 0);
485 scm_i_fraction2double (SCM z
)
487 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
488 SCM_FRACTION_DENOMINATOR (z
)));
491 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
493 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
495 #define FUNC_NAME s_scm_exact_p
501 if (SCM_FRACTIONP (x
))
505 SCM_WRONG_TYPE_ARG (1, x
);
510 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
512 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
514 #define FUNC_NAME s_scm_odd_p
518 long val
= SCM_I_INUM (n
);
519 return scm_from_bool ((val
& 1L) != 0);
521 else if (SCM_BIGP (n
))
523 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
524 scm_remember_upto_here_1 (n
);
525 return scm_from_bool (odd_p
);
527 else if (scm_is_true (scm_inf_p (n
)))
529 else if (SCM_REALP (n
))
531 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
537 SCM_WRONG_TYPE_ARG (1, n
);
540 SCM_WRONG_TYPE_ARG (1, n
);
545 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
547 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
549 #define FUNC_NAME s_scm_even_p
553 long val
= SCM_I_INUM (n
);
554 return scm_from_bool ((val
& 1L) == 0);
556 else if (SCM_BIGP (n
))
558 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
559 scm_remember_upto_here_1 (n
);
560 return scm_from_bool (even_p
);
562 else if (scm_is_true (scm_inf_p (n
)))
564 else if (SCM_REALP (n
))
566 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
572 SCM_WRONG_TYPE_ARG (1, n
);
575 SCM_WRONG_TYPE_ARG (1, n
);
579 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
581 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
582 "or @samp{-inf.0}, @code{#f} otherwise.")
583 #define FUNC_NAME s_scm_inf_p
586 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
587 else if (SCM_COMPLEXP (x
))
588 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
589 || xisinf (SCM_COMPLEX_IMAG (x
)));
595 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
597 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
599 #define FUNC_NAME s_scm_nan_p
602 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
603 else if (SCM_COMPLEXP (n
))
604 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
605 || xisnan (SCM_COMPLEX_IMAG (n
)));
611 /* Guile's idea of infinity. */
612 static double guile_Inf
;
614 /* Guile's idea of not a number. */
615 static double guile_NaN
;
618 guile_ieee_init (void)
620 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
622 /* Some version of gcc on some old version of Linux used to crash when
623 trying to make Inf and NaN. */
626 /* C99 INFINITY, when available.
627 FIXME: The standard allows for INFINITY to be something that overflows
628 at compile time. We ought to have a configure test to check for that
629 before trying to use it. (But in practice we believe this is not a
630 problem on any system guile is likely to target.) */
631 guile_Inf
= INFINITY
;
632 #elif defined HAVE_DINFINITY
634 extern unsigned int DINFINITY
[2];
635 guile_Inf
= (*((double *) (DINFINITY
)));
642 if (guile_Inf
== tmp
)
650 #if defined (HAVE_ISNAN)
653 /* C99 NAN, when available */
655 #elif defined HAVE_DQNAN
658 extern unsigned int DQNAN
[2];
659 guile_NaN
= (*((double *)(DQNAN
)));
662 guile_NaN
= guile_Inf
/ guile_Inf
;
668 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
671 #define FUNC_NAME s_scm_inf
673 static int initialized
= 0;
679 return scm_from_double (guile_Inf
);
683 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
686 #define FUNC_NAME s_scm_nan
688 static int initialized
= 0;
694 return scm_from_double (guile_NaN
);
699 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
701 "Return the absolute value of @var{x}.")
706 long int xx
= SCM_I_INUM (x
);
709 else if (SCM_POSFIXABLE (-xx
))
710 return SCM_I_MAKINUM (-xx
);
712 return scm_i_long2big (-xx
);
714 else if (SCM_BIGP (x
))
716 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
718 return scm_i_clonebig (x
, 0);
722 else if (SCM_REALP (x
))
724 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
725 double xx
= SCM_REAL_VALUE (x
);
727 return scm_from_double (-xx
);
731 else if (SCM_FRACTIONP (x
))
733 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
735 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
736 SCM_FRACTION_DENOMINATOR (x
));
739 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
744 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
745 /* "Return the quotient of the numbers @var{x} and @var{y}."
748 scm_quotient (SCM x
, SCM y
)
752 long xx
= SCM_I_INUM (x
);
755 long yy
= SCM_I_INUM (y
);
757 scm_num_overflow (s_quotient
);
762 return SCM_I_MAKINUM (z
);
764 return scm_i_long2big (z
);
767 else if (SCM_BIGP (y
))
769 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
770 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
771 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
773 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
774 scm_remember_upto_here_1 (y
);
775 return SCM_I_MAKINUM (-1);
778 return SCM_I_MAKINUM (0);
781 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
783 else if (SCM_BIGP (x
))
787 long yy
= SCM_I_INUM (y
);
789 scm_num_overflow (s_quotient
);
794 SCM result
= scm_i_mkbig ();
797 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
800 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
803 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
804 scm_remember_upto_here_1 (x
);
805 return scm_i_normbig (result
);
808 else if (SCM_BIGP (y
))
810 SCM result
= scm_i_mkbig ();
811 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
814 scm_remember_upto_here_2 (x
, y
);
815 return scm_i_normbig (result
);
818 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
821 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
824 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
825 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
827 * "(remainder 13 4) @result{} 1\n"
828 * "(remainder -13 4) @result{} -1\n"
832 scm_remainder (SCM x
, SCM y
)
838 long yy
= SCM_I_INUM (y
);
840 scm_num_overflow (s_remainder
);
843 long z
= SCM_I_INUM (x
) % yy
;
844 return SCM_I_MAKINUM (z
);
847 else if (SCM_BIGP (y
))
849 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
850 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
851 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
853 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
854 scm_remember_upto_here_1 (y
);
855 return SCM_I_MAKINUM (0);
861 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
863 else if (SCM_BIGP (x
))
867 long yy
= SCM_I_INUM (y
);
869 scm_num_overflow (s_remainder
);
872 SCM result
= scm_i_mkbig ();
875 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
876 scm_remember_upto_here_1 (x
);
877 return scm_i_normbig (result
);
880 else if (SCM_BIGP (y
))
882 SCM result
= scm_i_mkbig ();
883 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
886 scm_remember_upto_here_2 (x
, y
);
887 return scm_i_normbig (result
);
890 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
893 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
897 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
898 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
900 * "(modulo 13 4) @result{} 1\n"
901 * "(modulo -13 4) @result{} 3\n"
905 scm_modulo (SCM x
, SCM y
)
909 long xx
= SCM_I_INUM (x
);
912 long yy
= SCM_I_INUM (y
);
914 scm_num_overflow (s_modulo
);
917 /* C99 specifies that "%" is the remainder corresponding to a
918 quotient rounded towards zero, and that's also traditional
919 for machine division, so z here should be well defined. */
937 return SCM_I_MAKINUM (result
);
940 else if (SCM_BIGP (y
))
942 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
949 SCM pos_y
= scm_i_clonebig (y
, 0);
950 /* do this after the last scm_op */
951 mpz_init_set_si (z_x
, xx
);
952 result
= pos_y
; /* re-use this bignum */
953 mpz_mod (SCM_I_BIG_MPZ (result
),
955 SCM_I_BIG_MPZ (pos_y
));
956 scm_remember_upto_here_1 (pos_y
);
960 result
= scm_i_mkbig ();
961 /* do this after the last scm_op */
962 mpz_init_set_si (z_x
, xx
);
963 mpz_mod (SCM_I_BIG_MPZ (result
),
966 scm_remember_upto_here_1 (y
);
969 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
970 mpz_add (SCM_I_BIG_MPZ (result
),
972 SCM_I_BIG_MPZ (result
));
973 scm_remember_upto_here_1 (y
);
974 /* and do this before the next one */
976 return scm_i_normbig (result
);
980 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
982 else if (SCM_BIGP (x
))
986 long yy
= SCM_I_INUM (y
);
988 scm_num_overflow (s_modulo
);
991 SCM result
= scm_i_mkbig ();
992 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
994 (yy
< 0) ? - yy
: yy
);
995 scm_remember_upto_here_1 (x
);
996 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
997 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
998 SCM_I_BIG_MPZ (result
),
1000 return scm_i_normbig (result
);
1003 else if (SCM_BIGP (y
))
1006 SCM result
= scm_i_mkbig ();
1007 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1008 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1009 mpz_mod (SCM_I_BIG_MPZ (result
),
1011 SCM_I_BIG_MPZ (pos_y
));
1013 scm_remember_upto_here_1 (x
);
1014 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1015 mpz_add (SCM_I_BIG_MPZ (result
),
1017 SCM_I_BIG_MPZ (result
));
1018 scm_remember_upto_here_2 (y
, pos_y
);
1019 return scm_i_normbig (result
);
1023 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1026 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1029 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1030 (SCM x
, SCM y
, SCM rest
),
1031 "Return the greatest common divisor of all parameter values.\n"
1032 "If called without arguments, 0 is returned.")
1033 #define FUNC_NAME s_scm_i_gcd
1035 while (!scm_is_null (rest
))
1036 { x
= scm_gcd (x
, y
);
1038 rest
= scm_cdr (rest
);
1040 return scm_gcd (x
, y
);
1044 #define s_gcd s_scm_i_gcd
1045 #define g_gcd g_scm_i_gcd
1048 scm_gcd (SCM x
, SCM y
)
1051 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1053 if (SCM_I_INUMP (x
))
1055 if (SCM_I_INUMP (y
))
1057 long xx
= SCM_I_INUM (x
);
1058 long yy
= SCM_I_INUM (y
);
1059 long u
= xx
< 0 ? -xx
: xx
;
1060 long v
= yy
< 0 ? -yy
: yy
;
1070 /* Determine a common factor 2^k */
1071 while (!(1 & (u
| v
)))
1077 /* Now, any factor 2^n can be eliminated */
1097 return (SCM_POSFIXABLE (result
)
1098 ? SCM_I_MAKINUM (result
)
1099 : scm_i_long2big (result
));
1101 else if (SCM_BIGP (y
))
1107 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1109 else if (SCM_BIGP (x
))
1111 if (SCM_I_INUMP (y
))
1113 unsigned long result
;
1116 yy
= SCM_I_INUM (y
);
1121 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1122 scm_remember_upto_here_1 (x
);
1123 return (SCM_POSFIXABLE (result
)
1124 ? SCM_I_MAKINUM (result
)
1125 : scm_from_ulong (result
));
1127 else if (SCM_BIGP (y
))
1129 SCM result
= scm_i_mkbig ();
1130 mpz_gcd (SCM_I_BIG_MPZ (result
),
1133 scm_remember_upto_here_2 (x
, y
);
1134 return scm_i_normbig (result
);
1137 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1140 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1143 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1144 (SCM x
, SCM y
, SCM rest
),
1145 "Return the least common multiple of the arguments.\n"
1146 "If called without arguments, 1 is returned.")
1147 #define FUNC_NAME s_scm_i_lcm
1149 while (!scm_is_null (rest
))
1150 { x
= scm_lcm (x
, y
);
1152 rest
= scm_cdr (rest
);
1154 return scm_lcm (x
, y
);
1158 #define s_lcm s_scm_i_lcm
1159 #define g_lcm g_scm_i_lcm
1162 scm_lcm (SCM n1
, SCM n2
)
1164 if (SCM_UNBNDP (n2
))
1166 if (SCM_UNBNDP (n1
))
1167 return SCM_I_MAKINUM (1L);
1168 n2
= SCM_I_MAKINUM (1L);
1171 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1172 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1173 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1174 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1176 if (SCM_I_INUMP (n1
))
1178 if (SCM_I_INUMP (n2
))
1180 SCM d
= scm_gcd (n1
, n2
);
1181 if (scm_is_eq (d
, SCM_INUM0
))
1184 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1188 /* inum n1, big n2 */
1191 SCM result
= scm_i_mkbig ();
1192 long nn1
= SCM_I_INUM (n1
);
1193 if (nn1
== 0) return SCM_INUM0
;
1194 if (nn1
< 0) nn1
= - nn1
;
1195 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1196 scm_remember_upto_here_1 (n2
);
1204 if (SCM_I_INUMP (n2
))
1211 SCM result
= scm_i_mkbig ();
1212 mpz_lcm(SCM_I_BIG_MPZ (result
),
1214 SCM_I_BIG_MPZ (n2
));
1215 scm_remember_upto_here_2(n1
, n2
);
1216 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1222 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1227 + + + x (map digit:logand X Y)
1228 + - + x (map digit:logand X (lognot (+ -1 Y)))
1229 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1230 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1235 + + + (map digit:logior X Y)
1236 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1237 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1238 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1243 + + + (map digit:logxor X Y)
1244 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1245 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1246 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1251 + + (any digit:logand X Y)
1252 + - (any digit:logand X (lognot (+ -1 Y)))
1253 - + (any digit:logand (lognot (+ -1 X)) Y)
1258 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1259 (SCM x
, SCM y
, SCM rest
),
1260 "Return the bitwise AND of the integer arguments.\n\n"
1262 "(logand) @result{} -1\n"
1263 "(logand 7) @result{} 7\n"
1264 "(logand #b111 #b011 #b001) @result{} 1\n"
1266 #define FUNC_NAME s_scm_i_logand
1268 while (!scm_is_null (rest
))
1269 { x
= scm_logand (x
, y
);
1271 rest
= scm_cdr (rest
);
1273 return scm_logand (x
, y
);
1277 #define s_scm_logand s_scm_i_logand
1279 SCM
scm_logand (SCM n1
, SCM n2
)
1280 #define FUNC_NAME s_scm_logand
1284 if (SCM_UNBNDP (n2
))
1286 if (SCM_UNBNDP (n1
))
1287 return SCM_I_MAKINUM (-1);
1288 else if (!SCM_NUMBERP (n1
))
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1290 else if (SCM_NUMBERP (n1
))
1293 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1296 if (SCM_I_INUMP (n1
))
1298 nn1
= SCM_I_INUM (n1
);
1299 if (SCM_I_INUMP (n2
))
1301 long nn2
= SCM_I_INUM (n2
);
1302 return SCM_I_MAKINUM (nn1
& nn2
);
1304 else if SCM_BIGP (n2
)
1310 SCM result_z
= scm_i_mkbig ();
1312 mpz_init_set_si (nn1_z
, nn1
);
1313 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1314 scm_remember_upto_here_1 (n2
);
1316 return scm_i_normbig (result_z
);
1320 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1322 else if (SCM_BIGP (n1
))
1324 if (SCM_I_INUMP (n2
))
1327 nn1
= SCM_I_INUM (n1
);
1330 else if (SCM_BIGP (n2
))
1332 SCM result_z
= scm_i_mkbig ();
1333 mpz_and (SCM_I_BIG_MPZ (result_z
),
1335 SCM_I_BIG_MPZ (n2
));
1336 scm_remember_upto_here_2 (n1
, n2
);
1337 return scm_i_normbig (result_z
);
1340 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1343 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1348 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1349 (SCM x
, SCM y
, SCM rest
),
1350 "Return the bitwise OR of the integer arguments.\n\n"
1352 "(logior) @result{} 0\n"
1353 "(logior 7) @result{} 7\n"
1354 "(logior #b000 #b001 #b011) @result{} 3\n"
1356 #define FUNC_NAME s_scm_i_logior
1358 while (!scm_is_null (rest
))
1359 { x
= scm_logior (x
, y
);
1361 rest
= scm_cdr (rest
);
1363 return scm_logior (x
, y
);
1367 #define s_scm_logior s_scm_i_logior
1369 SCM
scm_logior (SCM n1
, SCM n2
)
1370 #define FUNC_NAME s_scm_logior
1374 if (SCM_UNBNDP (n2
))
1376 if (SCM_UNBNDP (n1
))
1378 else if (SCM_NUMBERP (n1
))
1381 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1384 if (SCM_I_INUMP (n1
))
1386 nn1
= SCM_I_INUM (n1
);
1387 if (SCM_I_INUMP (n2
))
1389 long nn2
= SCM_I_INUM (n2
);
1390 return SCM_I_MAKINUM (nn1
| nn2
);
1392 else if (SCM_BIGP (n2
))
1398 SCM result_z
= scm_i_mkbig ();
1400 mpz_init_set_si (nn1_z
, nn1
);
1401 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1402 scm_remember_upto_here_1 (n2
);
1404 return scm_i_normbig (result_z
);
1408 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1410 else if (SCM_BIGP (n1
))
1412 if (SCM_I_INUMP (n2
))
1415 nn1
= SCM_I_INUM (n1
);
1418 else if (SCM_BIGP (n2
))
1420 SCM result_z
= scm_i_mkbig ();
1421 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1423 SCM_I_BIG_MPZ (n2
));
1424 scm_remember_upto_here_2 (n1
, n2
);
1425 return scm_i_normbig (result_z
);
1428 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1431 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1436 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1437 (SCM x
, SCM y
, SCM rest
),
1438 "Return the bitwise XOR of the integer arguments. A bit is\n"
1439 "set in the result if it is set in an odd number of arguments.\n"
1441 "(logxor) @result{} 0\n"
1442 "(logxor 7) @result{} 7\n"
1443 "(logxor #b000 #b001 #b011) @result{} 2\n"
1444 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1446 #define FUNC_NAME s_scm_i_logxor
1448 while (!scm_is_null (rest
))
1449 { x
= scm_logxor (x
, y
);
1451 rest
= scm_cdr (rest
);
1453 return scm_logxor (x
, y
);
1457 #define s_scm_logxor s_scm_i_logxor
1459 SCM
scm_logxor (SCM n1
, SCM n2
)
1460 #define FUNC_NAME s_scm_logxor
1464 if (SCM_UNBNDP (n2
))
1466 if (SCM_UNBNDP (n1
))
1468 else if (SCM_NUMBERP (n1
))
1471 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1474 if (SCM_I_INUMP (n1
))
1476 nn1
= SCM_I_INUM (n1
);
1477 if (SCM_I_INUMP (n2
))
1479 long nn2
= SCM_I_INUM (n2
);
1480 return SCM_I_MAKINUM (nn1
^ nn2
);
1482 else if (SCM_BIGP (n2
))
1486 SCM result_z
= scm_i_mkbig ();
1488 mpz_init_set_si (nn1_z
, nn1
);
1489 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1490 scm_remember_upto_here_1 (n2
);
1492 return scm_i_normbig (result_z
);
1496 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1498 else if (SCM_BIGP (n1
))
1500 if (SCM_I_INUMP (n2
))
1503 nn1
= SCM_I_INUM (n1
);
1506 else if (SCM_BIGP (n2
))
1508 SCM result_z
= scm_i_mkbig ();
1509 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1511 SCM_I_BIG_MPZ (n2
));
1512 scm_remember_upto_here_2 (n1
, n2
);
1513 return scm_i_normbig (result_z
);
1516 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1519 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1524 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1526 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1527 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1528 "without actually calculating the @code{logand}, just testing\n"
1532 "(logtest #b0100 #b1011) @result{} #f\n"
1533 "(logtest #b0100 #b0111) @result{} #t\n"
1535 #define FUNC_NAME s_scm_logtest
1539 if (SCM_I_INUMP (j
))
1541 nj
= SCM_I_INUM (j
);
1542 if (SCM_I_INUMP (k
))
1544 long nk
= SCM_I_INUM (k
);
1545 return scm_from_bool (nj
& nk
);
1547 else if (SCM_BIGP (k
))
1555 mpz_init_set_si (nj_z
, nj
);
1556 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1557 scm_remember_upto_here_1 (k
);
1558 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1564 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1566 else if (SCM_BIGP (j
))
1568 if (SCM_I_INUMP (k
))
1571 nj
= SCM_I_INUM (j
);
1574 else if (SCM_BIGP (k
))
1578 mpz_init (result_z
);
1582 scm_remember_upto_here_2 (j
, k
);
1583 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1584 mpz_clear (result_z
);
1588 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1591 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1596 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1598 "Test whether bit number @var{index} in @var{j} is set.\n"
1599 "@var{index} starts from 0 for the least significant bit.\n"
1602 "(logbit? 0 #b1101) @result{} #t\n"
1603 "(logbit? 1 #b1101) @result{} #f\n"
1604 "(logbit? 2 #b1101) @result{} #t\n"
1605 "(logbit? 3 #b1101) @result{} #t\n"
1606 "(logbit? 4 #b1101) @result{} #f\n"
1608 #define FUNC_NAME s_scm_logbit_p
1610 unsigned long int iindex
;
1611 iindex
= scm_to_ulong (index
);
1613 if (SCM_I_INUMP (j
))
1615 /* bits above what's in an inum follow the sign bit */
1616 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1617 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1619 else if (SCM_BIGP (j
))
1621 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1622 scm_remember_upto_here_1 (j
);
1623 return scm_from_bool (val
);
1626 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1631 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1633 "Return the integer which is the ones-complement of the integer\n"
1637 "(number->string (lognot #b10000000) 2)\n"
1638 " @result{} \"-10000001\"\n"
1639 "(number->string (lognot #b0) 2)\n"
1640 " @result{} \"-1\"\n"
1642 #define FUNC_NAME s_scm_lognot
1644 if (SCM_I_INUMP (n
)) {
1645 /* No overflow here, just need to toggle all the bits making up the inum.
1646 Enhancement: No need to strip the tag and add it back, could just xor
1647 a block of 1 bits, if that worked with the various debug versions of
1649 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1651 } else if (SCM_BIGP (n
)) {
1652 SCM result
= scm_i_mkbig ();
1653 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1654 scm_remember_upto_here_1 (n
);
1658 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1663 /* returns 0 if IN is not an integer. OUT must already be
1666 coerce_to_big (SCM in
, mpz_t out
)
1669 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1670 else if (SCM_I_INUMP (in
))
1671 mpz_set_si (out
, SCM_I_INUM (in
));
1678 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1679 (SCM n
, SCM k
, SCM m
),
1680 "Return @var{n} raised to the integer exponent\n"
1681 "@var{k}, modulo @var{m}.\n"
1684 "(modulo-expt 2 3 5)\n"
1687 #define FUNC_NAME s_scm_modulo_expt
1693 /* There are two classes of error we might encounter --
1694 1) Math errors, which we'll report by calling scm_num_overflow,
1696 2) wrong-type errors, which of course we'll report by calling
1698 We don't report those errors immediately, however; instead we do
1699 some cleanup first. These variables tell us which error (if
1700 any) we should report after cleaning up.
1702 int report_overflow
= 0;
1704 int position_of_wrong_type
= 0;
1705 SCM value_of_wrong_type
= SCM_INUM0
;
1707 SCM result
= SCM_UNDEFINED
;
1713 if (scm_is_eq (m
, SCM_INUM0
))
1715 report_overflow
= 1;
1719 if (!coerce_to_big (n
, n_tmp
))
1721 value_of_wrong_type
= n
;
1722 position_of_wrong_type
= 1;
1726 if (!coerce_to_big (k
, k_tmp
))
1728 value_of_wrong_type
= k
;
1729 position_of_wrong_type
= 2;
1733 if (!coerce_to_big (m
, m_tmp
))
1735 value_of_wrong_type
= m
;
1736 position_of_wrong_type
= 3;
1740 /* if the exponent K is negative, and we simply call mpz_powm, we
1741 will get a divide-by-zero exception when an inverse 1/n mod m
1742 doesn't exist (or is not unique). Since exceptions are hard to
1743 handle, we'll attempt the inversion "by hand" -- that way, we get
1744 a simple failure code, which is easy to handle. */
1746 if (-1 == mpz_sgn (k_tmp
))
1748 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1750 report_overflow
= 1;
1753 mpz_neg (k_tmp
, k_tmp
);
1756 result
= scm_i_mkbig ();
1757 mpz_powm (SCM_I_BIG_MPZ (result
),
1762 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1763 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1770 if (report_overflow
)
1771 scm_num_overflow (FUNC_NAME
);
1773 if (position_of_wrong_type
)
1774 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1775 value_of_wrong_type
);
1777 return scm_i_normbig (result
);
1781 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1783 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1784 "exact integer, @var{n} can be any number.\n"
1786 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1787 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1788 "includes @math{0^0} is 1.\n"
1791 "(integer-expt 2 5) @result{} 32\n"
1792 "(integer-expt -3 3) @result{} -27\n"
1793 "(integer-expt 5 -3) @result{} 1/125\n"
1794 "(integer-expt 0 0) @result{} 1\n"
1796 #define FUNC_NAME s_scm_integer_expt
1799 SCM z_i2
= SCM_BOOL_F
;
1801 SCM acc
= SCM_I_MAKINUM (1L);
1803 /* 0^0 == 1 according to R5RS */
1804 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1805 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1806 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1807 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1809 if (SCM_I_INUMP (k
))
1810 i2
= SCM_I_INUM (k
);
1811 else if (SCM_BIGP (k
))
1813 z_i2
= scm_i_clonebig (k
, 1);
1814 scm_remember_upto_here_1 (k
);
1818 SCM_WRONG_TYPE_ARG (2, k
);
1822 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1824 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1825 n
= scm_divide (n
, SCM_UNDEFINED
);
1829 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1833 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1835 return scm_product (acc
, n
);
1837 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1838 acc
= scm_product (acc
, n
);
1839 n
= scm_product (n
, n
);
1840 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1848 n
= scm_divide (n
, SCM_UNDEFINED
);
1855 return scm_product (acc
, n
);
1857 acc
= scm_product (acc
, n
);
1858 n
= scm_product (n
, n
);
1865 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1867 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1868 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1870 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1871 "@var{cnt} is negative it's a division, rounded towards negative\n"
1872 "infinity. (Note that this is not the same rounding as\n"
1873 "@code{quotient} does.)\n"
1875 "With @var{n} viewed as an infinite precision twos complement,\n"
1876 "@code{ash} means a left shift introducing zero bits, or a right\n"
1877 "shift dropping bits.\n"
1880 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1881 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1883 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1884 "(ash -23 -2) @result{} -6\n"
1886 #define FUNC_NAME s_scm_ash
1889 bits_to_shift
= scm_to_long (cnt
);
1891 if (SCM_I_INUMP (n
))
1893 long nn
= SCM_I_INUM (n
);
1895 if (bits_to_shift
> 0)
1897 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1898 overflow a non-zero fixnum. For smaller shifts we check the
1899 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1900 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1901 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1907 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1909 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1912 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1916 SCM result
= scm_i_long2big (nn
);
1917 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1924 bits_to_shift
= -bits_to_shift
;
1925 if (bits_to_shift
>= SCM_LONG_BIT
)
1926 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1928 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1932 else if (SCM_BIGP (n
))
1936 if (bits_to_shift
== 0)
1939 result
= scm_i_mkbig ();
1940 if (bits_to_shift
>= 0)
1942 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1948 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1949 we have to allocate a bignum even if the result is going to be a
1951 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1953 return scm_i_normbig (result
);
1959 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1965 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1966 (SCM n
, SCM start
, SCM end
),
1967 "Return the integer composed of the @var{start} (inclusive)\n"
1968 "through @var{end} (exclusive) bits of @var{n}. The\n"
1969 "@var{start}th bit becomes the 0-th bit in the result.\n"
1972 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1973 " @result{} \"1010\"\n"
1974 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1975 " @result{} \"10110\"\n"
1977 #define FUNC_NAME s_scm_bit_extract
1979 unsigned long int istart
, iend
, bits
;
1980 istart
= scm_to_ulong (start
);
1981 iend
= scm_to_ulong (end
);
1982 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1984 /* how many bits to keep */
1985 bits
= iend
- istart
;
1987 if (SCM_I_INUMP (n
))
1989 long int in
= SCM_I_INUM (n
);
1991 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1992 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1993 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1995 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1997 /* Since we emulate two's complement encoded numbers, this
1998 * special case requires us to produce a result that has
1999 * more bits than can be stored in a fixnum.
2001 SCM result
= scm_i_long2big (in
);
2002 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2007 /* mask down to requisite bits */
2008 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2009 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2011 else if (SCM_BIGP (n
))
2016 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2020 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2021 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2022 such bits into a ulong. */
2023 result
= scm_i_mkbig ();
2024 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2025 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2026 result
= scm_i_normbig (result
);
2028 scm_remember_upto_here_1 (n
);
2032 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2037 static const char scm_logtab
[] = {
2038 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2041 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2043 "Return the number of bits in integer @var{n}. If integer is\n"
2044 "positive, the 1-bits in its binary representation are counted.\n"
2045 "If negative, the 0-bits in its two's-complement binary\n"
2046 "representation are counted. If 0, 0 is returned.\n"
2049 "(logcount #b10101010)\n"
2056 #define FUNC_NAME s_scm_logcount
2058 if (SCM_I_INUMP (n
))
2060 unsigned long int c
= 0;
2061 long int nn
= SCM_I_INUM (n
);
2066 c
+= scm_logtab
[15 & nn
];
2069 return SCM_I_MAKINUM (c
);
2071 else if (SCM_BIGP (n
))
2073 unsigned long count
;
2074 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2075 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2077 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2078 scm_remember_upto_here_1 (n
);
2079 return SCM_I_MAKINUM (count
);
2082 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2087 static const char scm_ilentab
[] = {
2088 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2092 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2094 "Return the number of bits necessary to represent @var{n}.\n"
2097 "(integer-length #b10101010)\n"
2099 "(integer-length 0)\n"
2101 "(integer-length #b1111)\n"
2104 #define FUNC_NAME s_scm_integer_length
2106 if (SCM_I_INUMP (n
))
2108 unsigned long int c
= 0;
2110 long int nn
= SCM_I_INUM (n
);
2116 l
= scm_ilentab
[15 & nn
];
2119 return SCM_I_MAKINUM (c
- 4 + l
);
2121 else if (SCM_BIGP (n
))
2123 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2124 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2125 1 too big, so check for that and adjust. */
2126 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2127 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2128 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2129 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2131 scm_remember_upto_here_1 (n
);
2132 return SCM_I_MAKINUM (size
);
2135 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2139 /*** NUMBERS -> STRINGS ***/
2140 #define SCM_MAX_DBL_PREC 60
2141 #define SCM_MAX_DBL_RADIX 36
2143 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2144 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2145 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2148 void init_dblprec(int *prec
, int radix
) {
2149 /* determine floating point precision by adding successively
2150 smaller increments to 1.0 until it is considered == 1.0 */
2151 double f
= ((double)1.0)/radix
;
2152 double fsum
= 1.0 + f
;
2157 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2169 void init_fx_radix(double *fx_list
, int radix
)
2171 /* initialize a per-radix list of tolerances. When added
2172 to a number < 1.0, we can determine if we should raund
2173 up and quit converting a number to a string. */
2177 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2178 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2181 /* use this array as a way to generate a single digit */
2182 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2185 idbl2str (double f
, char *a
, int radix
)
2187 int efmt
, dpt
, d
, i
, wp
;
2189 #ifdef DBL_MIN_10_EXP
2192 #endif /* DBL_MIN_10_EXP */
2197 radix
> SCM_MAX_DBL_RADIX
)
2199 /* revert to existing behavior */
2203 wp
= scm_dblprec
[radix
-2];
2204 fx
= fx_per_radix
[radix
-2];
2208 #ifdef HAVE_COPYSIGN
2209 double sgn
= copysign (1.0, f
);
2214 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2220 strcpy (a
, "-inf.0");
2222 strcpy (a
, "+inf.0");
2225 else if (xisnan (f
))
2227 strcpy (a
, "+nan.0");
2237 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2238 make-uniform-vector, from causing infinite loops. */
2239 /* just do the checking...if it passes, we do the conversion for our
2240 radix again below */
2247 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2255 while (f_cpy
> 10.0)
2258 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2279 if (f
+ fx
[wp
] >= radix
)
2286 /* adding 9999 makes this equivalent to abs(x) % 3 */
2287 dpt
= (exp
+ 9999) % 3;
2291 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2313 a
[ch
++] = number_chars
[d
];
2316 if (f
+ fx
[wp
] >= 1.0)
2318 a
[ch
- 1] = number_chars
[d
+1];
2330 if ((dpt
> 4) && (exp
> 6))
2332 d
= (a
[0] == '-' ? 2 : 1);
2333 for (i
= ch
++; i
> d
; i
--)
2346 if (a
[ch
- 1] == '.')
2347 a
[ch
++] = '0'; /* trailing zero */
2356 for (i
= radix
; i
<= exp
; i
*= radix
);
2357 for (i
/= radix
; i
; i
/= radix
)
2359 a
[ch
++] = number_chars
[exp
/ i
];
2368 icmplx2str (double real
, double imag
, char *str
, int radix
)
2372 i
= idbl2str (real
, str
, radix
);
2375 /* Don't output a '+' for negative numbers or for Inf and
2376 NaN. They will provide their own sign. */
2377 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2379 i
+= idbl2str (imag
, &str
[i
], radix
);
2386 iflo2str (SCM flt
, char *str
, int radix
)
2389 if (SCM_REALP (flt
))
2390 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2392 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2397 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2398 characters in the result.
2400 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2402 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2407 return scm_iuint2str (-num
, rad
, p
) + 1;
2410 return scm_iuint2str (num
, rad
, p
);
2413 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2414 characters in the result.
2416 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2418 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2422 scm_t_uintmax n
= num
;
2424 for (n
/= rad
; n
> 0; n
/= rad
)
2434 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2439 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2441 "Return a string holding the external representation of the\n"
2442 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2443 "inexact, a radix of 10 will be used.")
2444 #define FUNC_NAME s_scm_number_to_string
2448 if (SCM_UNBNDP (radix
))
2451 base
= scm_to_signed_integer (radix
, 2, 36);
2453 if (SCM_I_INUMP (n
))
2455 char num_buf
[SCM_INTBUFLEN
];
2456 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2457 return scm_from_locale_stringn (num_buf
, length
);
2459 else if (SCM_BIGP (n
))
2461 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2462 scm_remember_upto_here_1 (n
);
2463 return scm_take_locale_string (str
);
2465 else if (SCM_FRACTIONP (n
))
2467 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2468 scm_from_locale_string ("/"),
2469 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2471 else if (SCM_INEXACTP (n
))
2473 char num_buf
[FLOBUFLEN
];
2474 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2477 SCM_WRONG_TYPE_ARG (1, n
);
2482 /* These print routines used to be stubbed here so that scm_repl.c
2483 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2486 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2488 char num_buf
[FLOBUFLEN
];
2489 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2494 scm_i_print_double (double val
, SCM port
)
2496 char num_buf
[FLOBUFLEN
];
2497 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2501 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2504 char num_buf
[FLOBUFLEN
];
2505 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2510 scm_i_print_complex (double real
, double imag
, SCM port
)
2512 char num_buf
[FLOBUFLEN
];
2513 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2517 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2520 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2521 scm_lfwrite_str (str
, port
);
2522 scm_remember_upto_here_1 (str
);
2527 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2529 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2530 scm_remember_upto_here_1 (exp
);
2531 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2535 /*** END nums->strs ***/
2538 /*** STRINGS -> NUMBERS ***/
2540 /* The following functions implement the conversion from strings to numbers.
2541 * The implementation somehow follows the grammar for numbers as it is given
2542 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2543 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2544 * points should be noted about the implementation:
2545 * * Each function keeps a local index variable 'idx' that points at the
2546 * current position within the parsed string. The global index is only
2547 * updated if the function could parse the corresponding syntactic unit
2549 * * Similarly, the functions keep track of indicators of inexactness ('#',
2550 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2551 * global exactness information is only updated after each part has been
2552 * successfully parsed.
2553 * * Sequences of digits are parsed into temporary variables holding fixnums.
2554 * Only if these fixnums would overflow, the result variables are updated
2555 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2556 * the temporary variables holding the fixnums are cleared, and the process
2557 * starts over again. If for example fixnums were able to store five decimal
2558 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2559 * and the result was computed as 12345 * 100000 + 67890. In other words,
2560 * only every five digits two bignum operations were performed.
2563 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2565 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2567 /* In non ASCII-style encodings the following macro might not work. */
2568 #define XDIGIT2UINT(d) \
2569 (uc_is_property_decimal_digit ((int) (unsigned char) d) \
2571 : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
2574 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2575 unsigned int radix
, enum t_exactness
*p_exactness
)
2577 unsigned int idx
= *p_idx
;
2578 unsigned int hash_seen
= 0;
2579 scm_t_bits shift
= 1;
2581 unsigned int digit_value
;
2584 size_t len
= scm_i_string_length (mem
);
2589 c
= scm_i_string_ref (mem
, idx
);
2590 if (!uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2592 digit_value
= XDIGIT2UINT (c
);
2593 if (digit_value
>= radix
)
2597 result
= SCM_I_MAKINUM (digit_value
);
2600 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2601 if (uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2605 digit_value
= XDIGIT2UINT (c
);
2606 if (digit_value
>= radix
)
2618 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2620 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2622 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2629 shift
= shift
* radix
;
2630 add
= add
* radix
+ digit_value
;
2635 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2637 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2641 *p_exactness
= INEXACT
;
2647 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2648 * covers the parts of the rules that start at a potential point. The value
2649 * of the digits up to the point have been parsed by the caller and are given
2650 * in variable result. The content of *p_exactness indicates, whether a hash
2651 * has already been seen in the digits before the point.
2654 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2657 mem2decimal_from_point (SCM result
, SCM mem
,
2658 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2660 unsigned int idx
= *p_idx
;
2661 enum t_exactness x
= *p_exactness
;
2662 size_t len
= scm_i_string_length (mem
);
2667 if (scm_i_string_ref (mem
, idx
) == '.')
2669 scm_t_bits shift
= 1;
2671 unsigned int digit_value
;
2672 SCM big_shift
= SCM_I_MAKINUM (1);
2677 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2678 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2683 digit_value
= DIGIT2UINT (c
);
2694 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2696 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2697 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2699 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2707 add
= add
* 10 + digit_value
;
2713 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2714 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2715 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2718 result
= scm_divide (result
, big_shift
);
2720 /* We've seen a decimal point, thus the value is implicitly inexact. */
2732 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2734 switch (scm_i_string_ref (mem
, idx
))
2746 c
= scm_i_string_ref (mem
, idx
);
2754 c
= scm_i_string_ref (mem
, idx
);
2763 c
= scm_i_string_ref (mem
, idx
);
2768 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2772 exponent
= DIGIT2UINT (c
);
2775 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2776 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2779 if (exponent
<= SCM_MAXEXP
)
2780 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2786 if (exponent
> SCM_MAXEXP
)
2788 size_t exp_len
= idx
- start
;
2789 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2790 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2791 scm_out_of_range ("string->number", exp_num
);
2794 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2796 result
= scm_product (result
, e
);
2798 result
= scm_divide2real (result
, e
);
2800 /* We've seen an exponent, thus the value is implicitly inexact. */
2818 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2821 mem2ureal (SCM mem
, unsigned int *p_idx
,
2822 unsigned int radix
, enum t_exactness
*p_exactness
)
2824 unsigned int idx
= *p_idx
;
2826 size_t len
= scm_i_string_length (mem
);
2828 /* Start off believing that the number will be exact. This changes
2829 to INEXACT if we see a decimal point or a hash. */
2830 enum t_exactness x
= EXACT
;
2835 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2841 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2843 /* Cobble up the fractional part. We might want to set the
2844 NaN's mantissa from it. */
2846 mem2uinteger (mem
, &idx
, 10, &x
);
2851 if (scm_i_string_ref (mem
, idx
) == '.')
2855 else if (idx
+ 1 == len
)
2857 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2860 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2867 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2868 if (scm_is_false (uinteger
))
2873 else if (scm_i_string_ref (mem
, idx
) == '/')
2881 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2882 if (scm_is_false (divisor
))
2885 /* both are int/big here, I assume */
2886 result
= scm_i_make_ratio (uinteger
, divisor
);
2888 else if (radix
== 10)
2890 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2891 if (scm_is_false (result
))
2900 /* Update *p_exactness if the number just read was inexact. This is
2901 important for complex numbers, so that a complex number is
2902 treated as inexact overall if either its real or imaginary part
2908 /* When returning an inexact zero, make sure it is represented as a
2909 floating point value so that we can change its sign.
2911 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2912 result
= scm_from_double (0.0);
2918 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2921 mem2complex (SCM mem
, unsigned int idx
,
2922 unsigned int radix
, enum t_exactness
*p_exactness
)
2927 size_t len
= scm_i_string_length (mem
);
2932 c
= scm_i_string_ref (mem
, idx
);
2947 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2948 if (scm_is_false (ureal
))
2950 /* input must be either +i or -i */
2955 if (scm_i_string_ref (mem
, idx
) == 'i'
2956 || scm_i_string_ref (mem
, idx
) == 'I')
2962 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2969 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2970 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2975 c
= scm_i_string_ref (mem
, idx
);
2979 /* either +<ureal>i or -<ureal>i */
2986 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2989 /* polar input: <real>@<real>. */
3000 c
= scm_i_string_ref (mem
, idx
);
3018 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3019 if (scm_is_false (angle
))
3024 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3025 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3027 result
= scm_make_polar (ureal
, angle
);
3032 /* expecting input matching <real>[+-]<ureal>?i */
3039 int sign
= (c
== '+') ? 1 : -1;
3040 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3042 if (scm_is_false (imag
))
3043 imag
= SCM_I_MAKINUM (sign
);
3044 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3045 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3049 if (scm_i_string_ref (mem
, idx
) != 'i'
3050 && scm_i_string_ref (mem
, idx
) != 'I')
3057 return scm_make_rectangular (ureal
, imag
);
3066 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3068 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3071 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3073 unsigned int idx
= 0;
3074 unsigned int radix
= NO_RADIX
;
3075 enum t_exactness forced_x
= NO_EXACTNESS
;
3076 enum t_exactness implicit_x
= EXACT
;
3078 size_t len
= scm_i_string_length (mem
);
3080 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3081 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3083 switch (scm_i_string_ref (mem
, idx
+ 1))
3086 if (radix
!= NO_RADIX
)
3091 if (radix
!= NO_RADIX
)
3096 if (forced_x
!= NO_EXACTNESS
)
3101 if (forced_x
!= NO_EXACTNESS
)
3106 if (radix
!= NO_RADIX
)
3111 if (radix
!= NO_RADIX
)
3121 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3122 if (radix
== NO_RADIX
)
3123 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3125 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3127 if (scm_is_false (result
))
3133 if (SCM_INEXACTP (result
))
3134 return scm_inexact_to_exact (result
);
3138 if (SCM_INEXACTP (result
))
3141 return scm_exact_to_inexact (result
);
3144 if (implicit_x
== INEXACT
)
3146 if (SCM_INEXACTP (result
))
3149 return scm_exact_to_inexact (result
);
3157 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3158 unsigned int default_radix
)
3160 SCM str
= scm_from_locale_stringn (mem
, len
);
3162 return scm_i_string_to_number (str
, default_radix
);
3166 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3167 (SCM string
, SCM radix
),
3168 "Return a number of the maximally precise representation\n"
3169 "expressed by the given @var{string}. @var{radix} must be an\n"
3170 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3171 "is a default radix that may be overridden by an explicit radix\n"
3172 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3173 "supplied, then the default radix is 10. If string is not a\n"
3174 "syntactically valid notation for a number, then\n"
3175 "@code{string->number} returns @code{#f}.")
3176 #define FUNC_NAME s_scm_string_to_number
3180 SCM_VALIDATE_STRING (1, string
);
3182 if (SCM_UNBNDP (radix
))
3185 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3187 answer
= scm_i_string_to_number (string
, base
);
3188 scm_remember_upto_here_1 (string
);
3194 /*** END strs->nums ***/
3198 scm_bigequal (SCM x
, SCM y
)
3200 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3201 scm_remember_upto_here_2 (x
, y
);
3202 return scm_from_bool (0 == result
);
3206 scm_real_equalp (SCM x
, SCM y
)
3208 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3212 scm_complex_equalp (SCM x
, SCM y
)
3214 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3215 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3219 scm_i_fraction_equalp (SCM x
, SCM y
)
3221 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3222 SCM_FRACTION_NUMERATOR (y
)))
3223 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3224 SCM_FRACTION_DENOMINATOR (y
))))
3231 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3233 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3235 #define FUNC_NAME s_scm_number_p
3237 return scm_from_bool (SCM_NUMBERP (x
));
3241 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3243 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3244 "otherwise. Note that the sets of real, rational and integer\n"
3245 "values form subsets of the set of complex numbers, i. e. the\n"
3246 "predicate will also be fulfilled if @var{x} is a real,\n"
3247 "rational or integer number.")
3248 #define FUNC_NAME s_scm_complex_p
3250 /* all numbers are complex. */
3251 return scm_number_p (x
);
3255 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3257 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3258 "otherwise. Note that the set of integer values forms a subset of\n"
3259 "the set of real numbers, i. e. the predicate will also be\n"
3260 "fulfilled if @var{x} is an integer number.")
3261 #define FUNC_NAME s_scm_real_p
3263 /* we can't represent irrational numbers. */
3264 return scm_rational_p (x
);
3268 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3270 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3271 "otherwise. Note that the set of integer values forms a subset of\n"
3272 "the set of rational numbers, i. e. the predicate will also be\n"
3273 "fulfilled if @var{x} is an integer number.")
3274 #define FUNC_NAME s_scm_rational_p
3276 if (SCM_I_INUMP (x
))
3278 else if (SCM_IMP (x
))
3280 else if (SCM_BIGP (x
))
3282 else if (SCM_FRACTIONP (x
))
3284 else if (SCM_REALP (x
))
3285 /* due to their limited precision, all floating point numbers are
3286 rational as well. */
3293 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3295 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3297 #define FUNC_NAME s_scm_integer_p
3300 if (SCM_I_INUMP (x
))
3306 if (!SCM_INEXACTP (x
))
3308 if (SCM_COMPLEXP (x
))
3310 r
= SCM_REAL_VALUE (x
);
3311 /* +/-inf passes r==floor(r), making those #t */
3319 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3321 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3323 #define FUNC_NAME s_scm_inexact_p
3325 if (SCM_INEXACTP (x
))
3327 if (SCM_NUMBERP (x
))
3329 SCM_WRONG_TYPE_ARG (1, x
);
3334 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3335 /* "Return @code{#t} if all parameters are numerically equal." */
3337 scm_num_eq_p (SCM x
, SCM y
)
3340 if (SCM_I_INUMP (x
))
3342 long xx
= SCM_I_INUM (x
);
3343 if (SCM_I_INUMP (y
))
3345 long yy
= SCM_I_INUM (y
);
3346 return scm_from_bool (xx
== yy
);
3348 else if (SCM_BIGP (y
))
3350 else if (SCM_REALP (y
))
3352 /* On a 32-bit system an inum fits a double, we can cast the inum
3353 to a double and compare.
3355 But on a 64-bit system an inum is bigger than a double and
3356 casting it to a double (call that dxx) will round. dxx is at
3357 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3358 an integer and fits a long. So we cast yy to a long and
3359 compare with plain xx.
3361 An alternative (for any size system actually) would be to check
3362 yy is an integer (with floor) and is in range of an inum
3363 (compare against appropriate powers of 2) then test
3364 xx==(long)yy. It's just a matter of which casts/comparisons
3365 might be fastest or easiest for the cpu. */
3367 double yy
= SCM_REAL_VALUE (y
);
3368 return scm_from_bool ((double) xx
== yy
3369 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3370 || xx
== (long) yy
));
3372 else if (SCM_COMPLEXP (y
))
3373 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3374 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3375 else if (SCM_FRACTIONP (y
))
3378 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3380 else if (SCM_BIGP (x
))
3382 if (SCM_I_INUMP (y
))
3384 else if (SCM_BIGP (y
))
3386 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3387 scm_remember_upto_here_2 (x
, y
);
3388 return scm_from_bool (0 == cmp
);
3390 else if (SCM_REALP (y
))
3393 if (xisnan (SCM_REAL_VALUE (y
)))
3395 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3396 scm_remember_upto_here_1 (x
);
3397 return scm_from_bool (0 == cmp
);
3399 else if (SCM_COMPLEXP (y
))
3402 if (0.0 != SCM_COMPLEX_IMAG (y
))
3404 if (xisnan (SCM_COMPLEX_REAL (y
)))
3406 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3407 scm_remember_upto_here_1 (x
);
3408 return scm_from_bool (0 == cmp
);
3410 else if (SCM_FRACTIONP (y
))
3413 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3415 else if (SCM_REALP (x
))
3417 double xx
= SCM_REAL_VALUE (x
);
3418 if (SCM_I_INUMP (y
))
3420 /* see comments with inum/real above */
3421 long yy
= SCM_I_INUM (y
);
3422 return scm_from_bool (xx
== (double) yy
3423 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3424 || (long) xx
== yy
));
3426 else if (SCM_BIGP (y
))
3429 if (xisnan (SCM_REAL_VALUE (x
)))
3431 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3432 scm_remember_upto_here_1 (y
);
3433 return scm_from_bool (0 == cmp
);
3435 else if (SCM_REALP (y
))
3436 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3437 else if (SCM_COMPLEXP (y
))
3438 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3439 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3440 else if (SCM_FRACTIONP (y
))
3442 double xx
= SCM_REAL_VALUE (x
);
3446 return scm_from_bool (xx
< 0.0);
3447 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3451 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3453 else if (SCM_COMPLEXP (x
))
3455 if (SCM_I_INUMP (y
))
3456 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3457 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3458 else if (SCM_BIGP (y
))
3461 if (0.0 != SCM_COMPLEX_IMAG (x
))
3463 if (xisnan (SCM_COMPLEX_REAL (x
)))
3465 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3466 scm_remember_upto_here_1 (y
);
3467 return scm_from_bool (0 == cmp
);
3469 else if (SCM_REALP (y
))
3470 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3471 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3472 else if (SCM_COMPLEXP (y
))
3473 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3474 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3475 else if (SCM_FRACTIONP (y
))
3478 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3480 xx
= SCM_COMPLEX_REAL (x
);
3484 return scm_from_bool (xx
< 0.0);
3485 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3489 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3491 else if (SCM_FRACTIONP (x
))
3493 if (SCM_I_INUMP (y
))
3495 else if (SCM_BIGP (y
))
3497 else if (SCM_REALP (y
))
3499 double yy
= SCM_REAL_VALUE (y
);
3503 return scm_from_bool (0.0 < yy
);
3504 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3507 else if (SCM_COMPLEXP (y
))
3510 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3512 yy
= SCM_COMPLEX_REAL (y
);
3516 return scm_from_bool (0.0 < yy
);
3517 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3520 else if (SCM_FRACTIONP (y
))
3521 return scm_i_fraction_equalp (x
, y
);
3523 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3526 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3530 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3531 done are good for inums, but for bignums an answer can almost always be
3532 had by just examining a few high bits of the operands, as done by GMP in
3533 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3534 of the float exponent to take into account. */
3536 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3537 /* "Return @code{#t} if the list of parameters is monotonically\n"
3541 scm_less_p (SCM x
, SCM y
)
3544 if (SCM_I_INUMP (x
))
3546 long xx
= SCM_I_INUM (x
);
3547 if (SCM_I_INUMP (y
))
3549 long yy
= SCM_I_INUM (y
);
3550 return scm_from_bool (xx
< yy
);
3552 else if (SCM_BIGP (y
))
3554 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3555 scm_remember_upto_here_1 (y
);
3556 return scm_from_bool (sgn
> 0);
3558 else if (SCM_REALP (y
))
3559 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3560 else if (SCM_FRACTIONP (y
))
3562 /* "x < a/b" becomes "x*b < a" */
3564 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3565 y
= SCM_FRACTION_NUMERATOR (y
);
3569 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3571 else if (SCM_BIGP (x
))
3573 if (SCM_I_INUMP (y
))
3575 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3576 scm_remember_upto_here_1 (x
);
3577 return scm_from_bool (sgn
< 0);
3579 else if (SCM_BIGP (y
))
3581 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3582 scm_remember_upto_here_2 (x
, y
);
3583 return scm_from_bool (cmp
< 0);
3585 else if (SCM_REALP (y
))
3588 if (xisnan (SCM_REAL_VALUE (y
)))
3590 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3591 scm_remember_upto_here_1 (x
);
3592 return scm_from_bool (cmp
< 0);
3594 else if (SCM_FRACTIONP (y
))
3597 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3599 else if (SCM_REALP (x
))
3601 if (SCM_I_INUMP (y
))
3602 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3603 else if (SCM_BIGP (y
))
3606 if (xisnan (SCM_REAL_VALUE (x
)))
3608 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3609 scm_remember_upto_here_1 (y
);
3610 return scm_from_bool (cmp
> 0);
3612 else if (SCM_REALP (y
))
3613 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3614 else if (SCM_FRACTIONP (y
))
3616 double xx
= SCM_REAL_VALUE (x
);
3620 return scm_from_bool (xx
< 0.0);
3621 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3625 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3627 else if (SCM_FRACTIONP (x
))
3629 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3631 /* "a/b < y" becomes "a < y*b" */
3632 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3633 x
= SCM_FRACTION_NUMERATOR (x
);
3636 else if (SCM_REALP (y
))
3638 double yy
= SCM_REAL_VALUE (y
);
3642 return scm_from_bool (0.0 < yy
);
3643 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3646 else if (SCM_FRACTIONP (y
))
3648 /* "a/b < c/d" becomes "a*d < c*b" */
3649 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3650 SCM_FRACTION_DENOMINATOR (y
));
3651 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3652 SCM_FRACTION_DENOMINATOR (x
));
3658 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3661 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3665 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3666 /* "Return @code{#t} if the list of parameters is monotonically\n"
3669 #define FUNC_NAME s_scm_gr_p
3671 scm_gr_p (SCM x
, SCM y
)
3673 if (!SCM_NUMBERP (x
))
3674 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3675 else if (!SCM_NUMBERP (y
))
3676 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3678 return scm_less_p (y
, x
);
3683 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3684 /* "Return @code{#t} if the list of parameters is monotonically\n"
3687 #define FUNC_NAME s_scm_leq_p
3689 scm_leq_p (SCM x
, SCM y
)
3691 if (!SCM_NUMBERP (x
))
3692 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3693 else if (!SCM_NUMBERP (y
))
3694 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3695 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3698 return scm_not (scm_less_p (y
, x
));
3703 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3704 /* "Return @code{#t} if the list of parameters is monotonically\n"
3707 #define FUNC_NAME s_scm_geq_p
3709 scm_geq_p (SCM x
, SCM y
)
3711 if (!SCM_NUMBERP (x
))
3712 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3713 else if (!SCM_NUMBERP (y
))
3714 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3715 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3718 return scm_not (scm_less_p (x
, y
));
3723 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3724 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3730 if (SCM_I_INUMP (z
))
3731 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3732 else if (SCM_BIGP (z
))
3734 else if (SCM_REALP (z
))
3735 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3736 else if (SCM_COMPLEXP (z
))
3737 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3738 && SCM_COMPLEX_IMAG (z
) == 0.0);
3739 else if (SCM_FRACTIONP (z
))
3742 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3746 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3747 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3751 scm_positive_p (SCM x
)
3753 if (SCM_I_INUMP (x
))
3754 return scm_from_bool (SCM_I_INUM (x
) > 0);
3755 else if (SCM_BIGP (x
))
3757 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3758 scm_remember_upto_here_1 (x
);
3759 return scm_from_bool (sgn
> 0);
3761 else if (SCM_REALP (x
))
3762 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3763 else if (SCM_FRACTIONP (x
))
3764 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3766 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3770 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3771 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3775 scm_negative_p (SCM x
)
3777 if (SCM_I_INUMP (x
))
3778 return scm_from_bool (SCM_I_INUM (x
) < 0);
3779 else if (SCM_BIGP (x
))
3781 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3782 scm_remember_upto_here_1 (x
);
3783 return scm_from_bool (sgn
< 0);
3785 else if (SCM_REALP (x
))
3786 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3787 else if (SCM_FRACTIONP (x
))
3788 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3790 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3794 /* scm_min and scm_max return an inexact when either argument is inexact, as
3795 required by r5rs. On that basis, for exact/inexact combinations the
3796 exact is converted to inexact to compare and possibly return. This is
3797 unlike scm_less_p above which takes some trouble to preserve all bits in
3798 its test, such trouble is not required for min and max. */
3800 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3801 (SCM x
, SCM y
, SCM rest
),
3802 "Return the maximum of all parameter values.")
3803 #define FUNC_NAME s_scm_i_max
3805 while (!scm_is_null (rest
))
3806 { x
= scm_max (x
, y
);
3808 rest
= scm_cdr (rest
);
3810 return scm_max (x
, y
);
3814 #define s_max s_scm_i_max
3815 #define g_max g_scm_i_max
3818 scm_max (SCM x
, SCM y
)
3823 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3824 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3827 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3830 if (SCM_I_INUMP (x
))
3832 long xx
= SCM_I_INUM (x
);
3833 if (SCM_I_INUMP (y
))
3835 long yy
= SCM_I_INUM (y
);
3836 return (xx
< yy
) ? y
: x
;
3838 else if (SCM_BIGP (y
))
3840 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3841 scm_remember_upto_here_1 (y
);
3842 return (sgn
< 0) ? x
: y
;
3844 else if (SCM_REALP (y
))
3847 /* if y==NaN then ">" is false and we return NaN */
3848 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3850 else if (SCM_FRACTIONP (y
))
3853 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3856 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3858 else if (SCM_BIGP (x
))
3860 if (SCM_I_INUMP (y
))
3862 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3863 scm_remember_upto_here_1 (x
);
3864 return (sgn
< 0) ? y
: x
;
3866 else if (SCM_BIGP (y
))
3868 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3869 scm_remember_upto_here_2 (x
, y
);
3870 return (cmp
> 0) ? x
: y
;
3872 else if (SCM_REALP (y
))
3874 /* if y==NaN then xx>yy is false, so we return the NaN y */
3877 xx
= scm_i_big2dbl (x
);
3878 yy
= SCM_REAL_VALUE (y
);
3879 return (xx
> yy
? scm_from_double (xx
) : y
);
3881 else if (SCM_FRACTIONP (y
))
3886 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3888 else if (SCM_REALP (x
))
3890 if (SCM_I_INUMP (y
))
3892 double z
= SCM_I_INUM (y
);
3893 /* if x==NaN then "<" is false and we return NaN */
3894 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3896 else if (SCM_BIGP (y
))
3901 else if (SCM_REALP (y
))
3903 /* if x==NaN then our explicit check means we return NaN
3904 if y==NaN then ">" is false and we return NaN
3905 calling isnan is unavoidable, since it's the only way to know
3906 which of x or y causes any compares to be false */
3907 double xx
= SCM_REAL_VALUE (x
);
3908 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3910 else if (SCM_FRACTIONP (y
))
3912 double yy
= scm_i_fraction2double (y
);
3913 double xx
= SCM_REAL_VALUE (x
);
3914 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3917 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3919 else if (SCM_FRACTIONP (x
))
3921 if (SCM_I_INUMP (y
))
3925 else if (SCM_BIGP (y
))
3929 else if (SCM_REALP (y
))
3931 double xx
= scm_i_fraction2double (x
);
3932 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3934 else if (SCM_FRACTIONP (y
))
3939 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3942 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3946 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
3947 (SCM x
, SCM y
, SCM rest
),
3948 "Return the minimum of all parameter values.")
3949 #define FUNC_NAME s_scm_i_min
3951 while (!scm_is_null (rest
))
3952 { x
= scm_min (x
, y
);
3954 rest
= scm_cdr (rest
);
3956 return scm_min (x
, y
);
3960 #define s_min s_scm_i_min
3961 #define g_min g_scm_i_min
3964 scm_min (SCM x
, SCM y
)
3969 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3970 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3973 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3976 if (SCM_I_INUMP (x
))
3978 long xx
= SCM_I_INUM (x
);
3979 if (SCM_I_INUMP (y
))
3981 long yy
= SCM_I_INUM (y
);
3982 return (xx
< yy
) ? x
: y
;
3984 else if (SCM_BIGP (y
))
3986 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3987 scm_remember_upto_here_1 (y
);
3988 return (sgn
< 0) ? y
: x
;
3990 else if (SCM_REALP (y
))
3993 /* if y==NaN then "<" is false and we return NaN */
3994 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3996 else if (SCM_FRACTIONP (y
))
3999 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4002 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4004 else if (SCM_BIGP (x
))
4006 if (SCM_I_INUMP (y
))
4008 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4009 scm_remember_upto_here_1 (x
);
4010 return (sgn
< 0) ? x
: y
;
4012 else if (SCM_BIGP (y
))
4014 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4015 scm_remember_upto_here_2 (x
, y
);
4016 return (cmp
> 0) ? y
: x
;
4018 else if (SCM_REALP (y
))
4020 /* if y==NaN then xx<yy is false, so we return the NaN y */
4023 xx
= scm_i_big2dbl (x
);
4024 yy
= SCM_REAL_VALUE (y
);
4025 return (xx
< yy
? scm_from_double (xx
) : y
);
4027 else if (SCM_FRACTIONP (y
))
4032 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4034 else if (SCM_REALP (x
))
4036 if (SCM_I_INUMP (y
))
4038 double z
= SCM_I_INUM (y
);
4039 /* if x==NaN then "<" is false and we return NaN */
4040 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4042 else if (SCM_BIGP (y
))
4047 else if (SCM_REALP (y
))
4049 /* if x==NaN then our explicit check means we return NaN
4050 if y==NaN then "<" is false and we return NaN
4051 calling isnan is unavoidable, since it's the only way to know
4052 which of x or y causes any compares to be false */
4053 double xx
= SCM_REAL_VALUE (x
);
4054 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4056 else if (SCM_FRACTIONP (y
))
4058 double yy
= scm_i_fraction2double (y
);
4059 double xx
= SCM_REAL_VALUE (x
);
4060 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4063 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4065 else if (SCM_FRACTIONP (x
))
4067 if (SCM_I_INUMP (y
))
4071 else if (SCM_BIGP (y
))
4075 else if (SCM_REALP (y
))
4077 double xx
= scm_i_fraction2double (x
);
4078 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4080 else if (SCM_FRACTIONP (y
))
4085 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4088 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4092 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4093 (SCM x
, SCM y
, SCM rest
),
4094 "Return the sum of all parameter values. Return 0 if called without\n"
4096 #define FUNC_NAME s_scm_i_sum
4098 while (!scm_is_null (rest
))
4099 { x
= scm_sum (x
, y
);
4101 rest
= scm_cdr (rest
);
4103 return scm_sum (x
, y
);
4107 #define s_sum s_scm_i_sum
4108 #define g_sum g_scm_i_sum
4111 scm_sum (SCM x
, SCM y
)
4113 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4115 if (SCM_NUMBERP (x
)) return x
;
4116 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4117 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4120 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4122 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4124 long xx
= SCM_I_INUM (x
);
4125 long yy
= SCM_I_INUM (y
);
4126 long int z
= xx
+ yy
;
4127 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4129 else if (SCM_BIGP (y
))
4134 else if (SCM_REALP (y
))
4136 long int xx
= SCM_I_INUM (x
);
4137 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4139 else if (SCM_COMPLEXP (y
))
4141 long int xx
= SCM_I_INUM (x
);
4142 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4143 SCM_COMPLEX_IMAG (y
));
4145 else if (SCM_FRACTIONP (y
))
4146 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4147 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4148 SCM_FRACTION_DENOMINATOR (y
));
4150 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4151 } else if (SCM_BIGP (x
))
4153 if (SCM_I_INUMP (y
))
4158 inum
= SCM_I_INUM (y
);
4161 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4164 SCM result
= scm_i_mkbig ();
4165 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4166 scm_remember_upto_here_1 (x
);
4167 /* we know the result will have to be a bignum */
4170 return scm_i_normbig (result
);
4174 SCM result
= scm_i_mkbig ();
4175 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4176 scm_remember_upto_here_1 (x
);
4177 /* we know the result will have to be a bignum */
4180 return scm_i_normbig (result
);
4183 else if (SCM_BIGP (y
))
4185 SCM result
= scm_i_mkbig ();
4186 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4187 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4188 mpz_add (SCM_I_BIG_MPZ (result
),
4191 scm_remember_upto_here_2 (x
, y
);
4192 /* we know the result will have to be a bignum */
4195 return scm_i_normbig (result
);
4197 else if (SCM_REALP (y
))
4199 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4200 scm_remember_upto_here_1 (x
);
4201 return scm_from_double (result
);
4203 else if (SCM_COMPLEXP (y
))
4205 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4206 + SCM_COMPLEX_REAL (y
));
4207 scm_remember_upto_here_1 (x
);
4208 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4210 else if (SCM_FRACTIONP (y
))
4211 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4212 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4213 SCM_FRACTION_DENOMINATOR (y
));
4215 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4217 else if (SCM_REALP (x
))
4219 if (SCM_I_INUMP (y
))
4220 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4221 else if (SCM_BIGP (y
))
4223 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4224 scm_remember_upto_here_1 (y
);
4225 return scm_from_double (result
);
4227 else if (SCM_REALP (y
))
4228 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4229 else if (SCM_COMPLEXP (y
))
4230 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4231 SCM_COMPLEX_IMAG (y
));
4232 else if (SCM_FRACTIONP (y
))
4233 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4235 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4237 else if (SCM_COMPLEXP (x
))
4239 if (SCM_I_INUMP (y
))
4240 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4241 SCM_COMPLEX_IMAG (x
));
4242 else if (SCM_BIGP (y
))
4244 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4245 + SCM_COMPLEX_REAL (x
));
4246 scm_remember_upto_here_1 (y
);
4247 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4249 else if (SCM_REALP (y
))
4250 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4251 SCM_COMPLEX_IMAG (x
));
4252 else if (SCM_COMPLEXP (y
))
4253 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4254 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4255 else if (SCM_FRACTIONP (y
))
4256 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4257 SCM_COMPLEX_IMAG (x
));
4259 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4261 else if (SCM_FRACTIONP (x
))
4263 if (SCM_I_INUMP (y
))
4264 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4265 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4266 SCM_FRACTION_DENOMINATOR (x
));
4267 else if (SCM_BIGP (y
))
4268 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4269 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4270 SCM_FRACTION_DENOMINATOR (x
));
4271 else if (SCM_REALP (y
))
4272 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4273 else if (SCM_COMPLEXP (y
))
4274 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4275 SCM_COMPLEX_IMAG (y
));
4276 else if (SCM_FRACTIONP (y
))
4277 /* a/b + c/d = (ad + bc) / bd */
4278 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4279 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4280 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4282 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4285 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4289 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4291 "Return @math{@var{x}+1}.")
4292 #define FUNC_NAME s_scm_oneplus
4294 return scm_sum (x
, SCM_I_MAKINUM (1));
4299 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4300 (SCM x
, SCM y
, SCM rest
),
4301 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4302 "the sum of all but the first argument are subtracted from the first\n"
4304 #define FUNC_NAME s_scm_i_difference
4306 while (!scm_is_null (rest
))
4307 { x
= scm_difference (x
, y
);
4309 rest
= scm_cdr (rest
);
4311 return scm_difference (x
, y
);
4315 #define s_difference s_scm_i_difference
4316 #define g_difference g_scm_i_difference
4319 scm_difference (SCM x
, SCM y
)
4320 #define FUNC_NAME s_difference
4322 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4325 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4327 if (SCM_I_INUMP (x
))
4329 long xx
= -SCM_I_INUM (x
);
4330 if (SCM_FIXABLE (xx
))
4331 return SCM_I_MAKINUM (xx
);
4333 return scm_i_long2big (xx
);
4335 else if (SCM_BIGP (x
))
4336 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4337 bignum, but negating that gives a fixnum. */
4338 return scm_i_normbig (scm_i_clonebig (x
, 0));
4339 else if (SCM_REALP (x
))
4340 return scm_from_double (-SCM_REAL_VALUE (x
));
4341 else if (SCM_COMPLEXP (x
))
4342 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4343 -SCM_COMPLEX_IMAG (x
));
4344 else if (SCM_FRACTIONP (x
))
4345 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4346 SCM_FRACTION_DENOMINATOR (x
));
4348 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4351 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4353 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4355 long int xx
= SCM_I_INUM (x
);
4356 long int yy
= SCM_I_INUM (y
);
4357 long int z
= xx
- yy
;
4358 if (SCM_FIXABLE (z
))
4359 return SCM_I_MAKINUM (z
);
4361 return scm_i_long2big (z
);
4363 else if (SCM_BIGP (y
))
4365 /* inum-x - big-y */
4366 long xx
= SCM_I_INUM (x
);
4369 return scm_i_clonebig (y
, 0);
4372 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4373 SCM result
= scm_i_mkbig ();
4376 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4379 /* x - y == -(y + -x) */
4380 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4381 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4383 scm_remember_upto_here_1 (y
);
4385 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4386 /* we know the result will have to be a bignum */
4389 return scm_i_normbig (result
);
4392 else if (SCM_REALP (y
))
4394 long int xx
= SCM_I_INUM (x
);
4395 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4397 else if (SCM_COMPLEXP (y
))
4399 long int xx
= SCM_I_INUM (x
);
4400 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4401 - SCM_COMPLEX_IMAG (y
));
4403 else if (SCM_FRACTIONP (y
))
4404 /* a - b/c = (ac - b) / c */
4405 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4406 SCM_FRACTION_NUMERATOR (y
)),
4407 SCM_FRACTION_DENOMINATOR (y
));
4409 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4411 else if (SCM_BIGP (x
))
4413 if (SCM_I_INUMP (y
))
4415 /* big-x - inum-y */
4416 long yy
= SCM_I_INUM (y
);
4417 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4419 scm_remember_upto_here_1 (x
);
4421 return (SCM_FIXABLE (-yy
) ?
4422 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4425 SCM result
= scm_i_mkbig ();
4428 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4430 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4431 scm_remember_upto_here_1 (x
);
4433 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4434 /* we know the result will have to be a bignum */
4437 return scm_i_normbig (result
);
4440 else if (SCM_BIGP (y
))
4442 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4443 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4444 SCM result
= scm_i_mkbig ();
4445 mpz_sub (SCM_I_BIG_MPZ (result
),
4448 scm_remember_upto_here_2 (x
, y
);
4449 /* we know the result will have to be a bignum */
4450 if ((sgn_x
== 1) && (sgn_y
== -1))
4452 if ((sgn_x
== -1) && (sgn_y
== 1))
4454 return scm_i_normbig (result
);
4456 else if (SCM_REALP (y
))
4458 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4459 scm_remember_upto_here_1 (x
);
4460 return scm_from_double (result
);
4462 else if (SCM_COMPLEXP (y
))
4464 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4465 - SCM_COMPLEX_REAL (y
));
4466 scm_remember_upto_here_1 (x
);
4467 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4469 else if (SCM_FRACTIONP (y
))
4470 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4471 SCM_FRACTION_NUMERATOR (y
)),
4472 SCM_FRACTION_DENOMINATOR (y
));
4473 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4475 else if (SCM_REALP (x
))
4477 if (SCM_I_INUMP (y
))
4478 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4479 else if (SCM_BIGP (y
))
4481 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4482 scm_remember_upto_here_1 (x
);
4483 return scm_from_double (result
);
4485 else if (SCM_REALP (y
))
4486 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4487 else if (SCM_COMPLEXP (y
))
4488 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4489 -SCM_COMPLEX_IMAG (y
));
4490 else if (SCM_FRACTIONP (y
))
4491 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4493 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4495 else if (SCM_COMPLEXP (x
))
4497 if (SCM_I_INUMP (y
))
4498 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4499 SCM_COMPLEX_IMAG (x
));
4500 else if (SCM_BIGP (y
))
4502 double real_part
= (SCM_COMPLEX_REAL (x
)
4503 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4504 scm_remember_upto_here_1 (x
);
4505 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4507 else if (SCM_REALP (y
))
4508 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4509 SCM_COMPLEX_IMAG (x
));
4510 else if (SCM_COMPLEXP (y
))
4511 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4512 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4513 else if (SCM_FRACTIONP (y
))
4514 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4515 SCM_COMPLEX_IMAG (x
));
4517 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4519 else if (SCM_FRACTIONP (x
))
4521 if (SCM_I_INUMP (y
))
4522 /* a/b - c = (a - cb) / b */
4523 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4524 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4525 SCM_FRACTION_DENOMINATOR (x
));
4526 else if (SCM_BIGP (y
))
4527 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4528 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4529 SCM_FRACTION_DENOMINATOR (x
));
4530 else if (SCM_REALP (y
))
4531 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4532 else if (SCM_COMPLEXP (y
))
4533 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4534 -SCM_COMPLEX_IMAG (y
));
4535 else if (SCM_FRACTIONP (y
))
4536 /* a/b - c/d = (ad - bc) / bd */
4537 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4538 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4539 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4541 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4544 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4549 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4551 "Return @math{@var{x}-1}.")
4552 #define FUNC_NAME s_scm_oneminus
4554 return scm_difference (x
, SCM_I_MAKINUM (1));
4559 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4560 (SCM x
, SCM y
, SCM rest
),
4561 "Return the product of all arguments. If called without arguments,\n"
4563 #define FUNC_NAME s_scm_i_product
4565 while (!scm_is_null (rest
))
4566 { x
= scm_product (x
, y
);
4568 rest
= scm_cdr (rest
);
4570 return scm_product (x
, y
);
4574 #define s_product s_scm_i_product
4575 #define g_product g_scm_i_product
4578 scm_product (SCM x
, SCM y
)
4580 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4583 return SCM_I_MAKINUM (1L);
4584 else if (SCM_NUMBERP (x
))
4587 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4590 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4595 xx
= SCM_I_INUM (x
);
4599 case 0: return x
; break;
4600 case 1: return y
; break;
4603 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4605 long yy
= SCM_I_INUM (y
);
4607 SCM k
= SCM_I_MAKINUM (kk
);
4608 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4612 SCM result
= scm_i_long2big (xx
);
4613 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4614 return scm_i_normbig (result
);
4617 else if (SCM_BIGP (y
))
4619 SCM result
= scm_i_mkbig ();
4620 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4621 scm_remember_upto_here_1 (y
);
4624 else if (SCM_REALP (y
))
4625 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4626 else if (SCM_COMPLEXP (y
))
4627 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4628 xx
* SCM_COMPLEX_IMAG (y
));
4629 else if (SCM_FRACTIONP (y
))
4630 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4631 SCM_FRACTION_DENOMINATOR (y
));
4633 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4635 else if (SCM_BIGP (x
))
4637 if (SCM_I_INUMP (y
))
4642 else if (SCM_BIGP (y
))
4644 SCM result
= scm_i_mkbig ();
4645 mpz_mul (SCM_I_BIG_MPZ (result
),
4648 scm_remember_upto_here_2 (x
, y
);
4651 else if (SCM_REALP (y
))
4653 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4654 scm_remember_upto_here_1 (x
);
4655 return scm_from_double (result
);
4657 else if (SCM_COMPLEXP (y
))
4659 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4660 scm_remember_upto_here_1 (x
);
4661 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4662 z
* SCM_COMPLEX_IMAG (y
));
4664 else if (SCM_FRACTIONP (y
))
4665 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4666 SCM_FRACTION_DENOMINATOR (y
));
4668 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4670 else if (SCM_REALP (x
))
4672 if (SCM_I_INUMP (y
))
4674 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4675 if (scm_is_eq (y
, SCM_INUM0
))
4677 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4679 else if (SCM_BIGP (y
))
4681 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4682 scm_remember_upto_here_1 (y
);
4683 return scm_from_double (result
);
4685 else if (SCM_REALP (y
))
4686 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4687 else if (SCM_COMPLEXP (y
))
4688 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4689 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4690 else if (SCM_FRACTIONP (y
))
4691 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4693 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4695 else if (SCM_COMPLEXP (x
))
4697 if (SCM_I_INUMP (y
))
4699 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4700 if (scm_is_eq (y
, SCM_INUM0
))
4702 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4703 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4705 else if (SCM_BIGP (y
))
4707 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4708 scm_remember_upto_here_1 (y
);
4709 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4710 z
* SCM_COMPLEX_IMAG (x
));
4712 else if (SCM_REALP (y
))
4713 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4714 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4715 else if (SCM_COMPLEXP (y
))
4717 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4718 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4719 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4720 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4722 else if (SCM_FRACTIONP (y
))
4724 double yy
= scm_i_fraction2double (y
);
4725 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4726 yy
* SCM_COMPLEX_IMAG (x
));
4729 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4731 else if (SCM_FRACTIONP (x
))
4733 if (SCM_I_INUMP (y
))
4734 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4735 SCM_FRACTION_DENOMINATOR (x
));
4736 else if (SCM_BIGP (y
))
4737 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4738 SCM_FRACTION_DENOMINATOR (x
));
4739 else if (SCM_REALP (y
))
4740 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4741 else if (SCM_COMPLEXP (y
))
4743 double xx
= scm_i_fraction2double (x
);
4744 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4745 xx
* SCM_COMPLEX_IMAG (y
));
4747 else if (SCM_FRACTIONP (y
))
4748 /* a/b * c/d = ac / bd */
4749 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4750 SCM_FRACTION_NUMERATOR (y
)),
4751 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4752 SCM_FRACTION_DENOMINATOR (y
)));
4754 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4757 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4760 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4761 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4762 #define ALLOW_DIVIDE_BY_ZERO
4763 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4766 /* The code below for complex division is adapted from the GNU
4767 libstdc++, which adapted it from f2c's libF77, and is subject to
4770 /****************************************************************
4771 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4773 Permission to use, copy, modify, and distribute this software
4774 and its documentation for any purpose and without fee is hereby
4775 granted, provided that the above copyright notice appear in all
4776 copies and that both that the copyright notice and this
4777 permission notice and warranty disclaimer appear in supporting
4778 documentation, and that the names of AT&T Bell Laboratories or
4779 Bellcore or any of their entities not be used in advertising or
4780 publicity pertaining to distribution of the software without
4781 specific, written prior permission.
4783 AT&T and Bellcore disclaim all warranties with regard to this
4784 software, including all implied warranties of merchantability
4785 and fitness. In no event shall AT&T or Bellcore be liable for
4786 any special, indirect or consequential damages or any damages
4787 whatsoever resulting from loss of use, data or profits, whether
4788 in an action of contract, negligence or other tortious action,
4789 arising out of or in connection with the use or performance of
4791 ****************************************************************/
4793 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4794 (SCM x
, SCM y
, SCM rest
),
4795 "Divide the first argument by the product of the remaining\n"
4796 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4798 #define FUNC_NAME s_scm_i_divide
4800 while (!scm_is_null (rest
))
4801 { x
= scm_divide (x
, y
);
4803 rest
= scm_cdr (rest
);
4805 return scm_divide (x
, y
);
4809 #define s_divide s_scm_i_divide
4810 #define g_divide g_scm_i_divide
4813 do_divide (SCM x
, SCM y
, int inexact
)
4814 #define FUNC_NAME s_divide
4818 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4821 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4822 else if (SCM_I_INUMP (x
))
4824 long xx
= SCM_I_INUM (x
);
4825 if (xx
== 1 || xx
== -1)
4827 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4829 scm_num_overflow (s_divide
);
4834 return scm_from_double (1.0 / (double) xx
);
4835 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4838 else if (SCM_BIGP (x
))
4841 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4842 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4844 else if (SCM_REALP (x
))
4846 double xx
= SCM_REAL_VALUE (x
);
4847 #ifndef ALLOW_DIVIDE_BY_ZERO
4849 scm_num_overflow (s_divide
);
4852 return scm_from_double (1.0 / xx
);
4854 else if (SCM_COMPLEXP (x
))
4856 double r
= SCM_COMPLEX_REAL (x
);
4857 double i
= SCM_COMPLEX_IMAG (x
);
4858 if (fabs(r
) <= fabs(i
))
4861 double d
= i
* (1.0 + t
* t
);
4862 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4867 double d
= r
* (1.0 + t
* t
);
4868 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4871 else if (SCM_FRACTIONP (x
))
4872 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4873 SCM_FRACTION_NUMERATOR (x
));
4875 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4878 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4880 long xx
= SCM_I_INUM (x
);
4881 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4883 long yy
= SCM_I_INUM (y
);
4886 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4887 scm_num_overflow (s_divide
);
4889 return scm_from_double ((double) xx
/ (double) yy
);
4892 else if (xx
% yy
!= 0)
4895 return scm_from_double ((double) xx
/ (double) yy
);
4896 else return scm_i_make_ratio (x
, y
);
4901 if (SCM_FIXABLE (z
))
4902 return SCM_I_MAKINUM (z
);
4904 return scm_i_long2big (z
);
4907 else if (SCM_BIGP (y
))
4910 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4911 else return scm_i_make_ratio (x
, y
);
4913 else if (SCM_REALP (y
))
4915 double yy
= SCM_REAL_VALUE (y
);
4916 #ifndef ALLOW_DIVIDE_BY_ZERO
4918 scm_num_overflow (s_divide
);
4921 return scm_from_double ((double) xx
/ yy
);
4923 else if (SCM_COMPLEXP (y
))
4926 complex_div
: /* y _must_ be a complex number */
4928 double r
= SCM_COMPLEX_REAL (y
);
4929 double i
= SCM_COMPLEX_IMAG (y
);
4930 if (fabs(r
) <= fabs(i
))
4933 double d
= i
* (1.0 + t
* t
);
4934 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4939 double d
= r
* (1.0 + t
* t
);
4940 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4944 else if (SCM_FRACTIONP (y
))
4945 /* a / b/c = ac / b */
4946 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4947 SCM_FRACTION_NUMERATOR (y
));
4949 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4951 else if (SCM_BIGP (x
))
4953 if (SCM_I_INUMP (y
))
4955 long int yy
= SCM_I_INUM (y
);
4958 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4959 scm_num_overflow (s_divide
);
4961 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4962 scm_remember_upto_here_1 (x
);
4963 return (sgn
== 0) ? scm_nan () : scm_inf ();
4970 /* FIXME: HMM, what are the relative performance issues here?
4971 We need to test. Is it faster on average to test
4972 divisible_p, then perform whichever operation, or is it
4973 faster to perform the integer div opportunistically and
4974 switch to real if there's a remainder? For now we take the
4975 middle ground: test, then if divisible, use the faster div
4978 long abs_yy
= yy
< 0 ? -yy
: yy
;
4979 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4983 SCM result
= scm_i_mkbig ();
4984 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4985 scm_remember_upto_here_1 (x
);
4987 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4988 return scm_i_normbig (result
);
4993 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4994 else return scm_i_make_ratio (x
, y
);
4998 else if (SCM_BIGP (y
))
5000 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5003 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5004 scm_num_overflow (s_divide
);
5006 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5007 scm_remember_upto_here_1 (x
);
5008 return (sgn
== 0) ? scm_nan () : scm_inf ();
5016 /* It's easily possible for the ratio x/y to fit a double
5017 but one or both x and y be too big to fit a double,
5018 hence the use of mpq_get_d rather than converting and
5021 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5022 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5023 return scm_from_double (mpq_get_d (q
));
5027 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5031 SCM result
= scm_i_mkbig ();
5032 mpz_divexact (SCM_I_BIG_MPZ (result
),
5035 scm_remember_upto_here_2 (x
, y
);
5036 return scm_i_normbig (result
);
5039 return scm_i_make_ratio (x
, y
);
5043 else if (SCM_REALP (y
))
5045 double yy
= SCM_REAL_VALUE (y
);
5046 #ifndef ALLOW_DIVIDE_BY_ZERO
5048 scm_num_overflow (s_divide
);
5051 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5053 else if (SCM_COMPLEXP (y
))
5055 a
= scm_i_big2dbl (x
);
5058 else if (SCM_FRACTIONP (y
))
5059 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5060 SCM_FRACTION_NUMERATOR (y
));
5062 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5064 else if (SCM_REALP (x
))
5066 double rx
= SCM_REAL_VALUE (x
);
5067 if (SCM_I_INUMP (y
))
5069 long int yy
= SCM_I_INUM (y
);
5070 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5072 scm_num_overflow (s_divide
);
5075 return scm_from_double (rx
/ (double) yy
);
5077 else if (SCM_BIGP (y
))
5079 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5080 scm_remember_upto_here_1 (y
);
5081 return scm_from_double (rx
/ dby
);
5083 else if (SCM_REALP (y
))
5085 double yy
= SCM_REAL_VALUE (y
);
5086 #ifndef ALLOW_DIVIDE_BY_ZERO
5088 scm_num_overflow (s_divide
);
5091 return scm_from_double (rx
/ yy
);
5093 else if (SCM_COMPLEXP (y
))
5098 else if (SCM_FRACTIONP (y
))
5099 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5101 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5103 else if (SCM_COMPLEXP (x
))
5105 double rx
= SCM_COMPLEX_REAL (x
);
5106 double ix
= SCM_COMPLEX_IMAG (x
);
5107 if (SCM_I_INUMP (y
))
5109 long int yy
= SCM_I_INUM (y
);
5110 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5112 scm_num_overflow (s_divide
);
5117 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5120 else if (SCM_BIGP (y
))
5122 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5123 scm_remember_upto_here_1 (y
);
5124 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5126 else if (SCM_REALP (y
))
5128 double yy
= SCM_REAL_VALUE (y
);
5129 #ifndef ALLOW_DIVIDE_BY_ZERO
5131 scm_num_overflow (s_divide
);
5134 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5136 else if (SCM_COMPLEXP (y
))
5138 double ry
= SCM_COMPLEX_REAL (y
);
5139 double iy
= SCM_COMPLEX_IMAG (y
);
5140 if (fabs(ry
) <= fabs(iy
))
5143 double d
= iy
* (1.0 + t
* t
);
5144 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5149 double d
= ry
* (1.0 + t
* t
);
5150 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5153 else if (SCM_FRACTIONP (y
))
5155 double yy
= scm_i_fraction2double (y
);
5156 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5159 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5161 else if (SCM_FRACTIONP (x
))
5163 if (SCM_I_INUMP (y
))
5165 long int yy
= SCM_I_INUM (y
);
5166 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5168 scm_num_overflow (s_divide
);
5171 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5172 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5174 else if (SCM_BIGP (y
))
5176 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5177 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5179 else if (SCM_REALP (y
))
5181 double yy
= SCM_REAL_VALUE (y
);
5182 #ifndef ALLOW_DIVIDE_BY_ZERO
5184 scm_num_overflow (s_divide
);
5187 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5189 else if (SCM_COMPLEXP (y
))
5191 a
= scm_i_fraction2double (x
);
5194 else if (SCM_FRACTIONP (y
))
5195 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5196 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5198 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5201 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5205 scm_divide (SCM x
, SCM y
)
5207 return do_divide (x
, y
, 0);
5210 static SCM
scm_divide2real (SCM x
, SCM y
)
5212 return do_divide (x
, y
, 1);
5218 scm_c_truncate (double x
)
5229 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5230 half-way case (ie. when x is an integer plus 0.5) going upwards.
5231 Then half-way cases are identified and adjusted down if the
5232 round-upwards didn't give the desired even integer.
5234 "plus_half == result" identifies a half-way case. If plus_half, which is
5235 x + 0.5, is an integer then x must be an integer plus 0.5.
5237 An odd "result" value is identified with result/2 != floor(result/2).
5238 This is done with plus_half, since that value is ready for use sooner in
5239 a pipelined cpu, and we're already requiring plus_half == result.
5241 Note however that we need to be careful when x is big and already an
5242 integer. In that case "x+0.5" may round to an adjacent integer, causing
5243 us to return such a value, incorrectly. For instance if the hardware is
5244 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5245 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5246 returned. Or if the hardware is in round-upwards mode, then other bigger
5247 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5248 representable value, 2^128+2^76 (or whatever), again incorrect.
5250 These bad roundings of x+0.5 are avoided by testing at the start whether
5251 x is already an integer. If it is then clearly that's the desired result
5252 already. And if it's not then the exponent must be small enough to allow
5253 an 0.5 to be represented, and hence added without a bad rounding. */
5256 scm_c_round (double x
)
5258 double plus_half
, result
;
5263 plus_half
= x
+ 0.5;
5264 result
= floor (plus_half
);
5265 /* Adjust so that the rounding is towards even. */
5266 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5271 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5273 "Round the number @var{x} towards zero.")
5274 #define FUNC_NAME s_scm_truncate_number
5276 if (scm_is_false (scm_negative_p (x
)))
5277 return scm_floor (x
);
5279 return scm_ceiling (x
);
5283 static SCM exactly_one_half
;
5285 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5287 "Round the number @var{x} towards the nearest integer. "
5288 "When it is exactly halfway between two integers, "
5289 "round towards the even one.")
5290 #define FUNC_NAME s_scm_round_number
5292 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5294 else if (SCM_REALP (x
))
5295 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5298 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5299 single quotient+remainder division then examining to see which way
5300 the rounding should go. */
5301 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5302 SCM result
= scm_floor (plus_half
);
5303 /* Adjust so that the rounding is towards even. */
5304 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5305 && scm_is_true (scm_odd_p (result
)))
5306 return scm_difference (result
, SCM_I_MAKINUM (1));
5313 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5315 "Round the number @var{x} towards minus infinity.")
5316 #define FUNC_NAME s_scm_floor
5318 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5320 else if (SCM_REALP (x
))
5321 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5322 else if (SCM_FRACTIONP (x
))
5324 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5325 SCM_FRACTION_DENOMINATOR (x
));
5326 if (scm_is_false (scm_negative_p (x
)))
5328 /* For positive x, rounding towards zero is correct. */
5333 /* For negative x, we need to return q-1 unless x is an
5334 integer. But fractions are never integer, per our
5336 return scm_difference (q
, SCM_I_MAKINUM (1));
5340 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5344 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5346 "Round the number @var{x} towards infinity.")
5347 #define FUNC_NAME s_scm_ceiling
5349 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5351 else if (SCM_REALP (x
))
5352 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5353 else if (SCM_FRACTIONP (x
))
5355 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5356 SCM_FRACTION_DENOMINATOR (x
));
5357 if (scm_is_false (scm_positive_p (x
)))
5359 /* For negative x, rounding towards zero is correct. */
5364 /* For positive x, we need to return q+1 unless x is an
5365 integer. But fractions are never integer, per our
5367 return scm_sum (q
, SCM_I_MAKINUM (1));
5371 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5375 /* sin/cos/tan/asin/acos/atan
5376 sinh/cosh/tanh/asinh/acosh/atanh
5377 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5378 Written by Jerry D. Hedden, (C) FSF.
5379 See the file `COPYING' for terms applying to this program. */
5381 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5383 "Return @var{x} raised to the power of @var{y}.")
5384 #define FUNC_NAME s_scm_expt
5386 if (!SCM_INEXACTP (y
) && scm_is_integer (y
))
5387 return scm_integer_expt (x
, y
);
5388 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5390 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5393 return scm_exp (scm_product (scm_log (x
), y
));
5397 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5399 "Compute the sine of @var{z}.")
5400 #define FUNC_NAME s_scm_sin
5402 if (scm_is_real (z
))
5403 return scm_from_double (sin (scm_to_double (z
)));
5404 else if (SCM_COMPLEXP (z
))
5406 x
= SCM_COMPLEX_REAL (z
);
5407 y
= SCM_COMPLEX_IMAG (z
);
5408 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5409 cos (x
) * sinh (y
));
5412 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5416 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5418 "Compute the cosine of @var{z}.")
5419 #define FUNC_NAME s_scm_cos
5421 if (scm_is_real (z
))
5422 return scm_from_double (cos (scm_to_double (z
)));
5423 else if (SCM_COMPLEXP (z
))
5425 x
= SCM_COMPLEX_REAL (z
);
5426 y
= SCM_COMPLEX_IMAG (z
);
5427 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5428 -sin (x
) * sinh (y
));
5431 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5435 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5437 "Compute the tangent of @var{z}.")
5438 #define FUNC_NAME s_scm_tan
5440 if (scm_is_real (z
))
5441 return scm_from_double (tan (scm_to_double (z
)));
5442 else if (SCM_COMPLEXP (z
))
5444 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5445 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5446 w
= cos (x
) + cosh (y
);
5447 #ifndef ALLOW_DIVIDE_BY_ZERO
5449 scm_num_overflow (s_scm_tan
);
5451 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5454 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5458 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5460 "Compute the hyperbolic sine of @var{z}.")
5461 #define FUNC_NAME s_scm_sinh
5463 if (scm_is_real (z
))
5464 return scm_from_double (sinh (scm_to_double (z
)));
5465 else if (SCM_COMPLEXP (z
))
5467 x
= SCM_COMPLEX_REAL (z
);
5468 y
= SCM_COMPLEX_IMAG (z
);
5469 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5470 cosh (x
) * sin (y
));
5473 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5477 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5479 "Compute the hyperbolic cosine of @var{z}.")
5480 #define FUNC_NAME s_scm_cosh
5482 if (scm_is_real (z
))
5483 return scm_from_double (cosh (scm_to_double (z
)));
5484 else if (SCM_COMPLEXP (z
))
5486 x
= SCM_COMPLEX_REAL (z
);
5487 y
= SCM_COMPLEX_IMAG (z
);
5488 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5489 sinh (x
) * sin (y
));
5492 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5496 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5498 "Compute the hyperbolic tangent of @var{z}.")
5499 #define FUNC_NAME s_scm_tanh
5501 if (scm_is_real (z
))
5502 return scm_from_double (tanh (scm_to_double (z
)));
5503 else if (SCM_COMPLEXP (z
))
5505 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5506 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5507 w
= cosh (x
) + cos (y
);
5508 #ifndef ALLOW_DIVIDE_BY_ZERO
5510 scm_num_overflow (s_scm_tanh
);
5512 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5515 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5519 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5521 "Compute the arc sine of @var{z}.")
5522 #define FUNC_NAME s_scm_asin
5524 if (scm_is_real (z
))
5526 double w
= scm_to_double (z
);
5527 if (w
>= -1.0 && w
<= 1.0)
5528 return scm_from_double (asin (w
));
5530 return scm_product (scm_c_make_rectangular (0, -1),
5531 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5533 else if (SCM_COMPLEXP (z
))
5535 x
= SCM_COMPLEX_REAL (z
);
5536 y
= SCM_COMPLEX_IMAG (z
);
5537 return scm_product (scm_c_make_rectangular (0, -1),
5538 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5541 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5545 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5547 "Compute the arc cosine of @var{z}.")
5548 #define FUNC_NAME s_scm_acos
5550 if (scm_is_real (z
))
5552 double w
= scm_to_double (z
);
5553 if (w
>= -1.0 && w
<= 1.0)
5554 return scm_from_double (acos (w
));
5556 return scm_sum (scm_from_double (acos (0.0)),
5557 scm_product (scm_c_make_rectangular (0, 1),
5558 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5560 else if (SCM_COMPLEXP (z
))
5562 x
= SCM_COMPLEX_REAL (z
);
5563 y
= SCM_COMPLEX_IMAG (z
);
5564 return scm_sum (scm_from_double (acos (0.0)),
5565 scm_product (scm_c_make_rectangular (0, 1),
5566 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5569 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5573 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5575 "With one argument, compute the arc tangent of @var{z}.\n"
5576 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5577 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5578 #define FUNC_NAME s_scm_atan
5582 if (scm_is_real (z
))
5583 return scm_from_double (atan (scm_to_double (z
)));
5584 else if (SCM_COMPLEXP (z
))
5587 v
= SCM_COMPLEX_REAL (z
);
5588 w
= SCM_COMPLEX_IMAG (z
);
5589 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5590 scm_c_make_rectangular (v
, w
+ 1.0))),
5591 scm_c_make_rectangular (0, 2));
5594 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5596 else if (scm_is_real (z
))
5598 if (scm_is_real (y
))
5599 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5601 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5604 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5608 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5610 "Compute the inverse hyperbolic sine of @var{z}.")
5611 #define FUNC_NAME s_scm_sys_asinh
5613 if (scm_is_real (z
))
5614 return scm_from_double (asinh (scm_to_double (z
)));
5615 else if (scm_is_number (z
))
5616 return scm_log (scm_sum (z
,
5617 scm_sqrt (scm_sum (scm_product (z
, z
),
5618 SCM_I_MAKINUM (1)))));
5620 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5624 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5626 "Compute the inverse hyperbolic cosine of @var{z}.")
5627 #define FUNC_NAME s_scm_sys_acosh
5629 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5630 return scm_from_double (acosh (scm_to_double (z
)));
5631 else if (scm_is_number (z
))
5632 return scm_log (scm_sum (z
,
5633 scm_sqrt (scm_difference (scm_product (z
, z
),
5634 SCM_I_MAKINUM (1)))));
5636 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5640 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5642 "Compute the inverse hyperbolic tangent of @var{z}.")
5643 #define FUNC_NAME s_scm_sys_atanh
5645 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5646 return scm_from_double (atanh (scm_to_double (z
)));
5647 else if (scm_is_number (z
))
5648 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5649 scm_difference (SCM_I_MAKINUM (1), z
))),
5652 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5657 scm_c_make_rectangular (double re
, double im
)
5660 return scm_from_double (re
);
5664 SCM_NEWSMOB (z
, scm_tc16_complex
,
5665 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5667 SCM_COMPLEX_REAL (z
) = re
;
5668 SCM_COMPLEX_IMAG (z
) = im
;
5673 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5674 (SCM real_part
, SCM imaginary_part
),
5675 "Return a complex number constructed of the given @var{real-part} "
5676 "and @var{imaginary-part} parts.")
5677 #define FUNC_NAME s_scm_make_rectangular
5679 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5680 SCM_ARG1
, FUNC_NAME
, "real");
5681 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5682 SCM_ARG2
, FUNC_NAME
, "real");
5683 return scm_c_make_rectangular (scm_to_double (real_part
),
5684 scm_to_double (imaginary_part
));
5689 scm_c_make_polar (double mag
, double ang
)
5693 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5694 use it on Glibc-based systems that have it (it's a GNU extension). See
5695 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5697 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5698 sincos (ang
, &s
, &c
);
5703 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5706 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5708 "Return the complex number @var{x} * e^(i * @var{y}).")
5709 #define FUNC_NAME s_scm_make_polar
5711 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5712 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5713 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5718 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5719 /* "Return the real part of the number @var{z}."
5722 scm_real_part (SCM z
)
5724 if (SCM_I_INUMP (z
))
5726 else if (SCM_BIGP (z
))
5728 else if (SCM_REALP (z
))
5730 else if (SCM_COMPLEXP (z
))
5731 return scm_from_double (SCM_COMPLEX_REAL (z
));
5732 else if (SCM_FRACTIONP (z
))
5735 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5739 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5740 /* "Return the imaginary part of the number @var{z}."
5743 scm_imag_part (SCM z
)
5745 if (SCM_I_INUMP (z
))
5747 else if (SCM_BIGP (z
))
5749 else if (SCM_REALP (z
))
5751 else if (SCM_COMPLEXP (z
))
5752 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5753 else if (SCM_FRACTIONP (z
))
5756 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5759 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5760 /* "Return the numerator of the number @var{z}."
5763 scm_numerator (SCM z
)
5765 if (SCM_I_INUMP (z
))
5767 else if (SCM_BIGP (z
))
5769 else if (SCM_FRACTIONP (z
))
5770 return SCM_FRACTION_NUMERATOR (z
);
5771 else if (SCM_REALP (z
))
5772 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5774 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5778 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5779 /* "Return the denominator of the number @var{z}."
5782 scm_denominator (SCM z
)
5784 if (SCM_I_INUMP (z
))
5785 return SCM_I_MAKINUM (1);
5786 else if (SCM_BIGP (z
))
5787 return SCM_I_MAKINUM (1);
5788 else if (SCM_FRACTIONP (z
))
5789 return SCM_FRACTION_DENOMINATOR (z
);
5790 else if (SCM_REALP (z
))
5791 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5793 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5796 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5797 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5798 * "@code{abs} for real arguments, but also allows complex numbers."
5801 scm_magnitude (SCM z
)
5803 if (SCM_I_INUMP (z
))
5805 long int zz
= SCM_I_INUM (z
);
5808 else if (SCM_POSFIXABLE (-zz
))
5809 return SCM_I_MAKINUM (-zz
);
5811 return scm_i_long2big (-zz
);
5813 else if (SCM_BIGP (z
))
5815 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5816 scm_remember_upto_here_1 (z
);
5818 return scm_i_clonebig (z
, 0);
5822 else if (SCM_REALP (z
))
5823 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5824 else if (SCM_COMPLEXP (z
))
5825 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5826 else if (SCM_FRACTIONP (z
))
5828 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5830 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5831 SCM_FRACTION_DENOMINATOR (z
));
5834 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5838 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5839 /* "Return the angle of the complex number @var{z}."
5844 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5845 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5846 But if atan2 follows the floating point rounding mode, then the value
5847 is not a constant. Maybe it'd be close enough though. */
5848 if (SCM_I_INUMP (z
))
5850 if (SCM_I_INUM (z
) >= 0)
5853 return scm_from_double (atan2 (0.0, -1.0));
5855 else if (SCM_BIGP (z
))
5857 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5858 scm_remember_upto_here_1 (z
);
5860 return scm_from_double (atan2 (0.0, -1.0));
5864 else if (SCM_REALP (z
))
5866 if (SCM_REAL_VALUE (z
) >= 0)
5869 return scm_from_double (atan2 (0.0, -1.0));
5871 else if (SCM_COMPLEXP (z
))
5872 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5873 else if (SCM_FRACTIONP (z
))
5875 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5877 else return scm_from_double (atan2 (0.0, -1.0));
5880 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5884 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5885 /* Convert the number @var{x} to its inexact representation.\n"
5888 scm_exact_to_inexact (SCM z
)
5890 if (SCM_I_INUMP (z
))
5891 return scm_from_double ((double) SCM_I_INUM (z
));
5892 else if (SCM_BIGP (z
))
5893 return scm_from_double (scm_i_big2dbl (z
));
5894 else if (SCM_FRACTIONP (z
))
5895 return scm_from_double (scm_i_fraction2double (z
));
5896 else if (SCM_INEXACTP (z
))
5899 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5903 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5905 "Return an exact number that is numerically closest to @var{z}.")
5906 #define FUNC_NAME s_scm_inexact_to_exact
5908 if (SCM_I_INUMP (z
))
5910 else if (SCM_BIGP (z
))
5912 else if (SCM_REALP (z
))
5914 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5915 SCM_OUT_OF_RANGE (1, z
);
5922 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5923 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5924 scm_i_mpz2num (mpq_denref (frac
)));
5926 /* When scm_i_make_ratio throws, we leak the memory allocated
5933 else if (SCM_FRACTIONP (z
))
5936 SCM_WRONG_TYPE_ARG (1, z
);
5940 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5942 "Returns the @emph{simplest} rational number differing\n"
5943 "from @var{x} by no more than @var{eps}.\n"
5945 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5946 "exact result when both its arguments are exact. Thus, you might need\n"
5947 "to use @code{inexact->exact} on the arguments.\n"
5950 "(rationalize (inexact->exact 1.2) 1/100)\n"
5953 #define FUNC_NAME s_scm_rationalize
5955 if (SCM_I_INUMP (x
))
5957 else if (SCM_BIGP (x
))
5959 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5961 /* Use continued fractions to find closest ratio. All
5962 arithmetic is done with exact numbers.
5965 SCM ex
= scm_inexact_to_exact (x
);
5966 SCM int_part
= scm_floor (ex
);
5967 SCM tt
= SCM_I_MAKINUM (1);
5968 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5969 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5973 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5976 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5977 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5979 /* We stop after a million iterations just to be absolutely sure
5980 that we don't go into an infinite loop. The process normally
5981 converges after less than a dozen iterations.
5984 eps
= scm_abs (eps
);
5985 while (++i
< 1000000)
5987 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5988 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5989 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5991 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5992 eps
))) /* abs(x-a/b) <= eps */
5994 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5995 if (scm_is_false (scm_exact_p (x
))
5996 || scm_is_false (scm_exact_p (eps
)))
5997 return scm_exact_to_inexact (res
);
6001 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6003 tt
= scm_floor (rx
); /* tt = floor (rx) */
6009 scm_num_overflow (s_scm_rationalize
);
6012 SCM_WRONG_TYPE_ARG (1, x
);
6016 /* conversion functions */
6019 scm_is_integer (SCM val
)
6021 return scm_is_true (scm_integer_p (val
));
6025 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6027 if (SCM_I_INUMP (val
))
6029 scm_t_signed_bits n
= SCM_I_INUM (val
);
6030 return n
>= min
&& n
<= max
;
6032 else if (SCM_BIGP (val
))
6034 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6036 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6038 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6040 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6041 return n
>= min
&& n
<= max
;
6051 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6052 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6055 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6056 SCM_I_BIG_MPZ (val
));
6058 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6070 return n
>= min
&& n
<= max
;
6078 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6080 if (SCM_I_INUMP (val
))
6082 scm_t_signed_bits n
= SCM_I_INUM (val
);
6083 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6085 else if (SCM_BIGP (val
))
6087 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6089 else if (max
<= ULONG_MAX
)
6091 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6093 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6094 return n
>= min
&& n
<= max
;
6104 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6107 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6108 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6111 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6112 SCM_I_BIG_MPZ (val
));
6114 return n
>= min
&& n
<= max
;
6122 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6124 scm_error (scm_out_of_range_key
,
6126 "Value out of range ~S to ~S: ~S",
6127 scm_list_3 (min
, max
, bad_val
),
6128 scm_list_1 (bad_val
));
6131 #define TYPE scm_t_intmax
6132 #define TYPE_MIN min
6133 #define TYPE_MAX max
6134 #define SIZEOF_TYPE 0
6135 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6136 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6137 #include "libguile/conv-integer.i.c"
6139 #define TYPE scm_t_uintmax
6140 #define TYPE_MIN min
6141 #define TYPE_MAX max
6142 #define SIZEOF_TYPE 0
6143 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6144 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6145 #include "libguile/conv-uinteger.i.c"
6147 #define TYPE scm_t_int8
6148 #define TYPE_MIN SCM_T_INT8_MIN
6149 #define TYPE_MAX SCM_T_INT8_MAX
6150 #define SIZEOF_TYPE 1
6151 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6152 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6153 #include "libguile/conv-integer.i.c"
6155 #define TYPE scm_t_uint8
6157 #define TYPE_MAX SCM_T_UINT8_MAX
6158 #define SIZEOF_TYPE 1
6159 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6160 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6161 #include "libguile/conv-uinteger.i.c"
6163 #define TYPE scm_t_int16
6164 #define TYPE_MIN SCM_T_INT16_MIN
6165 #define TYPE_MAX SCM_T_INT16_MAX
6166 #define SIZEOF_TYPE 2
6167 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6168 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6169 #include "libguile/conv-integer.i.c"
6171 #define TYPE scm_t_uint16
6173 #define TYPE_MAX SCM_T_UINT16_MAX
6174 #define SIZEOF_TYPE 2
6175 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6176 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6177 #include "libguile/conv-uinteger.i.c"
6179 #define TYPE scm_t_int32
6180 #define TYPE_MIN SCM_T_INT32_MIN
6181 #define TYPE_MAX SCM_T_INT32_MAX
6182 #define SIZEOF_TYPE 4
6183 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6184 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6185 #include "libguile/conv-integer.i.c"
6187 #define TYPE scm_t_uint32
6189 #define TYPE_MAX SCM_T_UINT32_MAX
6190 #define SIZEOF_TYPE 4
6191 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6192 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6193 #include "libguile/conv-uinteger.i.c"
6195 #define TYPE scm_t_wchar
6196 #define TYPE_MIN (scm_t_int32)-1
6197 #define TYPE_MAX (scm_t_int32)0x10ffff
6198 #define SIZEOF_TYPE 4
6199 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6200 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6201 #include "libguile/conv-integer.i.c"
6203 #if SCM_HAVE_T_INT64
6205 #define TYPE scm_t_int64
6206 #define TYPE_MIN SCM_T_INT64_MIN
6207 #define TYPE_MAX SCM_T_INT64_MAX
6208 #define SIZEOF_TYPE 8
6209 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6210 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6211 #include "libguile/conv-integer.i.c"
6213 #define TYPE scm_t_uint64
6215 #define TYPE_MAX SCM_T_UINT64_MAX
6216 #define SIZEOF_TYPE 8
6217 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6218 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6219 #include "libguile/conv-uinteger.i.c"
6224 scm_to_mpz (SCM val
, mpz_t rop
)
6226 if (SCM_I_INUMP (val
))
6227 mpz_set_si (rop
, SCM_I_INUM (val
));
6228 else if (SCM_BIGP (val
))
6229 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6231 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6235 scm_from_mpz (mpz_t val
)
6237 return scm_i_mpz2num (val
);
6241 scm_is_real (SCM val
)
6243 return scm_is_true (scm_real_p (val
));
6247 scm_is_rational (SCM val
)
6249 return scm_is_true (scm_rational_p (val
));
6253 scm_to_double (SCM val
)
6255 if (SCM_I_INUMP (val
))
6256 return SCM_I_INUM (val
);
6257 else if (SCM_BIGP (val
))
6258 return scm_i_big2dbl (val
);
6259 else if (SCM_FRACTIONP (val
))
6260 return scm_i_fraction2double (val
);
6261 else if (SCM_REALP (val
))
6262 return SCM_REAL_VALUE (val
);
6264 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6268 scm_from_double (double val
)
6270 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
6271 SCM_REAL_VALUE (z
) = val
;
6275 #if SCM_ENABLE_DISCOURAGED == 1
6278 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6282 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6286 scm_out_of_range (NULL
, num
);
6289 return scm_to_double (num
);
6293 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6297 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6301 scm_out_of_range (NULL
, num
);
6304 return scm_to_double (num
);
6310 scm_is_complex (SCM val
)
6312 return scm_is_true (scm_complex_p (val
));
6316 scm_c_real_part (SCM z
)
6318 if (SCM_COMPLEXP (z
))
6319 return SCM_COMPLEX_REAL (z
);
6322 /* Use the scm_real_part to get proper error checking and
6325 return scm_to_double (scm_real_part (z
));
6330 scm_c_imag_part (SCM z
)
6332 if (SCM_COMPLEXP (z
))
6333 return SCM_COMPLEX_IMAG (z
);
6336 /* Use the scm_imag_part to get proper error checking and
6337 dispatching. The result will almost always be 0.0, but not
6340 return scm_to_double (scm_imag_part (z
));
6345 scm_c_magnitude (SCM z
)
6347 return scm_to_double (scm_magnitude (z
));
6353 return scm_to_double (scm_angle (z
));
6357 scm_is_number (SCM z
)
6359 return scm_is_true (scm_number_p (z
));
6363 /* In the following functions we dispatch to the real-arg funcs like log()
6364 when we know the arg is real, instead of just handing everything to
6365 clog() for instance. This is in case clog() doesn't optimize for a
6366 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6367 well use it to go straight to the applicable C func. */
6369 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6371 "Return the natural logarithm of @var{z}.")
6372 #define FUNC_NAME s_scm_log
6374 if (SCM_COMPLEXP (z
))
6376 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6377 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6379 double re
= SCM_COMPLEX_REAL (z
);
6380 double im
= SCM_COMPLEX_IMAG (z
);
6381 return scm_c_make_rectangular (log (hypot (re
, im
)),
6387 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6388 although the value itself overflows. */
6389 double re
= scm_to_double (z
);
6390 double l
= log (fabs (re
));
6392 return scm_from_double (l
);
6394 return scm_c_make_rectangular (l
, M_PI
);
6400 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6402 "Return the base 10 logarithm of @var{z}.")
6403 #define FUNC_NAME s_scm_log10
6405 if (SCM_COMPLEXP (z
))
6407 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6408 clog() and a multiply by M_LOG10E, rather than the fallback
6409 log10+hypot+atan2.) */
6410 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6411 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6413 double re
= SCM_COMPLEX_REAL (z
);
6414 double im
= SCM_COMPLEX_IMAG (z
);
6415 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6416 M_LOG10E
* atan2 (im
, re
));
6421 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6422 although the value itself overflows. */
6423 double re
= scm_to_double (z
);
6424 double l
= log10 (fabs (re
));
6426 return scm_from_double (l
);
6428 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6434 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6436 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6437 "base of natural logarithms (2.71828@dots{}).")
6438 #define FUNC_NAME s_scm_exp
6440 if (SCM_COMPLEXP (z
))
6442 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6443 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6445 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6446 SCM_COMPLEX_IMAG (z
));
6451 /* When z is a negative bignum the conversion to double overflows,
6452 giving -infinity, but that's ok, the exp is still 0.0. */
6453 return scm_from_double (exp (scm_to_double (z
)));
6459 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6461 "Return the square root of @var{z}. Of the two possible roots\n"
6462 "(positive and negative), the one with the a positive real part\n"
6463 "is returned, or if that's zero then a positive imaginary part.\n"
6467 "(sqrt 9.0) @result{} 3.0\n"
6468 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6469 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6470 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6472 #define FUNC_NAME s_scm_sqrt
6474 if (SCM_COMPLEXP (x
))
6476 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6477 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6479 double re
= SCM_COMPLEX_REAL (x
);
6480 double im
= SCM_COMPLEX_IMAG (x
);
6481 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6482 0.5 * atan2 (im
, re
));
6487 double xx
= scm_to_double (x
);
6489 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6491 return scm_from_double (sqrt (xx
));
6503 mpz_init_set_si (z_negative_one
, -1);
6505 /* It may be possible to tune the performance of some algorithms by using
6506 * the following constants to avoid the creation of bignums. Please, before
6507 * using these values, remember the two rules of program optimization:
6508 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6509 scm_c_define ("most-positive-fixnum",
6510 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6511 scm_c_define ("most-negative-fixnum",
6512 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6514 scm_add_feature ("complex");
6515 scm_add_feature ("inexact");
6516 scm_flo0
= scm_from_double (0.0);
6518 /* determine floating point precision */
6519 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6521 init_dblprec(&scm_dblprec
[i
-2],i
);
6522 init_fx_radix(fx_per_radix
[i
-2],i
);
6525 /* hard code precision for base 10 if the preprocessor tells us to... */
6526 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6529 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6530 SCM_I_MAKINUM (2)));
6531 #include "libguile/numbers.x"