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
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3335 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3336 (SCM x
, SCM y
, SCM rest
),
3337 "Return @code{#t} if all parameters are numerically equal.")
3338 #define FUNC_NAME s_scm_i_num_eq_p
3340 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3342 while (!scm_is_null (rest
))
3344 if (scm_is_false (scm_num_eq_p (x
, y
)))
3348 rest
= scm_cdr (rest
);
3350 return scm_num_eq_p (x
, y
);
3354 scm_num_eq_p (SCM x
, SCM y
)
3357 if (SCM_I_INUMP (x
))
3359 long xx
= SCM_I_INUM (x
);
3360 if (SCM_I_INUMP (y
))
3362 long yy
= SCM_I_INUM (y
);
3363 return scm_from_bool (xx
== yy
);
3365 else if (SCM_BIGP (y
))
3367 else if (SCM_REALP (y
))
3369 /* On a 32-bit system an inum fits a double, we can cast the inum
3370 to a double and compare.
3372 But on a 64-bit system an inum is bigger than a double and
3373 casting it to a double (call that dxx) will round. dxx is at
3374 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3375 an integer and fits a long. So we cast yy to a long and
3376 compare with plain xx.
3378 An alternative (for any size system actually) would be to check
3379 yy is an integer (with floor) and is in range of an inum
3380 (compare against appropriate powers of 2) then test
3381 xx==(long)yy. It's just a matter of which casts/comparisons
3382 might be fastest or easiest for the cpu. */
3384 double yy
= SCM_REAL_VALUE (y
);
3385 return scm_from_bool ((double) xx
== yy
3386 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3387 || xx
== (long) yy
));
3389 else if (SCM_COMPLEXP (y
))
3390 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3391 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3392 else if (SCM_FRACTIONP (y
))
3395 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3397 else if (SCM_BIGP (x
))
3399 if (SCM_I_INUMP (y
))
3401 else if (SCM_BIGP (y
))
3403 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3404 scm_remember_upto_here_2 (x
, y
);
3405 return scm_from_bool (0 == cmp
);
3407 else if (SCM_REALP (y
))
3410 if (xisnan (SCM_REAL_VALUE (y
)))
3412 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3413 scm_remember_upto_here_1 (x
);
3414 return scm_from_bool (0 == cmp
);
3416 else if (SCM_COMPLEXP (y
))
3419 if (0.0 != SCM_COMPLEX_IMAG (y
))
3421 if (xisnan (SCM_COMPLEX_REAL (y
)))
3423 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3424 scm_remember_upto_here_1 (x
);
3425 return scm_from_bool (0 == cmp
);
3427 else if (SCM_FRACTIONP (y
))
3430 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3432 else if (SCM_REALP (x
))
3434 double xx
= SCM_REAL_VALUE (x
);
3435 if (SCM_I_INUMP (y
))
3437 /* see comments with inum/real above */
3438 long yy
= SCM_I_INUM (y
);
3439 return scm_from_bool (xx
== (double) yy
3440 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3441 || (long) xx
== yy
));
3443 else if (SCM_BIGP (y
))
3446 if (xisnan (SCM_REAL_VALUE (x
)))
3448 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3449 scm_remember_upto_here_1 (y
);
3450 return scm_from_bool (0 == cmp
);
3452 else if (SCM_REALP (y
))
3453 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3454 else if (SCM_COMPLEXP (y
))
3455 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3456 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3457 else if (SCM_FRACTIONP (y
))
3459 double xx
= SCM_REAL_VALUE (x
);
3463 return scm_from_bool (xx
< 0.0);
3464 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3468 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3470 else if (SCM_COMPLEXP (x
))
3472 if (SCM_I_INUMP (y
))
3473 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3474 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3475 else if (SCM_BIGP (y
))
3478 if (0.0 != SCM_COMPLEX_IMAG (x
))
3480 if (xisnan (SCM_COMPLEX_REAL (x
)))
3482 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3483 scm_remember_upto_here_1 (y
);
3484 return scm_from_bool (0 == cmp
);
3486 else if (SCM_REALP (y
))
3487 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3488 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3489 else if (SCM_COMPLEXP (y
))
3490 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3491 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3492 else if (SCM_FRACTIONP (y
))
3495 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3497 xx
= SCM_COMPLEX_REAL (x
);
3501 return scm_from_bool (xx
< 0.0);
3502 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3506 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3508 else if (SCM_FRACTIONP (x
))
3510 if (SCM_I_INUMP (y
))
3512 else if (SCM_BIGP (y
))
3514 else if (SCM_REALP (y
))
3516 double yy
= SCM_REAL_VALUE (y
);
3520 return scm_from_bool (0.0 < yy
);
3521 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3524 else if (SCM_COMPLEXP (y
))
3527 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3529 yy
= SCM_COMPLEX_REAL (y
);
3533 return scm_from_bool (0.0 < yy
);
3534 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3537 else if (SCM_FRACTIONP (y
))
3538 return scm_i_fraction_equalp (x
, y
);
3540 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3543 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3547 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3548 done are good for inums, but for bignums an answer can almost always be
3549 had by just examining a few high bits of the operands, as done by GMP in
3550 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3551 of the float exponent to take into account. */
3553 SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3554 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3555 (SCM x
, SCM y
, SCM rest
),
3556 "Return @code{#t} if the list of parameters is monotonically\n"
3558 #define FUNC_NAME s_scm_i_num_less_p
3560 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3562 while (!scm_is_null (rest
))
3564 if (scm_is_false (scm_less_p (x
, y
)))
3568 rest
= scm_cdr (rest
);
3570 return scm_less_p (x
, y
);
3574 scm_less_p (SCM x
, SCM y
)
3577 if (SCM_I_INUMP (x
))
3579 long xx
= SCM_I_INUM (x
);
3580 if (SCM_I_INUMP (y
))
3582 long yy
= SCM_I_INUM (y
);
3583 return scm_from_bool (xx
< yy
);
3585 else if (SCM_BIGP (y
))
3587 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3588 scm_remember_upto_here_1 (y
);
3589 return scm_from_bool (sgn
> 0);
3591 else if (SCM_REALP (y
))
3592 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3593 else if (SCM_FRACTIONP (y
))
3595 /* "x < a/b" becomes "x*b < a" */
3597 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3598 y
= SCM_FRACTION_NUMERATOR (y
);
3602 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3604 else if (SCM_BIGP (x
))
3606 if (SCM_I_INUMP (y
))
3608 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3609 scm_remember_upto_here_1 (x
);
3610 return scm_from_bool (sgn
< 0);
3612 else if (SCM_BIGP (y
))
3614 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3615 scm_remember_upto_here_2 (x
, y
);
3616 return scm_from_bool (cmp
< 0);
3618 else if (SCM_REALP (y
))
3621 if (xisnan (SCM_REAL_VALUE (y
)))
3623 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3624 scm_remember_upto_here_1 (x
);
3625 return scm_from_bool (cmp
< 0);
3627 else if (SCM_FRACTIONP (y
))
3630 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3632 else if (SCM_REALP (x
))
3634 if (SCM_I_INUMP (y
))
3635 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3636 else if (SCM_BIGP (y
))
3639 if (xisnan (SCM_REAL_VALUE (x
)))
3641 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3642 scm_remember_upto_here_1 (y
);
3643 return scm_from_bool (cmp
> 0);
3645 else if (SCM_REALP (y
))
3646 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3647 else if (SCM_FRACTIONP (y
))
3649 double xx
= SCM_REAL_VALUE (x
);
3653 return scm_from_bool (xx
< 0.0);
3654 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3658 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3660 else if (SCM_FRACTIONP (x
))
3662 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3664 /* "a/b < y" becomes "a < y*b" */
3665 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3666 x
= SCM_FRACTION_NUMERATOR (x
);
3669 else if (SCM_REALP (y
))
3671 double yy
= SCM_REAL_VALUE (y
);
3675 return scm_from_bool (0.0 < yy
);
3676 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3679 else if (SCM_FRACTIONP (y
))
3681 /* "a/b < c/d" becomes "a*d < c*b" */
3682 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3683 SCM_FRACTION_DENOMINATOR (y
));
3684 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3685 SCM_FRACTION_DENOMINATOR (x
));
3691 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3694 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3698 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3699 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3700 (SCM x
, SCM y
, SCM rest
),
3701 "Return @code{#t} if the list of parameters is monotonically\n"
3703 #define FUNC_NAME s_scm_i_num_gr_p
3705 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3707 while (!scm_is_null (rest
))
3709 if (scm_is_false (scm_gr_p (x
, y
)))
3713 rest
= scm_cdr (rest
);
3715 return scm_gr_p (x
, y
);
3718 #define FUNC_NAME s_scm_i_num_gr_p
3720 scm_gr_p (SCM x
, SCM y
)
3722 if (!SCM_NUMBERP (x
))
3723 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3724 else if (!SCM_NUMBERP (y
))
3725 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3727 return scm_less_p (y
, x
);
3732 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3733 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3734 (SCM x
, SCM y
, SCM rest
),
3735 "Return @code{#t} if the list of parameters is monotonically\n"
3737 #define FUNC_NAME s_scm_i_num_leq_p
3739 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3741 while (!scm_is_null (rest
))
3743 if (scm_is_false (scm_leq_p (x
, y
)))
3747 rest
= scm_cdr (rest
);
3749 return scm_leq_p (x
, y
);
3752 #define FUNC_NAME s_scm_i_num_leq_p
3754 scm_leq_p (SCM x
, SCM y
)
3756 if (!SCM_NUMBERP (x
))
3757 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3758 else if (!SCM_NUMBERP (y
))
3759 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3760 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3763 return scm_not (scm_less_p (y
, x
));
3768 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3769 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3770 (SCM x
, SCM y
, SCM rest
),
3771 "Return @code{#t} if the list of parameters is monotonically\n"
3773 #define FUNC_NAME s_scm_i_num_geq_p
3775 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3777 while (!scm_is_null (rest
))
3779 if (scm_is_false (scm_geq_p (x
, y
)))
3783 rest
= scm_cdr (rest
);
3785 return scm_geq_p (x
, y
);
3788 #define FUNC_NAME s_scm_i_num_geq_p
3790 scm_geq_p (SCM x
, SCM y
)
3792 if (!SCM_NUMBERP (x
))
3793 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3794 else if (!SCM_NUMBERP (y
))
3795 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3796 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3799 return scm_not (scm_less_p (x
, y
));
3804 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3805 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3811 if (SCM_I_INUMP (z
))
3812 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3813 else if (SCM_BIGP (z
))
3815 else if (SCM_REALP (z
))
3816 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3817 else if (SCM_COMPLEXP (z
))
3818 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3819 && SCM_COMPLEX_IMAG (z
) == 0.0);
3820 else if (SCM_FRACTIONP (z
))
3823 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3827 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3828 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3832 scm_positive_p (SCM x
)
3834 if (SCM_I_INUMP (x
))
3835 return scm_from_bool (SCM_I_INUM (x
) > 0);
3836 else if (SCM_BIGP (x
))
3838 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3839 scm_remember_upto_here_1 (x
);
3840 return scm_from_bool (sgn
> 0);
3842 else if (SCM_REALP (x
))
3843 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3844 else if (SCM_FRACTIONP (x
))
3845 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3847 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3851 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3852 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3856 scm_negative_p (SCM x
)
3858 if (SCM_I_INUMP (x
))
3859 return scm_from_bool (SCM_I_INUM (x
) < 0);
3860 else if (SCM_BIGP (x
))
3862 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3863 scm_remember_upto_here_1 (x
);
3864 return scm_from_bool (sgn
< 0);
3866 else if (SCM_REALP (x
))
3867 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3868 else if (SCM_FRACTIONP (x
))
3869 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3871 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3875 /* scm_min and scm_max return an inexact when either argument is inexact, as
3876 required by r5rs. On that basis, for exact/inexact combinations the
3877 exact is converted to inexact to compare and possibly return. This is
3878 unlike scm_less_p above which takes some trouble to preserve all bits in
3879 its test, such trouble is not required for min and max. */
3881 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3882 (SCM x
, SCM y
, SCM rest
),
3883 "Return the maximum of all parameter values.")
3884 #define FUNC_NAME s_scm_i_max
3886 while (!scm_is_null (rest
))
3887 { x
= scm_max (x
, y
);
3889 rest
= scm_cdr (rest
);
3891 return scm_max (x
, y
);
3895 #define s_max s_scm_i_max
3896 #define g_max g_scm_i_max
3899 scm_max (SCM x
, SCM y
)
3904 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3905 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3908 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3911 if (SCM_I_INUMP (x
))
3913 long xx
= SCM_I_INUM (x
);
3914 if (SCM_I_INUMP (y
))
3916 long yy
= SCM_I_INUM (y
);
3917 return (xx
< yy
) ? y
: x
;
3919 else if (SCM_BIGP (y
))
3921 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3922 scm_remember_upto_here_1 (y
);
3923 return (sgn
< 0) ? x
: y
;
3925 else if (SCM_REALP (y
))
3928 /* if y==NaN then ">" is false and we return NaN */
3929 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3931 else if (SCM_FRACTIONP (y
))
3934 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3937 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3939 else if (SCM_BIGP (x
))
3941 if (SCM_I_INUMP (y
))
3943 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3944 scm_remember_upto_here_1 (x
);
3945 return (sgn
< 0) ? y
: x
;
3947 else if (SCM_BIGP (y
))
3949 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3950 scm_remember_upto_here_2 (x
, y
);
3951 return (cmp
> 0) ? x
: y
;
3953 else if (SCM_REALP (y
))
3955 /* if y==NaN then xx>yy is false, so we return the NaN y */
3958 xx
= scm_i_big2dbl (x
);
3959 yy
= SCM_REAL_VALUE (y
);
3960 return (xx
> yy
? scm_from_double (xx
) : y
);
3962 else if (SCM_FRACTIONP (y
))
3967 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3969 else if (SCM_REALP (x
))
3971 if (SCM_I_INUMP (y
))
3973 double z
= SCM_I_INUM (y
);
3974 /* if x==NaN then "<" is false and we return NaN */
3975 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3977 else if (SCM_BIGP (y
))
3982 else if (SCM_REALP (y
))
3984 /* if x==NaN then our explicit check means we return NaN
3985 if y==NaN then ">" is false and we return NaN
3986 calling isnan is unavoidable, since it's the only way to know
3987 which of x or y causes any compares to be false */
3988 double xx
= SCM_REAL_VALUE (x
);
3989 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3991 else if (SCM_FRACTIONP (y
))
3993 double yy
= scm_i_fraction2double (y
);
3994 double xx
= SCM_REAL_VALUE (x
);
3995 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3998 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4000 else if (SCM_FRACTIONP (x
))
4002 if (SCM_I_INUMP (y
))
4006 else if (SCM_BIGP (y
))
4010 else if (SCM_REALP (y
))
4012 double xx
= scm_i_fraction2double (x
);
4013 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4015 else if (SCM_FRACTIONP (y
))
4020 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4023 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4027 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4028 (SCM x
, SCM y
, SCM rest
),
4029 "Return the minimum of all parameter values.")
4030 #define FUNC_NAME s_scm_i_min
4032 while (!scm_is_null (rest
))
4033 { x
= scm_min (x
, y
);
4035 rest
= scm_cdr (rest
);
4037 return scm_min (x
, y
);
4041 #define s_min s_scm_i_min
4042 #define g_min g_scm_i_min
4045 scm_min (SCM x
, SCM y
)
4050 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4051 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4054 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4057 if (SCM_I_INUMP (x
))
4059 long xx
= SCM_I_INUM (x
);
4060 if (SCM_I_INUMP (y
))
4062 long yy
= SCM_I_INUM (y
);
4063 return (xx
< yy
) ? x
: y
;
4065 else if (SCM_BIGP (y
))
4067 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4068 scm_remember_upto_here_1 (y
);
4069 return (sgn
< 0) ? y
: x
;
4071 else if (SCM_REALP (y
))
4074 /* if y==NaN then "<" is false and we return NaN */
4075 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4077 else if (SCM_FRACTIONP (y
))
4080 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4083 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4085 else if (SCM_BIGP (x
))
4087 if (SCM_I_INUMP (y
))
4089 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4090 scm_remember_upto_here_1 (x
);
4091 return (sgn
< 0) ? x
: y
;
4093 else if (SCM_BIGP (y
))
4095 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4096 scm_remember_upto_here_2 (x
, y
);
4097 return (cmp
> 0) ? y
: x
;
4099 else if (SCM_REALP (y
))
4101 /* if y==NaN then xx<yy is false, so we return the NaN y */
4104 xx
= scm_i_big2dbl (x
);
4105 yy
= SCM_REAL_VALUE (y
);
4106 return (xx
< yy
? scm_from_double (xx
) : y
);
4108 else if (SCM_FRACTIONP (y
))
4113 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4115 else if (SCM_REALP (x
))
4117 if (SCM_I_INUMP (y
))
4119 double z
= SCM_I_INUM (y
);
4120 /* if x==NaN then "<" is false and we return NaN */
4121 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4123 else if (SCM_BIGP (y
))
4128 else if (SCM_REALP (y
))
4130 /* if x==NaN then our explicit check means we return NaN
4131 if y==NaN then "<" is false and we return NaN
4132 calling isnan is unavoidable, since it's the only way to know
4133 which of x or y causes any compares to be false */
4134 double xx
= SCM_REAL_VALUE (x
);
4135 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4137 else if (SCM_FRACTIONP (y
))
4139 double yy
= scm_i_fraction2double (y
);
4140 double xx
= SCM_REAL_VALUE (x
);
4141 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4144 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4146 else if (SCM_FRACTIONP (x
))
4148 if (SCM_I_INUMP (y
))
4152 else if (SCM_BIGP (y
))
4156 else if (SCM_REALP (y
))
4158 double xx
= scm_i_fraction2double (x
);
4159 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4161 else if (SCM_FRACTIONP (y
))
4166 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4169 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4173 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4174 (SCM x
, SCM y
, SCM rest
),
4175 "Return the sum of all parameter values. Return 0 if called without\n"
4177 #define FUNC_NAME s_scm_i_sum
4179 while (!scm_is_null (rest
))
4180 { x
= scm_sum (x
, y
);
4182 rest
= scm_cdr (rest
);
4184 return scm_sum (x
, y
);
4188 #define s_sum s_scm_i_sum
4189 #define g_sum g_scm_i_sum
4192 scm_sum (SCM x
, SCM y
)
4194 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4196 if (SCM_NUMBERP (x
)) return x
;
4197 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4198 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4201 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4203 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4205 long xx
= SCM_I_INUM (x
);
4206 long yy
= SCM_I_INUM (y
);
4207 long int z
= xx
+ yy
;
4208 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4210 else if (SCM_BIGP (y
))
4215 else if (SCM_REALP (y
))
4217 long int xx
= SCM_I_INUM (x
);
4218 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4220 else if (SCM_COMPLEXP (y
))
4222 long int xx
= SCM_I_INUM (x
);
4223 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4224 SCM_COMPLEX_IMAG (y
));
4226 else if (SCM_FRACTIONP (y
))
4227 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4228 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4229 SCM_FRACTION_DENOMINATOR (y
));
4231 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4232 } else if (SCM_BIGP (x
))
4234 if (SCM_I_INUMP (y
))
4239 inum
= SCM_I_INUM (y
);
4242 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4245 SCM result
= scm_i_mkbig ();
4246 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4247 scm_remember_upto_here_1 (x
);
4248 /* we know the result will have to be a bignum */
4251 return scm_i_normbig (result
);
4255 SCM result
= scm_i_mkbig ();
4256 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4257 scm_remember_upto_here_1 (x
);
4258 /* we know the result will have to be a bignum */
4261 return scm_i_normbig (result
);
4264 else if (SCM_BIGP (y
))
4266 SCM result
= scm_i_mkbig ();
4267 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4268 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4269 mpz_add (SCM_I_BIG_MPZ (result
),
4272 scm_remember_upto_here_2 (x
, y
);
4273 /* we know the result will have to be a bignum */
4276 return scm_i_normbig (result
);
4278 else if (SCM_REALP (y
))
4280 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4281 scm_remember_upto_here_1 (x
);
4282 return scm_from_double (result
);
4284 else if (SCM_COMPLEXP (y
))
4286 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4287 + SCM_COMPLEX_REAL (y
));
4288 scm_remember_upto_here_1 (x
);
4289 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4291 else if (SCM_FRACTIONP (y
))
4292 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4293 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4294 SCM_FRACTION_DENOMINATOR (y
));
4296 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4298 else if (SCM_REALP (x
))
4300 if (SCM_I_INUMP (y
))
4301 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4302 else if (SCM_BIGP (y
))
4304 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4305 scm_remember_upto_here_1 (y
);
4306 return scm_from_double (result
);
4308 else if (SCM_REALP (y
))
4309 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4310 else if (SCM_COMPLEXP (y
))
4311 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4312 SCM_COMPLEX_IMAG (y
));
4313 else if (SCM_FRACTIONP (y
))
4314 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4316 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4318 else if (SCM_COMPLEXP (x
))
4320 if (SCM_I_INUMP (y
))
4321 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4322 SCM_COMPLEX_IMAG (x
));
4323 else if (SCM_BIGP (y
))
4325 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4326 + SCM_COMPLEX_REAL (x
));
4327 scm_remember_upto_here_1 (y
);
4328 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4330 else if (SCM_REALP (y
))
4331 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4332 SCM_COMPLEX_IMAG (x
));
4333 else if (SCM_COMPLEXP (y
))
4334 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4335 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4336 else if (SCM_FRACTIONP (y
))
4337 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4338 SCM_COMPLEX_IMAG (x
));
4340 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4342 else if (SCM_FRACTIONP (x
))
4344 if (SCM_I_INUMP (y
))
4345 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4346 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4347 SCM_FRACTION_DENOMINATOR (x
));
4348 else if (SCM_BIGP (y
))
4349 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4350 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4351 SCM_FRACTION_DENOMINATOR (x
));
4352 else if (SCM_REALP (y
))
4353 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4354 else if (SCM_COMPLEXP (y
))
4355 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4356 SCM_COMPLEX_IMAG (y
));
4357 else if (SCM_FRACTIONP (y
))
4358 /* a/b + c/d = (ad + bc) / bd */
4359 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4360 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4361 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4363 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4366 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4370 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4372 "Return @math{@var{x}+1}.")
4373 #define FUNC_NAME s_scm_oneplus
4375 return scm_sum (x
, SCM_I_MAKINUM (1));
4380 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4381 (SCM x
, SCM y
, SCM rest
),
4382 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4383 "the sum of all but the first argument are subtracted from the first\n"
4385 #define FUNC_NAME s_scm_i_difference
4387 while (!scm_is_null (rest
))
4388 { x
= scm_difference (x
, y
);
4390 rest
= scm_cdr (rest
);
4392 return scm_difference (x
, y
);
4396 #define s_difference s_scm_i_difference
4397 #define g_difference g_scm_i_difference
4400 scm_difference (SCM x
, SCM y
)
4401 #define FUNC_NAME s_difference
4403 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4406 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4408 if (SCM_I_INUMP (x
))
4410 long xx
= -SCM_I_INUM (x
);
4411 if (SCM_FIXABLE (xx
))
4412 return SCM_I_MAKINUM (xx
);
4414 return scm_i_long2big (xx
);
4416 else if (SCM_BIGP (x
))
4417 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4418 bignum, but negating that gives a fixnum. */
4419 return scm_i_normbig (scm_i_clonebig (x
, 0));
4420 else if (SCM_REALP (x
))
4421 return scm_from_double (-SCM_REAL_VALUE (x
));
4422 else if (SCM_COMPLEXP (x
))
4423 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4424 -SCM_COMPLEX_IMAG (x
));
4425 else if (SCM_FRACTIONP (x
))
4426 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4427 SCM_FRACTION_DENOMINATOR (x
));
4429 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4432 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4434 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4436 long int xx
= SCM_I_INUM (x
);
4437 long int yy
= SCM_I_INUM (y
);
4438 long int z
= xx
- yy
;
4439 if (SCM_FIXABLE (z
))
4440 return SCM_I_MAKINUM (z
);
4442 return scm_i_long2big (z
);
4444 else if (SCM_BIGP (y
))
4446 /* inum-x - big-y */
4447 long xx
= SCM_I_INUM (x
);
4450 return scm_i_clonebig (y
, 0);
4453 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4454 SCM result
= scm_i_mkbig ();
4457 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4460 /* x - y == -(y + -x) */
4461 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4462 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4464 scm_remember_upto_here_1 (y
);
4466 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4467 /* we know the result will have to be a bignum */
4470 return scm_i_normbig (result
);
4473 else if (SCM_REALP (y
))
4475 long int xx
= SCM_I_INUM (x
);
4476 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4478 else if (SCM_COMPLEXP (y
))
4480 long int xx
= SCM_I_INUM (x
);
4481 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4482 - SCM_COMPLEX_IMAG (y
));
4484 else if (SCM_FRACTIONP (y
))
4485 /* a - b/c = (ac - b) / c */
4486 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4487 SCM_FRACTION_NUMERATOR (y
)),
4488 SCM_FRACTION_DENOMINATOR (y
));
4490 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4492 else if (SCM_BIGP (x
))
4494 if (SCM_I_INUMP (y
))
4496 /* big-x - inum-y */
4497 long yy
= SCM_I_INUM (y
);
4498 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4500 scm_remember_upto_here_1 (x
);
4502 return (SCM_FIXABLE (-yy
) ?
4503 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4506 SCM result
= scm_i_mkbig ();
4509 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4511 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4512 scm_remember_upto_here_1 (x
);
4514 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4515 /* we know the result will have to be a bignum */
4518 return scm_i_normbig (result
);
4521 else if (SCM_BIGP (y
))
4523 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4524 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4525 SCM result
= scm_i_mkbig ();
4526 mpz_sub (SCM_I_BIG_MPZ (result
),
4529 scm_remember_upto_here_2 (x
, y
);
4530 /* we know the result will have to be a bignum */
4531 if ((sgn_x
== 1) && (sgn_y
== -1))
4533 if ((sgn_x
== -1) && (sgn_y
== 1))
4535 return scm_i_normbig (result
);
4537 else if (SCM_REALP (y
))
4539 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4540 scm_remember_upto_here_1 (x
);
4541 return scm_from_double (result
);
4543 else if (SCM_COMPLEXP (y
))
4545 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4546 - SCM_COMPLEX_REAL (y
));
4547 scm_remember_upto_here_1 (x
);
4548 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4550 else if (SCM_FRACTIONP (y
))
4551 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4552 SCM_FRACTION_NUMERATOR (y
)),
4553 SCM_FRACTION_DENOMINATOR (y
));
4554 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4556 else if (SCM_REALP (x
))
4558 if (SCM_I_INUMP (y
))
4559 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4560 else if (SCM_BIGP (y
))
4562 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4563 scm_remember_upto_here_1 (x
);
4564 return scm_from_double (result
);
4566 else if (SCM_REALP (y
))
4567 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4568 else if (SCM_COMPLEXP (y
))
4569 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4570 -SCM_COMPLEX_IMAG (y
));
4571 else if (SCM_FRACTIONP (y
))
4572 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4574 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4576 else if (SCM_COMPLEXP (x
))
4578 if (SCM_I_INUMP (y
))
4579 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4580 SCM_COMPLEX_IMAG (x
));
4581 else if (SCM_BIGP (y
))
4583 double real_part
= (SCM_COMPLEX_REAL (x
)
4584 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4585 scm_remember_upto_here_1 (x
);
4586 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4588 else if (SCM_REALP (y
))
4589 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4590 SCM_COMPLEX_IMAG (x
));
4591 else if (SCM_COMPLEXP (y
))
4592 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4593 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4594 else if (SCM_FRACTIONP (y
))
4595 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4596 SCM_COMPLEX_IMAG (x
));
4598 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4600 else if (SCM_FRACTIONP (x
))
4602 if (SCM_I_INUMP (y
))
4603 /* a/b - c = (a - cb) / b */
4604 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4605 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4606 SCM_FRACTION_DENOMINATOR (x
));
4607 else if (SCM_BIGP (y
))
4608 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4609 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4610 SCM_FRACTION_DENOMINATOR (x
));
4611 else if (SCM_REALP (y
))
4612 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4613 else if (SCM_COMPLEXP (y
))
4614 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4615 -SCM_COMPLEX_IMAG (y
));
4616 else if (SCM_FRACTIONP (y
))
4617 /* a/b - c/d = (ad - bc) / bd */
4618 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4619 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4620 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4622 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4625 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4630 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4632 "Return @math{@var{x}-1}.")
4633 #define FUNC_NAME s_scm_oneminus
4635 return scm_difference (x
, SCM_I_MAKINUM (1));
4640 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4641 (SCM x
, SCM y
, SCM rest
),
4642 "Return the product of all arguments. If called without arguments,\n"
4644 #define FUNC_NAME s_scm_i_product
4646 while (!scm_is_null (rest
))
4647 { x
= scm_product (x
, y
);
4649 rest
= scm_cdr (rest
);
4651 return scm_product (x
, y
);
4655 #define s_product s_scm_i_product
4656 #define g_product g_scm_i_product
4659 scm_product (SCM x
, SCM y
)
4661 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4664 return SCM_I_MAKINUM (1L);
4665 else if (SCM_NUMBERP (x
))
4668 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4671 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4676 xx
= SCM_I_INUM (x
);
4680 case 0: return x
; break;
4681 case 1: return y
; break;
4684 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4686 long yy
= SCM_I_INUM (y
);
4688 SCM k
= SCM_I_MAKINUM (kk
);
4689 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4693 SCM result
= scm_i_long2big (xx
);
4694 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4695 return scm_i_normbig (result
);
4698 else if (SCM_BIGP (y
))
4700 SCM result
= scm_i_mkbig ();
4701 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4702 scm_remember_upto_here_1 (y
);
4705 else if (SCM_REALP (y
))
4706 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4707 else if (SCM_COMPLEXP (y
))
4708 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4709 xx
* SCM_COMPLEX_IMAG (y
));
4710 else if (SCM_FRACTIONP (y
))
4711 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4712 SCM_FRACTION_DENOMINATOR (y
));
4714 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4716 else if (SCM_BIGP (x
))
4718 if (SCM_I_INUMP (y
))
4723 else if (SCM_BIGP (y
))
4725 SCM result
= scm_i_mkbig ();
4726 mpz_mul (SCM_I_BIG_MPZ (result
),
4729 scm_remember_upto_here_2 (x
, y
);
4732 else if (SCM_REALP (y
))
4734 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4735 scm_remember_upto_here_1 (x
);
4736 return scm_from_double (result
);
4738 else if (SCM_COMPLEXP (y
))
4740 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4741 scm_remember_upto_here_1 (x
);
4742 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4743 z
* SCM_COMPLEX_IMAG (y
));
4745 else if (SCM_FRACTIONP (y
))
4746 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4747 SCM_FRACTION_DENOMINATOR (y
));
4749 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4751 else if (SCM_REALP (x
))
4753 if (SCM_I_INUMP (y
))
4755 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4756 if (scm_is_eq (y
, SCM_INUM0
))
4758 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4760 else if (SCM_BIGP (y
))
4762 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4763 scm_remember_upto_here_1 (y
);
4764 return scm_from_double (result
);
4766 else if (SCM_REALP (y
))
4767 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4768 else if (SCM_COMPLEXP (y
))
4769 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4770 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4771 else if (SCM_FRACTIONP (y
))
4772 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4774 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4776 else if (SCM_COMPLEXP (x
))
4778 if (SCM_I_INUMP (y
))
4780 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4781 if (scm_is_eq (y
, SCM_INUM0
))
4783 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4784 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4786 else if (SCM_BIGP (y
))
4788 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4789 scm_remember_upto_here_1 (y
);
4790 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4791 z
* SCM_COMPLEX_IMAG (x
));
4793 else if (SCM_REALP (y
))
4794 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4795 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4796 else if (SCM_COMPLEXP (y
))
4798 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4799 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4800 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4801 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4803 else if (SCM_FRACTIONP (y
))
4805 double yy
= scm_i_fraction2double (y
);
4806 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4807 yy
* SCM_COMPLEX_IMAG (x
));
4810 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4812 else if (SCM_FRACTIONP (x
))
4814 if (SCM_I_INUMP (y
))
4815 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4816 SCM_FRACTION_DENOMINATOR (x
));
4817 else if (SCM_BIGP (y
))
4818 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4819 SCM_FRACTION_DENOMINATOR (x
));
4820 else if (SCM_REALP (y
))
4821 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4822 else if (SCM_COMPLEXP (y
))
4824 double xx
= scm_i_fraction2double (x
);
4825 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4826 xx
* SCM_COMPLEX_IMAG (y
));
4828 else if (SCM_FRACTIONP (y
))
4829 /* a/b * c/d = ac / bd */
4830 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4831 SCM_FRACTION_NUMERATOR (y
)),
4832 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4833 SCM_FRACTION_DENOMINATOR (y
)));
4835 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4838 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4841 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4842 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4843 #define ALLOW_DIVIDE_BY_ZERO
4844 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4847 /* The code below for complex division is adapted from the GNU
4848 libstdc++, which adapted it from f2c's libF77, and is subject to
4851 /****************************************************************
4852 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4854 Permission to use, copy, modify, and distribute this software
4855 and its documentation for any purpose and without fee is hereby
4856 granted, provided that the above copyright notice appear in all
4857 copies and that both that the copyright notice and this
4858 permission notice and warranty disclaimer appear in supporting
4859 documentation, and that the names of AT&T Bell Laboratories or
4860 Bellcore or any of their entities not be used in advertising or
4861 publicity pertaining to distribution of the software without
4862 specific, written prior permission.
4864 AT&T and Bellcore disclaim all warranties with regard to this
4865 software, including all implied warranties of merchantability
4866 and fitness. In no event shall AT&T or Bellcore be liable for
4867 any special, indirect or consequential damages or any damages
4868 whatsoever resulting from loss of use, data or profits, whether
4869 in an action of contract, negligence or other tortious action,
4870 arising out of or in connection with the use or performance of
4872 ****************************************************************/
4874 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4875 (SCM x
, SCM y
, SCM rest
),
4876 "Divide the first argument by the product of the remaining\n"
4877 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4879 #define FUNC_NAME s_scm_i_divide
4881 while (!scm_is_null (rest
))
4882 { x
= scm_divide (x
, y
);
4884 rest
= scm_cdr (rest
);
4886 return scm_divide (x
, y
);
4890 #define s_divide s_scm_i_divide
4891 #define g_divide g_scm_i_divide
4894 do_divide (SCM x
, SCM y
, int inexact
)
4895 #define FUNC_NAME s_divide
4899 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4902 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4903 else if (SCM_I_INUMP (x
))
4905 long xx
= SCM_I_INUM (x
);
4906 if (xx
== 1 || xx
== -1)
4908 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4910 scm_num_overflow (s_divide
);
4915 return scm_from_double (1.0 / (double) xx
);
4916 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4919 else if (SCM_BIGP (x
))
4922 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4923 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4925 else if (SCM_REALP (x
))
4927 double xx
= SCM_REAL_VALUE (x
);
4928 #ifndef ALLOW_DIVIDE_BY_ZERO
4930 scm_num_overflow (s_divide
);
4933 return scm_from_double (1.0 / xx
);
4935 else if (SCM_COMPLEXP (x
))
4937 double r
= SCM_COMPLEX_REAL (x
);
4938 double i
= SCM_COMPLEX_IMAG (x
);
4939 if (fabs(r
) <= fabs(i
))
4942 double d
= i
* (1.0 + t
* t
);
4943 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4948 double d
= r
* (1.0 + t
* t
);
4949 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4952 else if (SCM_FRACTIONP (x
))
4953 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4954 SCM_FRACTION_NUMERATOR (x
));
4956 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4959 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4961 long xx
= SCM_I_INUM (x
);
4962 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4964 long yy
= SCM_I_INUM (y
);
4967 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4968 scm_num_overflow (s_divide
);
4970 return scm_from_double ((double) xx
/ (double) yy
);
4973 else if (xx
% yy
!= 0)
4976 return scm_from_double ((double) xx
/ (double) yy
);
4977 else return scm_i_make_ratio (x
, y
);
4982 if (SCM_FIXABLE (z
))
4983 return SCM_I_MAKINUM (z
);
4985 return scm_i_long2big (z
);
4988 else if (SCM_BIGP (y
))
4991 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4992 else return scm_i_make_ratio (x
, y
);
4994 else if (SCM_REALP (y
))
4996 double yy
= SCM_REAL_VALUE (y
);
4997 #ifndef ALLOW_DIVIDE_BY_ZERO
4999 scm_num_overflow (s_divide
);
5002 return scm_from_double ((double) xx
/ yy
);
5004 else if (SCM_COMPLEXP (y
))
5007 complex_div
: /* y _must_ be a complex number */
5009 double r
= SCM_COMPLEX_REAL (y
);
5010 double i
= SCM_COMPLEX_IMAG (y
);
5011 if (fabs(r
) <= fabs(i
))
5014 double d
= i
* (1.0 + t
* t
);
5015 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5020 double d
= r
* (1.0 + t
* t
);
5021 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5025 else if (SCM_FRACTIONP (y
))
5026 /* a / b/c = ac / b */
5027 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5028 SCM_FRACTION_NUMERATOR (y
));
5030 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5032 else if (SCM_BIGP (x
))
5034 if (SCM_I_INUMP (y
))
5036 long int yy
= SCM_I_INUM (y
);
5039 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5040 scm_num_overflow (s_divide
);
5042 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5043 scm_remember_upto_here_1 (x
);
5044 return (sgn
== 0) ? scm_nan () : scm_inf ();
5051 /* FIXME: HMM, what are the relative performance issues here?
5052 We need to test. Is it faster on average to test
5053 divisible_p, then perform whichever operation, or is it
5054 faster to perform the integer div opportunistically and
5055 switch to real if there's a remainder? For now we take the
5056 middle ground: test, then if divisible, use the faster div
5059 long abs_yy
= yy
< 0 ? -yy
: yy
;
5060 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5064 SCM result
= scm_i_mkbig ();
5065 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5066 scm_remember_upto_here_1 (x
);
5068 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5069 return scm_i_normbig (result
);
5074 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5075 else return scm_i_make_ratio (x
, y
);
5079 else if (SCM_BIGP (y
))
5081 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5084 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5085 scm_num_overflow (s_divide
);
5087 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5088 scm_remember_upto_here_1 (x
);
5089 return (sgn
== 0) ? scm_nan () : scm_inf ();
5097 /* It's easily possible for the ratio x/y to fit a double
5098 but one or both x and y be too big to fit a double,
5099 hence the use of mpq_get_d rather than converting and
5102 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5103 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5104 return scm_from_double (mpq_get_d (q
));
5108 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5112 SCM result
= scm_i_mkbig ();
5113 mpz_divexact (SCM_I_BIG_MPZ (result
),
5116 scm_remember_upto_here_2 (x
, y
);
5117 return scm_i_normbig (result
);
5120 return scm_i_make_ratio (x
, y
);
5124 else if (SCM_REALP (y
))
5126 double yy
= SCM_REAL_VALUE (y
);
5127 #ifndef ALLOW_DIVIDE_BY_ZERO
5129 scm_num_overflow (s_divide
);
5132 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5134 else if (SCM_COMPLEXP (y
))
5136 a
= scm_i_big2dbl (x
);
5139 else if (SCM_FRACTIONP (y
))
5140 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5141 SCM_FRACTION_NUMERATOR (y
));
5143 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5145 else if (SCM_REALP (x
))
5147 double rx
= SCM_REAL_VALUE (x
);
5148 if (SCM_I_INUMP (y
))
5150 long int yy
= SCM_I_INUM (y
);
5151 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5153 scm_num_overflow (s_divide
);
5156 return scm_from_double (rx
/ (double) yy
);
5158 else if (SCM_BIGP (y
))
5160 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5161 scm_remember_upto_here_1 (y
);
5162 return scm_from_double (rx
/ dby
);
5164 else if (SCM_REALP (y
))
5166 double yy
= SCM_REAL_VALUE (y
);
5167 #ifndef ALLOW_DIVIDE_BY_ZERO
5169 scm_num_overflow (s_divide
);
5172 return scm_from_double (rx
/ yy
);
5174 else if (SCM_COMPLEXP (y
))
5179 else if (SCM_FRACTIONP (y
))
5180 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5182 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5184 else if (SCM_COMPLEXP (x
))
5186 double rx
= SCM_COMPLEX_REAL (x
);
5187 double ix
= SCM_COMPLEX_IMAG (x
);
5188 if (SCM_I_INUMP (y
))
5190 long int yy
= SCM_I_INUM (y
);
5191 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5193 scm_num_overflow (s_divide
);
5198 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5201 else if (SCM_BIGP (y
))
5203 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5204 scm_remember_upto_here_1 (y
);
5205 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5207 else if (SCM_REALP (y
))
5209 double yy
= SCM_REAL_VALUE (y
);
5210 #ifndef ALLOW_DIVIDE_BY_ZERO
5212 scm_num_overflow (s_divide
);
5215 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5217 else if (SCM_COMPLEXP (y
))
5219 double ry
= SCM_COMPLEX_REAL (y
);
5220 double iy
= SCM_COMPLEX_IMAG (y
);
5221 if (fabs(ry
) <= fabs(iy
))
5224 double d
= iy
* (1.0 + t
* t
);
5225 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5230 double d
= ry
* (1.0 + t
* t
);
5231 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5234 else if (SCM_FRACTIONP (y
))
5236 double yy
= scm_i_fraction2double (y
);
5237 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5240 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5242 else if (SCM_FRACTIONP (x
))
5244 if (SCM_I_INUMP (y
))
5246 long int yy
= SCM_I_INUM (y
);
5247 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5249 scm_num_overflow (s_divide
);
5252 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5253 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5255 else if (SCM_BIGP (y
))
5257 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5258 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5260 else if (SCM_REALP (y
))
5262 double yy
= SCM_REAL_VALUE (y
);
5263 #ifndef ALLOW_DIVIDE_BY_ZERO
5265 scm_num_overflow (s_divide
);
5268 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5270 else if (SCM_COMPLEXP (y
))
5272 a
= scm_i_fraction2double (x
);
5275 else if (SCM_FRACTIONP (y
))
5276 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5277 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5279 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5282 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5286 scm_divide (SCM x
, SCM y
)
5288 return do_divide (x
, y
, 0);
5291 static SCM
scm_divide2real (SCM x
, SCM y
)
5293 return do_divide (x
, y
, 1);
5299 scm_c_truncate (double x
)
5310 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5311 half-way case (ie. when x is an integer plus 0.5) going upwards.
5312 Then half-way cases are identified and adjusted down if the
5313 round-upwards didn't give the desired even integer.
5315 "plus_half == result" identifies a half-way case. If plus_half, which is
5316 x + 0.5, is an integer then x must be an integer plus 0.5.
5318 An odd "result" value is identified with result/2 != floor(result/2).
5319 This is done with plus_half, since that value is ready for use sooner in
5320 a pipelined cpu, and we're already requiring plus_half == result.
5322 Note however that we need to be careful when x is big and already an
5323 integer. In that case "x+0.5" may round to an adjacent integer, causing
5324 us to return such a value, incorrectly. For instance if the hardware is
5325 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5326 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5327 returned. Or if the hardware is in round-upwards mode, then other bigger
5328 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5329 representable value, 2^128+2^76 (or whatever), again incorrect.
5331 These bad roundings of x+0.5 are avoided by testing at the start whether
5332 x is already an integer. If it is then clearly that's the desired result
5333 already. And if it's not then the exponent must be small enough to allow
5334 an 0.5 to be represented, and hence added without a bad rounding. */
5337 scm_c_round (double x
)
5339 double plus_half
, result
;
5344 plus_half
= x
+ 0.5;
5345 result
= floor (plus_half
);
5346 /* Adjust so that the rounding is towards even. */
5347 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5352 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5354 "Round the number @var{x} towards zero.")
5355 #define FUNC_NAME s_scm_truncate_number
5357 if (scm_is_false (scm_negative_p (x
)))
5358 return scm_floor (x
);
5360 return scm_ceiling (x
);
5364 static SCM exactly_one_half
;
5366 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5368 "Round the number @var{x} towards the nearest integer. "
5369 "When it is exactly halfway between two integers, "
5370 "round towards the even one.")
5371 #define FUNC_NAME s_scm_round_number
5373 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5375 else if (SCM_REALP (x
))
5376 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5379 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5380 single quotient+remainder division then examining to see which way
5381 the rounding should go. */
5382 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5383 SCM result
= scm_floor (plus_half
);
5384 /* Adjust so that the rounding is towards even. */
5385 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5386 && scm_is_true (scm_odd_p (result
)))
5387 return scm_difference (result
, SCM_I_MAKINUM (1));
5394 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5396 "Round the number @var{x} towards minus infinity.")
5397 #define FUNC_NAME s_scm_floor
5399 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5401 else if (SCM_REALP (x
))
5402 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5403 else if (SCM_FRACTIONP (x
))
5405 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5406 SCM_FRACTION_DENOMINATOR (x
));
5407 if (scm_is_false (scm_negative_p (x
)))
5409 /* For positive x, rounding towards zero is correct. */
5414 /* For negative x, we need to return q-1 unless x is an
5415 integer. But fractions are never integer, per our
5417 return scm_difference (q
, SCM_I_MAKINUM (1));
5421 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5425 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5427 "Round the number @var{x} towards infinity.")
5428 #define FUNC_NAME s_scm_ceiling
5430 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5432 else if (SCM_REALP (x
))
5433 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5434 else if (SCM_FRACTIONP (x
))
5436 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5437 SCM_FRACTION_DENOMINATOR (x
));
5438 if (scm_is_false (scm_positive_p (x
)))
5440 /* For negative x, rounding towards zero is correct. */
5445 /* For positive x, we need to return q+1 unless x is an
5446 integer. But fractions are never integer, per our
5448 return scm_sum (q
, SCM_I_MAKINUM (1));
5452 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5456 /* sin/cos/tan/asin/acos/atan
5457 sinh/cosh/tanh/asinh/acosh/atanh
5458 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5459 Written by Jerry D. Hedden, (C) FSF.
5460 See the file `COPYING' for terms applying to this program. */
5462 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5464 "Return @var{x} raised to the power of @var{y}.")
5465 #define FUNC_NAME s_scm_expt
5467 if (!SCM_INEXACTP (y
) && scm_is_integer (y
))
5468 return scm_integer_expt (x
, y
);
5469 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5471 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5474 return scm_exp (scm_product (scm_log (x
), y
));
5478 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5480 "Compute the sine of @var{z}.")
5481 #define FUNC_NAME s_scm_sin
5483 if (scm_is_real (z
))
5484 return scm_from_double (sin (scm_to_double (z
)));
5485 else if (SCM_COMPLEXP (z
))
5487 x
= SCM_COMPLEX_REAL (z
);
5488 y
= SCM_COMPLEX_IMAG (z
);
5489 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5490 cos (x
) * sinh (y
));
5493 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5497 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5499 "Compute the cosine of @var{z}.")
5500 #define FUNC_NAME s_scm_cos
5502 if (scm_is_real (z
))
5503 return scm_from_double (cos (scm_to_double (z
)));
5504 else if (SCM_COMPLEXP (z
))
5506 x
= SCM_COMPLEX_REAL (z
);
5507 y
= SCM_COMPLEX_IMAG (z
);
5508 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5509 -sin (x
) * sinh (y
));
5512 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5516 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5518 "Compute the tangent of @var{z}.")
5519 #define FUNC_NAME s_scm_tan
5521 if (scm_is_real (z
))
5522 return scm_from_double (tan (scm_to_double (z
)));
5523 else if (SCM_COMPLEXP (z
))
5525 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5526 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5527 w
= cos (x
) + cosh (y
);
5528 #ifndef ALLOW_DIVIDE_BY_ZERO
5530 scm_num_overflow (s_scm_tan
);
5532 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5535 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5539 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5541 "Compute the hyperbolic sine of @var{z}.")
5542 #define FUNC_NAME s_scm_sinh
5544 if (scm_is_real (z
))
5545 return scm_from_double (sinh (scm_to_double (z
)));
5546 else if (SCM_COMPLEXP (z
))
5548 x
= SCM_COMPLEX_REAL (z
);
5549 y
= SCM_COMPLEX_IMAG (z
);
5550 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5551 cosh (x
) * sin (y
));
5554 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5558 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5560 "Compute the hyperbolic cosine of @var{z}.")
5561 #define FUNC_NAME s_scm_cosh
5563 if (scm_is_real (z
))
5564 return scm_from_double (cosh (scm_to_double (z
)));
5565 else if (SCM_COMPLEXP (z
))
5567 x
= SCM_COMPLEX_REAL (z
);
5568 y
= SCM_COMPLEX_IMAG (z
);
5569 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5570 sinh (x
) * sin (y
));
5573 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5577 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5579 "Compute the hyperbolic tangent of @var{z}.")
5580 #define FUNC_NAME s_scm_tanh
5582 if (scm_is_real (z
))
5583 return scm_from_double (tanh (scm_to_double (z
)));
5584 else if (SCM_COMPLEXP (z
))
5586 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5587 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5588 w
= cosh (x
) + cos (y
);
5589 #ifndef ALLOW_DIVIDE_BY_ZERO
5591 scm_num_overflow (s_scm_tanh
);
5593 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5596 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5600 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5602 "Compute the arc sine of @var{z}.")
5603 #define FUNC_NAME s_scm_asin
5605 if (scm_is_real (z
))
5607 double w
= scm_to_double (z
);
5608 if (w
>= -1.0 && w
<= 1.0)
5609 return scm_from_double (asin (w
));
5611 return scm_product (scm_c_make_rectangular (0, -1),
5612 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5614 else if (SCM_COMPLEXP (z
))
5616 x
= SCM_COMPLEX_REAL (z
);
5617 y
= SCM_COMPLEX_IMAG (z
);
5618 return scm_product (scm_c_make_rectangular (0, -1),
5619 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5622 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5626 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5628 "Compute the arc cosine of @var{z}.")
5629 #define FUNC_NAME s_scm_acos
5631 if (scm_is_real (z
))
5633 double w
= scm_to_double (z
);
5634 if (w
>= -1.0 && w
<= 1.0)
5635 return scm_from_double (acos (w
));
5637 return scm_sum (scm_from_double (acos (0.0)),
5638 scm_product (scm_c_make_rectangular (0, 1),
5639 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5641 else if (SCM_COMPLEXP (z
))
5643 x
= SCM_COMPLEX_REAL (z
);
5644 y
= SCM_COMPLEX_IMAG (z
);
5645 return scm_sum (scm_from_double (acos (0.0)),
5646 scm_product (scm_c_make_rectangular (0, 1),
5647 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5650 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5654 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5656 "With one argument, compute the arc tangent of @var{z}.\n"
5657 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5658 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5659 #define FUNC_NAME s_scm_atan
5663 if (scm_is_real (z
))
5664 return scm_from_double (atan (scm_to_double (z
)));
5665 else if (SCM_COMPLEXP (z
))
5668 v
= SCM_COMPLEX_REAL (z
);
5669 w
= SCM_COMPLEX_IMAG (z
);
5670 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5671 scm_c_make_rectangular (v
, w
+ 1.0))),
5672 scm_c_make_rectangular (0, 2));
5675 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5677 else if (scm_is_real (z
))
5679 if (scm_is_real (y
))
5680 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5682 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5685 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5689 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5691 "Compute the inverse hyperbolic sine of @var{z}.")
5692 #define FUNC_NAME s_scm_sys_asinh
5694 if (scm_is_real (z
))
5695 return scm_from_double (asinh (scm_to_double (z
)));
5696 else if (scm_is_number (z
))
5697 return scm_log (scm_sum (z
,
5698 scm_sqrt (scm_sum (scm_product (z
, z
),
5699 SCM_I_MAKINUM (1)))));
5701 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5705 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5707 "Compute the inverse hyperbolic cosine of @var{z}.")
5708 #define FUNC_NAME s_scm_sys_acosh
5710 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5711 return scm_from_double (acosh (scm_to_double (z
)));
5712 else if (scm_is_number (z
))
5713 return scm_log (scm_sum (z
,
5714 scm_sqrt (scm_difference (scm_product (z
, z
),
5715 SCM_I_MAKINUM (1)))));
5717 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5721 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5723 "Compute the inverse hyperbolic tangent of @var{z}.")
5724 #define FUNC_NAME s_scm_sys_atanh
5726 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5727 return scm_from_double (atanh (scm_to_double (z
)));
5728 else if (scm_is_number (z
))
5729 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5730 scm_difference (SCM_I_MAKINUM (1), z
))),
5733 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5738 scm_c_make_rectangular (double re
, double im
)
5741 return scm_from_double (re
);
5745 SCM_NEWSMOB (z
, scm_tc16_complex
,
5746 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5748 SCM_COMPLEX_REAL (z
) = re
;
5749 SCM_COMPLEX_IMAG (z
) = im
;
5754 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5755 (SCM real_part
, SCM imaginary_part
),
5756 "Return a complex number constructed of the given @var{real-part} "
5757 "and @var{imaginary-part} parts.")
5758 #define FUNC_NAME s_scm_make_rectangular
5760 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5761 SCM_ARG1
, FUNC_NAME
, "real");
5762 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5763 SCM_ARG2
, FUNC_NAME
, "real");
5764 return scm_c_make_rectangular (scm_to_double (real_part
),
5765 scm_to_double (imaginary_part
));
5770 scm_c_make_polar (double mag
, double ang
)
5774 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5775 use it on Glibc-based systems that have it (it's a GNU extension). See
5776 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5778 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5779 sincos (ang
, &s
, &c
);
5784 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5787 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5789 "Return the complex number @var{x} * e^(i * @var{y}).")
5790 #define FUNC_NAME s_scm_make_polar
5792 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5793 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5794 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5799 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5800 /* "Return the real part of the number @var{z}."
5803 scm_real_part (SCM z
)
5805 if (SCM_I_INUMP (z
))
5807 else if (SCM_BIGP (z
))
5809 else if (SCM_REALP (z
))
5811 else if (SCM_COMPLEXP (z
))
5812 return scm_from_double (SCM_COMPLEX_REAL (z
));
5813 else if (SCM_FRACTIONP (z
))
5816 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5820 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5821 /* "Return the imaginary part of the number @var{z}."
5824 scm_imag_part (SCM z
)
5826 if (SCM_I_INUMP (z
))
5828 else if (SCM_BIGP (z
))
5830 else if (SCM_REALP (z
))
5832 else if (SCM_COMPLEXP (z
))
5833 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5834 else if (SCM_FRACTIONP (z
))
5837 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5840 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5841 /* "Return the numerator of the number @var{z}."
5844 scm_numerator (SCM z
)
5846 if (SCM_I_INUMP (z
))
5848 else if (SCM_BIGP (z
))
5850 else if (SCM_FRACTIONP (z
))
5851 return SCM_FRACTION_NUMERATOR (z
);
5852 else if (SCM_REALP (z
))
5853 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5855 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5859 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5860 /* "Return the denominator of the number @var{z}."
5863 scm_denominator (SCM z
)
5865 if (SCM_I_INUMP (z
))
5866 return SCM_I_MAKINUM (1);
5867 else if (SCM_BIGP (z
))
5868 return SCM_I_MAKINUM (1);
5869 else if (SCM_FRACTIONP (z
))
5870 return SCM_FRACTION_DENOMINATOR (z
);
5871 else if (SCM_REALP (z
))
5872 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5874 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5877 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5878 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5879 * "@code{abs} for real arguments, but also allows complex numbers."
5882 scm_magnitude (SCM z
)
5884 if (SCM_I_INUMP (z
))
5886 long int zz
= SCM_I_INUM (z
);
5889 else if (SCM_POSFIXABLE (-zz
))
5890 return SCM_I_MAKINUM (-zz
);
5892 return scm_i_long2big (-zz
);
5894 else if (SCM_BIGP (z
))
5896 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5897 scm_remember_upto_here_1 (z
);
5899 return scm_i_clonebig (z
, 0);
5903 else if (SCM_REALP (z
))
5904 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5905 else if (SCM_COMPLEXP (z
))
5906 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5907 else if (SCM_FRACTIONP (z
))
5909 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5911 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5912 SCM_FRACTION_DENOMINATOR (z
));
5915 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5919 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5920 /* "Return the angle of the complex number @var{z}."
5925 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5926 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5927 But if atan2 follows the floating point rounding mode, then the value
5928 is not a constant. Maybe it'd be close enough though. */
5929 if (SCM_I_INUMP (z
))
5931 if (SCM_I_INUM (z
) >= 0)
5934 return scm_from_double (atan2 (0.0, -1.0));
5936 else if (SCM_BIGP (z
))
5938 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5939 scm_remember_upto_here_1 (z
);
5941 return scm_from_double (atan2 (0.0, -1.0));
5945 else if (SCM_REALP (z
))
5947 if (SCM_REAL_VALUE (z
) >= 0)
5950 return scm_from_double (atan2 (0.0, -1.0));
5952 else if (SCM_COMPLEXP (z
))
5953 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5954 else if (SCM_FRACTIONP (z
))
5956 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5958 else return scm_from_double (atan2 (0.0, -1.0));
5961 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5965 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5966 /* Convert the number @var{x} to its inexact representation.\n"
5969 scm_exact_to_inexact (SCM z
)
5971 if (SCM_I_INUMP (z
))
5972 return scm_from_double ((double) SCM_I_INUM (z
));
5973 else if (SCM_BIGP (z
))
5974 return scm_from_double (scm_i_big2dbl (z
));
5975 else if (SCM_FRACTIONP (z
))
5976 return scm_from_double (scm_i_fraction2double (z
));
5977 else if (SCM_INEXACTP (z
))
5980 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5984 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5986 "Return an exact number that is numerically closest to @var{z}.")
5987 #define FUNC_NAME s_scm_inexact_to_exact
5989 if (SCM_I_INUMP (z
))
5991 else if (SCM_BIGP (z
))
5993 else if (SCM_REALP (z
))
5995 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5996 SCM_OUT_OF_RANGE (1, z
);
6003 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6004 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6005 scm_i_mpz2num (mpq_denref (frac
)));
6007 /* When scm_i_make_ratio throws, we leak the memory allocated
6014 else if (SCM_FRACTIONP (z
))
6017 SCM_WRONG_TYPE_ARG (1, z
);
6021 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6023 "Returns the @emph{simplest} rational number differing\n"
6024 "from @var{x} by no more than @var{eps}.\n"
6026 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6027 "exact result when both its arguments are exact. Thus, you might need\n"
6028 "to use @code{inexact->exact} on the arguments.\n"
6031 "(rationalize (inexact->exact 1.2) 1/100)\n"
6034 #define FUNC_NAME s_scm_rationalize
6036 if (SCM_I_INUMP (x
))
6038 else if (SCM_BIGP (x
))
6040 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6042 /* Use continued fractions to find closest ratio. All
6043 arithmetic is done with exact numbers.
6046 SCM ex
= scm_inexact_to_exact (x
);
6047 SCM int_part
= scm_floor (ex
);
6048 SCM tt
= SCM_I_MAKINUM (1);
6049 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
6050 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
6054 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6057 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6058 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6060 /* We stop after a million iterations just to be absolutely sure
6061 that we don't go into an infinite loop. The process normally
6062 converges after less than a dozen iterations.
6065 eps
= scm_abs (eps
);
6066 while (++i
< 1000000)
6068 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6069 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6070 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6072 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6073 eps
))) /* abs(x-a/b) <= eps */
6075 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6076 if (scm_is_false (scm_exact_p (x
))
6077 || scm_is_false (scm_exact_p (eps
)))
6078 return scm_exact_to_inexact (res
);
6082 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6084 tt
= scm_floor (rx
); /* tt = floor (rx) */
6090 scm_num_overflow (s_scm_rationalize
);
6093 SCM_WRONG_TYPE_ARG (1, x
);
6097 /* conversion functions */
6100 scm_is_integer (SCM val
)
6102 return scm_is_true (scm_integer_p (val
));
6106 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6108 if (SCM_I_INUMP (val
))
6110 scm_t_signed_bits n
= SCM_I_INUM (val
);
6111 return n
>= min
&& n
<= max
;
6113 else if (SCM_BIGP (val
))
6115 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6117 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6119 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6121 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6122 return n
>= min
&& n
<= max
;
6132 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6133 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6136 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6137 SCM_I_BIG_MPZ (val
));
6139 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6151 return n
>= min
&& n
<= max
;
6159 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6161 if (SCM_I_INUMP (val
))
6163 scm_t_signed_bits n
= SCM_I_INUM (val
);
6164 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6166 else if (SCM_BIGP (val
))
6168 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6170 else if (max
<= ULONG_MAX
)
6172 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6174 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6175 return n
>= min
&& n
<= max
;
6185 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6188 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6189 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6192 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6193 SCM_I_BIG_MPZ (val
));
6195 return n
>= min
&& n
<= max
;
6203 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6205 scm_error (scm_out_of_range_key
,
6207 "Value out of range ~S to ~S: ~S",
6208 scm_list_3 (min
, max
, bad_val
),
6209 scm_list_1 (bad_val
));
6212 #define TYPE scm_t_intmax
6213 #define TYPE_MIN min
6214 #define TYPE_MAX max
6215 #define SIZEOF_TYPE 0
6216 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6217 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6218 #include "libguile/conv-integer.i.c"
6220 #define TYPE scm_t_uintmax
6221 #define TYPE_MIN min
6222 #define TYPE_MAX max
6223 #define SIZEOF_TYPE 0
6224 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6225 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6226 #include "libguile/conv-uinteger.i.c"
6228 #define TYPE scm_t_int8
6229 #define TYPE_MIN SCM_T_INT8_MIN
6230 #define TYPE_MAX SCM_T_INT8_MAX
6231 #define SIZEOF_TYPE 1
6232 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6233 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6234 #include "libguile/conv-integer.i.c"
6236 #define TYPE scm_t_uint8
6238 #define TYPE_MAX SCM_T_UINT8_MAX
6239 #define SIZEOF_TYPE 1
6240 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6241 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6242 #include "libguile/conv-uinteger.i.c"
6244 #define TYPE scm_t_int16
6245 #define TYPE_MIN SCM_T_INT16_MIN
6246 #define TYPE_MAX SCM_T_INT16_MAX
6247 #define SIZEOF_TYPE 2
6248 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6249 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6250 #include "libguile/conv-integer.i.c"
6252 #define TYPE scm_t_uint16
6254 #define TYPE_MAX SCM_T_UINT16_MAX
6255 #define SIZEOF_TYPE 2
6256 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6257 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6258 #include "libguile/conv-uinteger.i.c"
6260 #define TYPE scm_t_int32
6261 #define TYPE_MIN SCM_T_INT32_MIN
6262 #define TYPE_MAX SCM_T_INT32_MAX
6263 #define SIZEOF_TYPE 4
6264 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6265 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6266 #include "libguile/conv-integer.i.c"
6268 #define TYPE scm_t_uint32
6270 #define TYPE_MAX SCM_T_UINT32_MAX
6271 #define SIZEOF_TYPE 4
6272 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6273 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6274 #include "libguile/conv-uinteger.i.c"
6276 #define TYPE scm_t_wchar
6277 #define TYPE_MIN (scm_t_int32)-1
6278 #define TYPE_MAX (scm_t_int32)0x10ffff
6279 #define SIZEOF_TYPE 4
6280 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6281 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6282 #include "libguile/conv-integer.i.c"
6284 #if SCM_HAVE_T_INT64
6286 #define TYPE scm_t_int64
6287 #define TYPE_MIN SCM_T_INT64_MIN
6288 #define TYPE_MAX SCM_T_INT64_MAX
6289 #define SIZEOF_TYPE 8
6290 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6291 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6292 #include "libguile/conv-integer.i.c"
6294 #define TYPE scm_t_uint64
6296 #define TYPE_MAX SCM_T_UINT64_MAX
6297 #define SIZEOF_TYPE 8
6298 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6299 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6300 #include "libguile/conv-uinteger.i.c"
6305 scm_to_mpz (SCM val
, mpz_t rop
)
6307 if (SCM_I_INUMP (val
))
6308 mpz_set_si (rop
, SCM_I_INUM (val
));
6309 else if (SCM_BIGP (val
))
6310 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6312 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6316 scm_from_mpz (mpz_t val
)
6318 return scm_i_mpz2num (val
);
6322 scm_is_real (SCM val
)
6324 return scm_is_true (scm_real_p (val
));
6328 scm_is_rational (SCM val
)
6330 return scm_is_true (scm_rational_p (val
));
6334 scm_to_double (SCM val
)
6336 if (SCM_I_INUMP (val
))
6337 return SCM_I_INUM (val
);
6338 else if (SCM_BIGP (val
))
6339 return scm_i_big2dbl (val
);
6340 else if (SCM_FRACTIONP (val
))
6341 return scm_i_fraction2double (val
);
6342 else if (SCM_REALP (val
))
6343 return SCM_REAL_VALUE (val
);
6345 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6349 scm_from_double (double val
)
6351 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
6352 SCM_REAL_VALUE (z
) = val
;
6356 #if SCM_ENABLE_DISCOURAGED == 1
6359 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6363 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6367 scm_out_of_range (NULL
, num
);
6370 return scm_to_double (num
);
6374 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6378 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6382 scm_out_of_range (NULL
, num
);
6385 return scm_to_double (num
);
6391 scm_is_complex (SCM val
)
6393 return scm_is_true (scm_complex_p (val
));
6397 scm_c_real_part (SCM z
)
6399 if (SCM_COMPLEXP (z
))
6400 return SCM_COMPLEX_REAL (z
);
6403 /* Use the scm_real_part to get proper error checking and
6406 return scm_to_double (scm_real_part (z
));
6411 scm_c_imag_part (SCM z
)
6413 if (SCM_COMPLEXP (z
))
6414 return SCM_COMPLEX_IMAG (z
);
6417 /* Use the scm_imag_part to get proper error checking and
6418 dispatching. The result will almost always be 0.0, but not
6421 return scm_to_double (scm_imag_part (z
));
6426 scm_c_magnitude (SCM z
)
6428 return scm_to_double (scm_magnitude (z
));
6434 return scm_to_double (scm_angle (z
));
6438 scm_is_number (SCM z
)
6440 return scm_is_true (scm_number_p (z
));
6444 /* In the following functions we dispatch to the real-arg funcs like log()
6445 when we know the arg is real, instead of just handing everything to
6446 clog() for instance. This is in case clog() doesn't optimize for a
6447 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6448 well use it to go straight to the applicable C func. */
6450 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6452 "Return the natural logarithm of @var{z}.")
6453 #define FUNC_NAME s_scm_log
6455 if (SCM_COMPLEXP (z
))
6457 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6458 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6460 double re
= SCM_COMPLEX_REAL (z
);
6461 double im
= SCM_COMPLEX_IMAG (z
);
6462 return scm_c_make_rectangular (log (hypot (re
, im
)),
6468 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6469 although the value itself overflows. */
6470 double re
= scm_to_double (z
);
6471 double l
= log (fabs (re
));
6473 return scm_from_double (l
);
6475 return scm_c_make_rectangular (l
, M_PI
);
6481 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6483 "Return the base 10 logarithm of @var{z}.")
6484 #define FUNC_NAME s_scm_log10
6486 if (SCM_COMPLEXP (z
))
6488 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6489 clog() and a multiply by M_LOG10E, rather than the fallback
6490 log10+hypot+atan2.) */
6491 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6492 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6494 double re
= SCM_COMPLEX_REAL (z
);
6495 double im
= SCM_COMPLEX_IMAG (z
);
6496 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6497 M_LOG10E
* atan2 (im
, re
));
6502 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6503 although the value itself overflows. */
6504 double re
= scm_to_double (z
);
6505 double l
= log10 (fabs (re
));
6507 return scm_from_double (l
);
6509 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6515 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6517 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6518 "base of natural logarithms (2.71828@dots{}).")
6519 #define FUNC_NAME s_scm_exp
6521 if (SCM_COMPLEXP (z
))
6523 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6524 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6526 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6527 SCM_COMPLEX_IMAG (z
));
6532 /* When z is a negative bignum the conversion to double overflows,
6533 giving -infinity, but that's ok, the exp is still 0.0. */
6534 return scm_from_double (exp (scm_to_double (z
)));
6540 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6542 "Return the square root of @var{z}. Of the two possible roots\n"
6543 "(positive and negative), the one with the a positive real part\n"
6544 "is returned, or if that's zero then a positive imaginary part.\n"
6548 "(sqrt 9.0) @result{} 3.0\n"
6549 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6550 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6551 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6553 #define FUNC_NAME s_scm_sqrt
6555 if (SCM_COMPLEXP (x
))
6557 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6558 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6560 double re
= SCM_COMPLEX_REAL (x
);
6561 double im
= SCM_COMPLEX_IMAG (x
);
6562 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6563 0.5 * atan2 (im
, re
));
6568 double xx
= scm_to_double (x
);
6570 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6572 return scm_from_double (sqrt (xx
));
6584 mpz_init_set_si (z_negative_one
, -1);
6586 /* It may be possible to tune the performance of some algorithms by using
6587 * the following constants to avoid the creation of bignums. Please, before
6588 * using these values, remember the two rules of program optimization:
6589 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6590 scm_c_define ("most-positive-fixnum",
6591 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6592 scm_c_define ("most-negative-fixnum",
6593 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6595 scm_add_feature ("complex");
6596 scm_add_feature ("inexact");
6597 scm_flo0
= scm_from_double (0.0);
6599 /* determine floating point precision */
6600 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6602 init_dblprec(&scm_dblprec
[i
-2],i
);
6603 init_fx_radix(fx_per_radix
[i
-2],i
);
6606 /* hard code precision for base 10 if the preprocessor tells us to... */
6607 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6610 exactly_one_half
= scm_divide (SCM_I_MAKINUM (1), SCM_I_MAKINUM (2));
6611 #include "libguile/numbers.x"