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_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
1030 /* "Return the greatest common divisor of all arguments.\n"
1031 * "If called without arguments, 0 is returned."
1034 scm_gcd (SCM x
, SCM y
)
1037 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1039 if (SCM_I_INUMP (x
))
1041 if (SCM_I_INUMP (y
))
1043 long xx
= SCM_I_INUM (x
);
1044 long yy
= SCM_I_INUM (y
);
1045 long u
= xx
< 0 ? -xx
: xx
;
1046 long v
= yy
< 0 ? -yy
: yy
;
1056 /* Determine a common factor 2^k */
1057 while (!(1 & (u
| v
)))
1063 /* Now, any factor 2^n can be eliminated */
1083 return (SCM_POSFIXABLE (result
)
1084 ? SCM_I_MAKINUM (result
)
1085 : scm_i_long2big (result
));
1087 else if (SCM_BIGP (y
))
1093 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1095 else if (SCM_BIGP (x
))
1097 if (SCM_I_INUMP (y
))
1099 unsigned long result
;
1102 yy
= SCM_I_INUM (y
);
1107 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1108 scm_remember_upto_here_1 (x
);
1109 return (SCM_POSFIXABLE (result
)
1110 ? SCM_I_MAKINUM (result
)
1111 : scm_from_ulong (result
));
1113 else if (SCM_BIGP (y
))
1115 SCM result
= scm_i_mkbig ();
1116 mpz_gcd (SCM_I_BIG_MPZ (result
),
1119 scm_remember_upto_here_2 (x
, y
);
1120 return scm_i_normbig (result
);
1123 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1126 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1129 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1130 /* "Return the least common multiple of the arguments.\n"
1131 * "If called without arguments, 1 is returned."
1134 scm_lcm (SCM n1
, SCM n2
)
1136 if (SCM_UNBNDP (n2
))
1138 if (SCM_UNBNDP (n1
))
1139 return SCM_I_MAKINUM (1L);
1140 n2
= SCM_I_MAKINUM (1L);
1143 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1144 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1145 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1146 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1148 if (SCM_I_INUMP (n1
))
1150 if (SCM_I_INUMP (n2
))
1152 SCM d
= scm_gcd (n1
, n2
);
1153 if (scm_is_eq (d
, SCM_INUM0
))
1156 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1160 /* inum n1, big n2 */
1163 SCM result
= scm_i_mkbig ();
1164 long nn1
= SCM_I_INUM (n1
);
1165 if (nn1
== 0) return SCM_INUM0
;
1166 if (nn1
< 0) nn1
= - nn1
;
1167 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1168 scm_remember_upto_here_1 (n2
);
1176 if (SCM_I_INUMP (n2
))
1183 SCM result
= scm_i_mkbig ();
1184 mpz_lcm(SCM_I_BIG_MPZ (result
),
1186 SCM_I_BIG_MPZ (n2
));
1187 scm_remember_upto_here_2(n1
, n2
);
1188 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1194 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1199 + + + x (map digit:logand X Y)
1200 + - + x (map digit:logand X (lognot (+ -1 Y)))
1201 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1202 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1207 + + + (map digit:logior X Y)
1208 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1209 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1210 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1215 + + + (map digit:logxor X Y)
1216 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1217 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1218 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1223 + + (any digit:logand X Y)
1224 + - (any digit:logand X (lognot (+ -1 Y)))
1225 - + (any digit:logand (lognot (+ -1 X)) Y)
1230 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1232 "Return the bitwise AND of the integer arguments.\n\n"
1234 "(logand) @result{} -1\n"
1235 "(logand 7) @result{} 7\n"
1236 "(logand #b111 #b011 #b001) @result{} 1\n"
1238 #define FUNC_NAME s_scm_logand
1242 if (SCM_UNBNDP (n2
))
1244 if (SCM_UNBNDP (n1
))
1245 return SCM_I_MAKINUM (-1);
1246 else if (!SCM_NUMBERP (n1
))
1247 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1248 else if (SCM_NUMBERP (n1
))
1251 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1254 if (SCM_I_INUMP (n1
))
1256 nn1
= SCM_I_INUM (n1
);
1257 if (SCM_I_INUMP (n2
))
1259 long nn2
= SCM_I_INUM (n2
);
1260 return SCM_I_MAKINUM (nn1
& nn2
);
1262 else if SCM_BIGP (n2
)
1268 SCM result_z
= scm_i_mkbig ();
1270 mpz_init_set_si (nn1_z
, nn1
);
1271 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1272 scm_remember_upto_here_1 (n2
);
1274 return scm_i_normbig (result_z
);
1278 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1280 else if (SCM_BIGP (n1
))
1282 if (SCM_I_INUMP (n2
))
1285 nn1
= SCM_I_INUM (n1
);
1288 else if (SCM_BIGP (n2
))
1290 SCM result_z
= scm_i_mkbig ();
1291 mpz_and (SCM_I_BIG_MPZ (result_z
),
1293 SCM_I_BIG_MPZ (n2
));
1294 scm_remember_upto_here_2 (n1
, n2
);
1295 return scm_i_normbig (result_z
);
1298 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1301 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1306 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1308 "Return the bitwise OR of the integer arguments.\n\n"
1310 "(logior) @result{} 0\n"
1311 "(logior 7) @result{} 7\n"
1312 "(logior #b000 #b001 #b011) @result{} 3\n"
1314 #define FUNC_NAME s_scm_logior
1318 if (SCM_UNBNDP (n2
))
1320 if (SCM_UNBNDP (n1
))
1322 else if (SCM_NUMBERP (n1
))
1325 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1328 if (SCM_I_INUMP (n1
))
1330 nn1
= SCM_I_INUM (n1
);
1331 if (SCM_I_INUMP (n2
))
1333 long nn2
= SCM_I_INUM (n2
);
1334 return SCM_I_MAKINUM (nn1
| nn2
);
1336 else if (SCM_BIGP (n2
))
1342 SCM result_z
= scm_i_mkbig ();
1344 mpz_init_set_si (nn1_z
, nn1
);
1345 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1346 scm_remember_upto_here_1 (n2
);
1348 return scm_i_normbig (result_z
);
1352 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1354 else if (SCM_BIGP (n1
))
1356 if (SCM_I_INUMP (n2
))
1359 nn1
= SCM_I_INUM (n1
);
1362 else if (SCM_BIGP (n2
))
1364 SCM result_z
= scm_i_mkbig ();
1365 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1367 SCM_I_BIG_MPZ (n2
));
1368 scm_remember_upto_here_2 (n1
, n2
);
1369 return scm_i_normbig (result_z
);
1372 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1375 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1380 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1382 "Return the bitwise XOR of the integer arguments. A bit is\n"
1383 "set in the result if it is set in an odd number of arguments.\n"
1385 "(logxor) @result{} 0\n"
1386 "(logxor 7) @result{} 7\n"
1387 "(logxor #b000 #b001 #b011) @result{} 2\n"
1388 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1390 #define FUNC_NAME s_scm_logxor
1394 if (SCM_UNBNDP (n2
))
1396 if (SCM_UNBNDP (n1
))
1398 else if (SCM_NUMBERP (n1
))
1401 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1404 if (SCM_I_INUMP (n1
))
1406 nn1
= SCM_I_INUM (n1
);
1407 if (SCM_I_INUMP (n2
))
1409 long nn2
= SCM_I_INUM (n2
);
1410 return SCM_I_MAKINUM (nn1
^ nn2
);
1412 else if (SCM_BIGP (n2
))
1416 SCM result_z
= scm_i_mkbig ();
1418 mpz_init_set_si (nn1_z
, nn1
);
1419 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1420 scm_remember_upto_here_1 (n2
);
1422 return scm_i_normbig (result_z
);
1426 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1428 else if (SCM_BIGP (n1
))
1430 if (SCM_I_INUMP (n2
))
1433 nn1
= SCM_I_INUM (n1
);
1436 else if (SCM_BIGP (n2
))
1438 SCM result_z
= scm_i_mkbig ();
1439 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1441 SCM_I_BIG_MPZ (n2
));
1442 scm_remember_upto_here_2 (n1
, n2
);
1443 return scm_i_normbig (result_z
);
1446 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1449 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1454 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1456 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1457 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1458 "without actually calculating the @code{logand}, just testing\n"
1462 "(logtest #b0100 #b1011) @result{} #f\n"
1463 "(logtest #b0100 #b0111) @result{} #t\n"
1465 #define FUNC_NAME s_scm_logtest
1469 if (SCM_I_INUMP (j
))
1471 nj
= SCM_I_INUM (j
);
1472 if (SCM_I_INUMP (k
))
1474 long nk
= SCM_I_INUM (k
);
1475 return scm_from_bool (nj
& nk
);
1477 else if (SCM_BIGP (k
))
1485 mpz_init_set_si (nj_z
, nj
);
1486 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1487 scm_remember_upto_here_1 (k
);
1488 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1494 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1496 else if (SCM_BIGP (j
))
1498 if (SCM_I_INUMP (k
))
1501 nj
= SCM_I_INUM (j
);
1504 else if (SCM_BIGP (k
))
1508 mpz_init (result_z
);
1512 scm_remember_upto_here_2 (j
, k
);
1513 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1514 mpz_clear (result_z
);
1518 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1521 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1526 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1528 "Test whether bit number @var{index} in @var{j} is set.\n"
1529 "@var{index} starts from 0 for the least significant bit.\n"
1532 "(logbit? 0 #b1101) @result{} #t\n"
1533 "(logbit? 1 #b1101) @result{} #f\n"
1534 "(logbit? 2 #b1101) @result{} #t\n"
1535 "(logbit? 3 #b1101) @result{} #t\n"
1536 "(logbit? 4 #b1101) @result{} #f\n"
1538 #define FUNC_NAME s_scm_logbit_p
1540 unsigned long int iindex
;
1541 iindex
= scm_to_ulong (index
);
1543 if (SCM_I_INUMP (j
))
1545 /* bits above what's in an inum follow the sign bit */
1546 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1547 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1549 else if (SCM_BIGP (j
))
1551 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1552 scm_remember_upto_here_1 (j
);
1553 return scm_from_bool (val
);
1556 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1561 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1563 "Return the integer which is the ones-complement of the integer\n"
1567 "(number->string (lognot #b10000000) 2)\n"
1568 " @result{} \"-10000001\"\n"
1569 "(number->string (lognot #b0) 2)\n"
1570 " @result{} \"-1\"\n"
1572 #define FUNC_NAME s_scm_lognot
1574 if (SCM_I_INUMP (n
)) {
1575 /* No overflow here, just need to toggle all the bits making up the inum.
1576 Enhancement: No need to strip the tag and add it back, could just xor
1577 a block of 1 bits, if that worked with the various debug versions of
1579 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1581 } else if (SCM_BIGP (n
)) {
1582 SCM result
= scm_i_mkbig ();
1583 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1584 scm_remember_upto_here_1 (n
);
1588 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1593 /* returns 0 if IN is not an integer. OUT must already be
1596 coerce_to_big (SCM in
, mpz_t out
)
1599 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1600 else if (SCM_I_INUMP (in
))
1601 mpz_set_si (out
, SCM_I_INUM (in
));
1608 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1609 (SCM n
, SCM k
, SCM m
),
1610 "Return @var{n} raised to the integer exponent\n"
1611 "@var{k}, modulo @var{m}.\n"
1614 "(modulo-expt 2 3 5)\n"
1617 #define FUNC_NAME s_scm_modulo_expt
1623 /* There are two classes of error we might encounter --
1624 1) Math errors, which we'll report by calling scm_num_overflow,
1626 2) wrong-type errors, which of course we'll report by calling
1628 We don't report those errors immediately, however; instead we do
1629 some cleanup first. These variables tell us which error (if
1630 any) we should report after cleaning up.
1632 int report_overflow
= 0;
1634 int position_of_wrong_type
= 0;
1635 SCM value_of_wrong_type
= SCM_INUM0
;
1637 SCM result
= SCM_UNDEFINED
;
1643 if (scm_is_eq (m
, SCM_INUM0
))
1645 report_overflow
= 1;
1649 if (!coerce_to_big (n
, n_tmp
))
1651 value_of_wrong_type
= n
;
1652 position_of_wrong_type
= 1;
1656 if (!coerce_to_big (k
, k_tmp
))
1658 value_of_wrong_type
= k
;
1659 position_of_wrong_type
= 2;
1663 if (!coerce_to_big (m
, m_tmp
))
1665 value_of_wrong_type
= m
;
1666 position_of_wrong_type
= 3;
1670 /* if the exponent K is negative, and we simply call mpz_powm, we
1671 will get a divide-by-zero exception when an inverse 1/n mod m
1672 doesn't exist (or is not unique). Since exceptions are hard to
1673 handle, we'll attempt the inversion "by hand" -- that way, we get
1674 a simple failure code, which is easy to handle. */
1676 if (-1 == mpz_sgn (k_tmp
))
1678 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1680 report_overflow
= 1;
1683 mpz_neg (k_tmp
, k_tmp
);
1686 result
= scm_i_mkbig ();
1687 mpz_powm (SCM_I_BIG_MPZ (result
),
1692 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1693 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1700 if (report_overflow
)
1701 scm_num_overflow (FUNC_NAME
);
1703 if (position_of_wrong_type
)
1704 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1705 value_of_wrong_type
);
1707 return scm_i_normbig (result
);
1711 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1713 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1714 "exact integer, @var{n} can be any number.\n"
1716 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1717 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1718 "includes @math{0^0} is 1.\n"
1721 "(integer-expt 2 5) @result{} 32\n"
1722 "(integer-expt -3 3) @result{} -27\n"
1723 "(integer-expt 5 -3) @result{} 1/125\n"
1724 "(integer-expt 0 0) @result{} 1\n"
1726 #define FUNC_NAME s_scm_integer_expt
1729 SCM z_i2
= SCM_BOOL_F
;
1731 SCM acc
= SCM_I_MAKINUM (1L);
1733 /* 0^0 == 1 according to R5RS */
1734 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1735 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1736 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1737 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1739 if (SCM_I_INUMP (k
))
1740 i2
= SCM_I_INUM (k
);
1741 else if (SCM_BIGP (k
))
1743 z_i2
= scm_i_clonebig (k
, 1);
1744 scm_remember_upto_here_1 (k
);
1748 SCM_WRONG_TYPE_ARG (2, k
);
1752 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1754 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1755 n
= scm_divide (n
, SCM_UNDEFINED
);
1759 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1763 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1765 return scm_product (acc
, n
);
1767 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1768 acc
= scm_product (acc
, n
);
1769 n
= scm_product (n
, n
);
1770 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1778 n
= scm_divide (n
, SCM_UNDEFINED
);
1785 return scm_product (acc
, n
);
1787 acc
= scm_product (acc
, n
);
1788 n
= scm_product (n
, n
);
1795 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1797 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1798 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1800 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1801 "@var{cnt} is negative it's a division, rounded towards negative\n"
1802 "infinity. (Note that this is not the same rounding as\n"
1803 "@code{quotient} does.)\n"
1805 "With @var{n} viewed as an infinite precision twos complement,\n"
1806 "@code{ash} means a left shift introducing zero bits, or a right\n"
1807 "shift dropping bits.\n"
1810 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1811 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1813 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1814 "(ash -23 -2) @result{} -6\n"
1816 #define FUNC_NAME s_scm_ash
1819 bits_to_shift
= scm_to_long (cnt
);
1821 if (SCM_I_INUMP (n
))
1823 long nn
= SCM_I_INUM (n
);
1825 if (bits_to_shift
> 0)
1827 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1828 overflow a non-zero fixnum. For smaller shifts we check the
1829 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1830 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1831 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1837 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1839 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1842 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1846 SCM result
= scm_i_long2big (nn
);
1847 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1854 bits_to_shift
= -bits_to_shift
;
1855 if (bits_to_shift
>= SCM_LONG_BIT
)
1856 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1858 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1862 else if (SCM_BIGP (n
))
1866 if (bits_to_shift
== 0)
1869 result
= scm_i_mkbig ();
1870 if (bits_to_shift
>= 0)
1872 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1878 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1879 we have to allocate a bignum even if the result is going to be a
1881 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1883 return scm_i_normbig (result
);
1889 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1895 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1896 (SCM n
, SCM start
, SCM end
),
1897 "Return the integer composed of the @var{start} (inclusive)\n"
1898 "through @var{end} (exclusive) bits of @var{n}. The\n"
1899 "@var{start}th bit becomes the 0-th bit in the result.\n"
1902 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1903 " @result{} \"1010\"\n"
1904 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1905 " @result{} \"10110\"\n"
1907 #define FUNC_NAME s_scm_bit_extract
1909 unsigned long int istart
, iend
, bits
;
1910 istart
= scm_to_ulong (start
);
1911 iend
= scm_to_ulong (end
);
1912 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1914 /* how many bits to keep */
1915 bits
= iend
- istart
;
1917 if (SCM_I_INUMP (n
))
1919 long int in
= SCM_I_INUM (n
);
1921 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1922 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1923 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1925 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1927 /* Since we emulate two's complement encoded numbers, this
1928 * special case requires us to produce a result that has
1929 * more bits than can be stored in a fixnum.
1931 SCM result
= scm_i_long2big (in
);
1932 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1937 /* mask down to requisite bits */
1938 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1939 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1941 else if (SCM_BIGP (n
))
1946 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1950 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1951 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1952 such bits into a ulong. */
1953 result
= scm_i_mkbig ();
1954 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1955 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1956 result
= scm_i_normbig (result
);
1958 scm_remember_upto_here_1 (n
);
1962 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1967 static const char scm_logtab
[] = {
1968 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1971 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1973 "Return the number of bits in integer @var{n}. If integer is\n"
1974 "positive, the 1-bits in its binary representation are counted.\n"
1975 "If negative, the 0-bits in its two's-complement binary\n"
1976 "representation are counted. If 0, 0 is returned.\n"
1979 "(logcount #b10101010)\n"
1986 #define FUNC_NAME s_scm_logcount
1988 if (SCM_I_INUMP (n
))
1990 unsigned long int c
= 0;
1991 long int nn
= SCM_I_INUM (n
);
1996 c
+= scm_logtab
[15 & nn
];
1999 return SCM_I_MAKINUM (c
);
2001 else if (SCM_BIGP (n
))
2003 unsigned long count
;
2004 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2005 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2007 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2008 scm_remember_upto_here_1 (n
);
2009 return SCM_I_MAKINUM (count
);
2012 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2017 static const char scm_ilentab
[] = {
2018 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2022 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2024 "Return the number of bits necessary to represent @var{n}.\n"
2027 "(integer-length #b10101010)\n"
2029 "(integer-length 0)\n"
2031 "(integer-length #b1111)\n"
2034 #define FUNC_NAME s_scm_integer_length
2036 if (SCM_I_INUMP (n
))
2038 unsigned long int c
= 0;
2040 long int nn
= SCM_I_INUM (n
);
2046 l
= scm_ilentab
[15 & nn
];
2049 return SCM_I_MAKINUM (c
- 4 + l
);
2051 else if (SCM_BIGP (n
))
2053 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2054 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2055 1 too big, so check for that and adjust. */
2056 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2057 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2058 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2059 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2061 scm_remember_upto_here_1 (n
);
2062 return SCM_I_MAKINUM (size
);
2065 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2069 /*** NUMBERS -> STRINGS ***/
2070 #define SCM_MAX_DBL_PREC 60
2071 #define SCM_MAX_DBL_RADIX 36
2073 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2074 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2075 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2078 void init_dblprec(int *prec
, int radix
) {
2079 /* determine floating point precision by adding successively
2080 smaller increments to 1.0 until it is considered == 1.0 */
2081 double f
= ((double)1.0)/radix
;
2082 double fsum
= 1.0 + f
;
2087 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2099 void init_fx_radix(double *fx_list
, int radix
)
2101 /* initialize a per-radix list of tolerances. When added
2102 to a number < 1.0, we can determine if we should raund
2103 up and quit converting a number to a string. */
2107 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2108 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2111 /* use this array as a way to generate a single digit */
2112 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2115 idbl2str (double f
, char *a
, int radix
)
2117 int efmt
, dpt
, d
, i
, wp
;
2119 #ifdef DBL_MIN_10_EXP
2122 #endif /* DBL_MIN_10_EXP */
2127 radix
> SCM_MAX_DBL_RADIX
)
2129 /* revert to existing behavior */
2133 wp
= scm_dblprec
[radix
-2];
2134 fx
= fx_per_radix
[radix
-2];
2138 #ifdef HAVE_COPYSIGN
2139 double sgn
= copysign (1.0, f
);
2144 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2150 strcpy (a
, "-inf.0");
2152 strcpy (a
, "+inf.0");
2155 else if (xisnan (f
))
2157 strcpy (a
, "+nan.0");
2167 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2168 make-uniform-vector, from causing infinite loops. */
2169 /* just do the checking...if it passes, we do the conversion for our
2170 radix again below */
2177 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2185 while (f_cpy
> 10.0)
2188 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2209 if (f
+ fx
[wp
] >= radix
)
2216 /* adding 9999 makes this equivalent to abs(x) % 3 */
2217 dpt
= (exp
+ 9999) % 3;
2221 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2243 a
[ch
++] = number_chars
[d
];
2246 if (f
+ fx
[wp
] >= 1.0)
2248 a
[ch
- 1] = number_chars
[d
+1];
2260 if ((dpt
> 4) && (exp
> 6))
2262 d
= (a
[0] == '-' ? 2 : 1);
2263 for (i
= ch
++; i
> d
; i
--)
2276 if (a
[ch
- 1] == '.')
2277 a
[ch
++] = '0'; /* trailing zero */
2286 for (i
= radix
; i
<= exp
; i
*= radix
);
2287 for (i
/= radix
; i
; i
/= radix
)
2289 a
[ch
++] = number_chars
[exp
/ i
];
2298 icmplx2str (double real
, double imag
, char *str
, int radix
)
2302 i
= idbl2str (real
, str
, radix
);
2305 /* Don't output a '+' for negative numbers or for Inf and
2306 NaN. They will provide their own sign. */
2307 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2309 i
+= idbl2str (imag
, &str
[i
], radix
);
2316 iflo2str (SCM flt
, char *str
, int radix
)
2319 if (SCM_REALP (flt
))
2320 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2322 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2327 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2328 characters in the result.
2330 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2332 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2337 return scm_iuint2str (-num
, rad
, p
) + 1;
2340 return scm_iuint2str (num
, rad
, p
);
2343 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2344 characters in the result.
2346 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2348 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2352 scm_t_uintmax n
= num
;
2354 for (n
/= rad
; n
> 0; n
/= rad
)
2364 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2369 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2371 "Return a string holding the external representation of the\n"
2372 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2373 "inexact, a radix of 10 will be used.")
2374 #define FUNC_NAME s_scm_number_to_string
2378 if (SCM_UNBNDP (radix
))
2381 base
= scm_to_signed_integer (radix
, 2, 36);
2383 if (SCM_I_INUMP (n
))
2385 char num_buf
[SCM_INTBUFLEN
];
2386 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2387 return scm_from_locale_stringn (num_buf
, length
);
2389 else if (SCM_BIGP (n
))
2391 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2392 scm_remember_upto_here_1 (n
);
2393 return scm_take_locale_string (str
);
2395 else if (SCM_FRACTIONP (n
))
2397 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2398 scm_from_locale_string ("/"),
2399 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2401 else if (SCM_INEXACTP (n
))
2403 char num_buf
[FLOBUFLEN
];
2404 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2407 SCM_WRONG_TYPE_ARG (1, n
);
2412 /* These print routines used to be stubbed here so that scm_repl.c
2413 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2416 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2418 char num_buf
[FLOBUFLEN
];
2419 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2424 scm_i_print_double (double val
, SCM port
)
2426 char num_buf
[FLOBUFLEN
];
2427 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2431 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2434 char num_buf
[FLOBUFLEN
];
2435 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2440 scm_i_print_complex (double real
, double imag
, SCM port
)
2442 char num_buf
[FLOBUFLEN
];
2443 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2447 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2450 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2451 scm_lfwrite_str (str
, port
);
2452 scm_remember_upto_here_1 (str
);
2457 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2459 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2460 scm_remember_upto_here_1 (exp
);
2461 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2465 /*** END nums->strs ***/
2468 /*** STRINGS -> NUMBERS ***/
2470 /* The following functions implement the conversion from strings to numbers.
2471 * The implementation somehow follows the grammar for numbers as it is given
2472 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2473 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2474 * points should be noted about the implementation:
2475 * * Each function keeps a local index variable 'idx' that points at the
2476 * current position within the parsed string. The global index is only
2477 * updated if the function could parse the corresponding syntactic unit
2479 * * Similarly, the functions keep track of indicators of inexactness ('#',
2480 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2481 * global exactness information is only updated after each part has been
2482 * successfully parsed.
2483 * * Sequences of digits are parsed into temporary variables holding fixnums.
2484 * Only if these fixnums would overflow, the result variables are updated
2485 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2486 * the temporary variables holding the fixnums are cleared, and the process
2487 * starts over again. If for example fixnums were able to store five decimal
2488 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2489 * and the result was computed as 12345 * 100000 + 67890. In other words,
2490 * only every five digits two bignum operations were performed.
2493 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2495 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2497 /* In non ASCII-style encodings the following macro might not work. */
2498 #define XDIGIT2UINT(d) \
2499 (uc_is_property_decimal_digit ((int) (unsigned char) d) \
2501 : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
2504 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2505 unsigned int radix
, enum t_exactness
*p_exactness
)
2507 unsigned int idx
= *p_idx
;
2508 unsigned int hash_seen
= 0;
2509 scm_t_bits shift
= 1;
2511 unsigned int digit_value
;
2514 size_t len
= scm_i_string_length (mem
);
2519 c
= scm_i_string_ref (mem
, idx
);
2520 if (!uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2522 digit_value
= XDIGIT2UINT (c
);
2523 if (digit_value
>= radix
)
2527 result
= SCM_I_MAKINUM (digit_value
);
2530 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2531 if (uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2535 digit_value
= XDIGIT2UINT (c
);
2536 if (digit_value
>= radix
)
2548 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2550 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2552 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2559 shift
= shift
* radix
;
2560 add
= add
* radix
+ digit_value
;
2565 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2567 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2571 *p_exactness
= INEXACT
;
2577 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2578 * covers the parts of the rules that start at a potential point. The value
2579 * of the digits up to the point have been parsed by the caller and are given
2580 * in variable result. The content of *p_exactness indicates, whether a hash
2581 * has already been seen in the digits before the point.
2584 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2587 mem2decimal_from_point (SCM result
, SCM mem
,
2588 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2590 unsigned int idx
= *p_idx
;
2591 enum t_exactness x
= *p_exactness
;
2592 size_t len
= scm_i_string_length (mem
);
2597 if (scm_i_string_ref (mem
, idx
) == '.')
2599 scm_t_bits shift
= 1;
2601 unsigned int digit_value
;
2602 SCM big_shift
= SCM_I_MAKINUM (1);
2607 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2608 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2613 digit_value
= DIGIT2UINT (c
);
2624 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2626 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2627 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2629 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2637 add
= add
* 10 + digit_value
;
2643 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2644 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2645 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2648 result
= scm_divide (result
, big_shift
);
2650 /* We've seen a decimal point, thus the value is implicitly inexact. */
2662 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2664 switch (scm_i_string_ref (mem
, idx
))
2676 c
= scm_i_string_ref (mem
, idx
);
2684 c
= scm_i_string_ref (mem
, idx
);
2693 c
= scm_i_string_ref (mem
, idx
);
2698 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2702 exponent
= DIGIT2UINT (c
);
2705 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2706 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2709 if (exponent
<= SCM_MAXEXP
)
2710 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2716 if (exponent
> SCM_MAXEXP
)
2718 size_t exp_len
= idx
- start
;
2719 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2720 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2721 scm_out_of_range ("string->number", exp_num
);
2724 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2726 result
= scm_product (result
, e
);
2728 result
= scm_divide2real (result
, e
);
2730 /* We've seen an exponent, thus the value is implicitly inexact. */
2748 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2751 mem2ureal (SCM mem
, unsigned int *p_idx
,
2752 unsigned int radix
, enum t_exactness
*p_exactness
)
2754 unsigned int idx
= *p_idx
;
2756 size_t len
= scm_i_string_length (mem
);
2758 /* Start off believing that the number will be exact. This changes
2759 to INEXACT if we see a decimal point or a hash. */
2760 enum t_exactness x
= EXACT
;
2765 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2771 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2773 /* Cobble up the fractional part. We might want to set the
2774 NaN's mantissa from it. */
2776 mem2uinteger (mem
, &idx
, 10, &x
);
2781 if (scm_i_string_ref (mem
, idx
) == '.')
2785 else if (idx
+ 1 == len
)
2787 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2790 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2797 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2798 if (scm_is_false (uinteger
))
2803 else if (scm_i_string_ref (mem
, idx
) == '/')
2811 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2812 if (scm_is_false (divisor
))
2815 /* both are int/big here, I assume */
2816 result
= scm_i_make_ratio (uinteger
, divisor
);
2818 else if (radix
== 10)
2820 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2821 if (scm_is_false (result
))
2830 /* Update *p_exactness if the number just read was inexact. This is
2831 important for complex numbers, so that a complex number is
2832 treated as inexact overall if either its real or imaginary part
2838 /* When returning an inexact zero, make sure it is represented as a
2839 floating point value so that we can change its sign.
2841 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2842 result
= scm_from_double (0.0);
2848 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2851 mem2complex (SCM mem
, unsigned int idx
,
2852 unsigned int radix
, enum t_exactness
*p_exactness
)
2857 size_t len
= scm_i_string_length (mem
);
2862 c
= scm_i_string_ref (mem
, idx
);
2877 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2878 if (scm_is_false (ureal
))
2880 /* input must be either +i or -i */
2885 if (scm_i_string_ref (mem
, idx
) == 'i'
2886 || scm_i_string_ref (mem
, idx
) == 'I')
2892 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2899 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2900 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2905 c
= scm_i_string_ref (mem
, idx
);
2909 /* either +<ureal>i or -<ureal>i */
2916 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2919 /* polar input: <real>@<real>. */
2930 c
= scm_i_string_ref (mem
, idx
);
2948 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2949 if (scm_is_false (angle
))
2954 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2955 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2957 result
= scm_make_polar (ureal
, angle
);
2962 /* expecting input matching <real>[+-]<ureal>?i */
2969 int sign
= (c
== '+') ? 1 : -1;
2970 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2972 if (scm_is_false (imag
))
2973 imag
= SCM_I_MAKINUM (sign
);
2974 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2975 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2979 if (scm_i_string_ref (mem
, idx
) != 'i'
2980 && scm_i_string_ref (mem
, idx
) != 'I')
2987 return scm_make_rectangular (ureal
, imag
);
2996 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2998 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3001 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3003 unsigned int idx
= 0;
3004 unsigned int radix
= NO_RADIX
;
3005 enum t_exactness forced_x
= NO_EXACTNESS
;
3006 enum t_exactness implicit_x
= EXACT
;
3008 size_t len
= scm_i_string_length (mem
);
3010 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3011 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3013 switch (scm_i_string_ref (mem
, idx
+ 1))
3016 if (radix
!= NO_RADIX
)
3021 if (radix
!= NO_RADIX
)
3026 if (forced_x
!= NO_EXACTNESS
)
3031 if (forced_x
!= NO_EXACTNESS
)
3036 if (radix
!= NO_RADIX
)
3041 if (radix
!= NO_RADIX
)
3051 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3052 if (radix
== NO_RADIX
)
3053 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3055 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3057 if (scm_is_false (result
))
3063 if (SCM_INEXACTP (result
))
3064 return scm_inexact_to_exact (result
);
3068 if (SCM_INEXACTP (result
))
3071 return scm_exact_to_inexact (result
);
3074 if (implicit_x
== INEXACT
)
3076 if (SCM_INEXACTP (result
))
3079 return scm_exact_to_inexact (result
);
3087 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3088 unsigned int default_radix
)
3090 SCM str
= scm_from_locale_stringn (mem
, len
);
3092 return scm_i_string_to_number (str
, default_radix
);
3096 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3097 (SCM string
, SCM radix
),
3098 "Return a number of the maximally precise representation\n"
3099 "expressed by the given @var{string}. @var{radix} must be an\n"
3100 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3101 "is a default radix that may be overridden by an explicit radix\n"
3102 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3103 "supplied, then the default radix is 10. If string is not a\n"
3104 "syntactically valid notation for a number, then\n"
3105 "@code{string->number} returns @code{#f}.")
3106 #define FUNC_NAME s_scm_string_to_number
3110 SCM_VALIDATE_STRING (1, string
);
3112 if (SCM_UNBNDP (radix
))
3115 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3117 answer
= scm_i_string_to_number (string
, base
);
3118 scm_remember_upto_here_1 (string
);
3124 /*** END strs->nums ***/
3128 scm_bigequal (SCM x
, SCM y
)
3130 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3131 scm_remember_upto_here_2 (x
, y
);
3132 return scm_from_bool (0 == result
);
3136 scm_real_equalp (SCM x
, SCM y
)
3138 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3142 scm_complex_equalp (SCM x
, SCM y
)
3144 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3145 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3149 scm_i_fraction_equalp (SCM x
, SCM y
)
3151 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3152 SCM_FRACTION_NUMERATOR (y
)))
3153 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3154 SCM_FRACTION_DENOMINATOR (y
))))
3161 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3163 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3165 #define FUNC_NAME s_scm_number_p
3167 return scm_from_bool (SCM_NUMBERP (x
));
3171 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3173 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3174 "otherwise. Note that the sets of real, rational and integer\n"
3175 "values form subsets of the set of complex numbers, i. e. the\n"
3176 "predicate will also be fulfilled if @var{x} is a real,\n"
3177 "rational or integer number.")
3178 #define FUNC_NAME s_scm_complex_p
3180 /* all numbers are complex. */
3181 return scm_number_p (x
);
3185 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3187 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3188 "otherwise. Note that the set of integer values forms a subset of\n"
3189 "the set of real numbers, i. e. the predicate will also be\n"
3190 "fulfilled if @var{x} is an integer number.")
3191 #define FUNC_NAME s_scm_real_p
3193 /* we can't represent irrational numbers. */
3194 return scm_rational_p (x
);
3198 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3200 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3201 "otherwise. Note that the set of integer values forms a subset of\n"
3202 "the set of rational numbers, i. e. the predicate will also be\n"
3203 "fulfilled if @var{x} is an integer number.")
3204 #define FUNC_NAME s_scm_rational_p
3206 if (SCM_I_INUMP (x
))
3208 else if (SCM_IMP (x
))
3210 else if (SCM_BIGP (x
))
3212 else if (SCM_FRACTIONP (x
))
3214 else if (SCM_REALP (x
))
3215 /* due to their limited precision, all floating point numbers are
3216 rational as well. */
3223 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3225 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3227 #define FUNC_NAME s_scm_integer_p
3230 if (SCM_I_INUMP (x
))
3236 if (!SCM_INEXACTP (x
))
3238 if (SCM_COMPLEXP (x
))
3240 r
= SCM_REAL_VALUE (x
);
3241 /* +/-inf passes r==floor(r), making those #t */
3249 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3251 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3253 #define FUNC_NAME s_scm_inexact_p
3255 if (SCM_INEXACTP (x
))
3257 if (SCM_NUMBERP (x
))
3259 SCM_WRONG_TYPE_ARG (1, x
);
3264 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3265 /* "Return @code{#t} if all parameters are numerically equal." */
3267 scm_num_eq_p (SCM x
, SCM y
)
3270 if (SCM_I_INUMP (x
))
3272 long xx
= SCM_I_INUM (x
);
3273 if (SCM_I_INUMP (y
))
3275 long yy
= SCM_I_INUM (y
);
3276 return scm_from_bool (xx
== yy
);
3278 else if (SCM_BIGP (y
))
3280 else if (SCM_REALP (y
))
3282 /* On a 32-bit system an inum fits a double, we can cast the inum
3283 to a double and compare.
3285 But on a 64-bit system an inum is bigger than a double and
3286 casting it to a double (call that dxx) will round. dxx is at
3287 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3288 an integer and fits a long. So we cast yy to a long and
3289 compare with plain xx.
3291 An alternative (for any size system actually) would be to check
3292 yy is an integer (with floor) and is in range of an inum
3293 (compare against appropriate powers of 2) then test
3294 xx==(long)yy. It's just a matter of which casts/comparisons
3295 might be fastest or easiest for the cpu. */
3297 double yy
= SCM_REAL_VALUE (y
);
3298 return scm_from_bool ((double) xx
== yy
3299 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3300 || xx
== (long) yy
));
3302 else if (SCM_COMPLEXP (y
))
3303 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3304 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3305 else if (SCM_FRACTIONP (y
))
3308 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3310 else if (SCM_BIGP (x
))
3312 if (SCM_I_INUMP (y
))
3314 else if (SCM_BIGP (y
))
3316 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3317 scm_remember_upto_here_2 (x
, y
);
3318 return scm_from_bool (0 == cmp
);
3320 else if (SCM_REALP (y
))
3323 if (xisnan (SCM_REAL_VALUE (y
)))
3325 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3326 scm_remember_upto_here_1 (x
);
3327 return scm_from_bool (0 == cmp
);
3329 else if (SCM_COMPLEXP (y
))
3332 if (0.0 != SCM_COMPLEX_IMAG (y
))
3334 if (xisnan (SCM_COMPLEX_REAL (y
)))
3336 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3337 scm_remember_upto_here_1 (x
);
3338 return scm_from_bool (0 == cmp
);
3340 else if (SCM_FRACTIONP (y
))
3343 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3345 else if (SCM_REALP (x
))
3347 double xx
= SCM_REAL_VALUE (x
);
3348 if (SCM_I_INUMP (y
))
3350 /* see comments with inum/real above */
3351 long yy
= SCM_I_INUM (y
);
3352 return scm_from_bool (xx
== (double) yy
3353 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3354 || (long) xx
== yy
));
3356 else if (SCM_BIGP (y
))
3359 if (xisnan (SCM_REAL_VALUE (x
)))
3361 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3362 scm_remember_upto_here_1 (y
);
3363 return scm_from_bool (0 == cmp
);
3365 else if (SCM_REALP (y
))
3366 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3367 else if (SCM_COMPLEXP (y
))
3368 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3369 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3370 else if (SCM_FRACTIONP (y
))
3372 double xx
= SCM_REAL_VALUE (x
);
3376 return scm_from_bool (xx
< 0.0);
3377 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3381 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3383 else if (SCM_COMPLEXP (x
))
3385 if (SCM_I_INUMP (y
))
3386 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3387 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3388 else if (SCM_BIGP (y
))
3391 if (0.0 != SCM_COMPLEX_IMAG (x
))
3393 if (xisnan (SCM_COMPLEX_REAL (x
)))
3395 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3396 scm_remember_upto_here_1 (y
);
3397 return scm_from_bool (0 == cmp
);
3399 else if (SCM_REALP (y
))
3400 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3401 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3402 else if (SCM_COMPLEXP (y
))
3403 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3404 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3405 else if (SCM_FRACTIONP (y
))
3408 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3410 xx
= SCM_COMPLEX_REAL (x
);
3414 return scm_from_bool (xx
< 0.0);
3415 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3419 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3421 else if (SCM_FRACTIONP (x
))
3423 if (SCM_I_INUMP (y
))
3425 else if (SCM_BIGP (y
))
3427 else if (SCM_REALP (y
))
3429 double yy
= SCM_REAL_VALUE (y
);
3433 return scm_from_bool (0.0 < yy
);
3434 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3437 else if (SCM_COMPLEXP (y
))
3440 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3442 yy
= SCM_COMPLEX_REAL (y
);
3446 return scm_from_bool (0.0 < yy
);
3447 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3450 else if (SCM_FRACTIONP (y
))
3451 return scm_i_fraction_equalp (x
, y
);
3453 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3456 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3460 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3461 done are good for inums, but for bignums an answer can almost always be
3462 had by just examining a few high bits of the operands, as done by GMP in
3463 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3464 of the float exponent to take into account. */
3466 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3467 /* "Return @code{#t} if the list of parameters is monotonically\n"
3471 scm_less_p (SCM x
, SCM y
)
3474 if (SCM_I_INUMP (x
))
3476 long xx
= SCM_I_INUM (x
);
3477 if (SCM_I_INUMP (y
))
3479 long yy
= SCM_I_INUM (y
);
3480 return scm_from_bool (xx
< yy
);
3482 else if (SCM_BIGP (y
))
3484 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3485 scm_remember_upto_here_1 (y
);
3486 return scm_from_bool (sgn
> 0);
3488 else if (SCM_REALP (y
))
3489 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3490 else if (SCM_FRACTIONP (y
))
3492 /* "x < a/b" becomes "x*b < a" */
3494 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3495 y
= SCM_FRACTION_NUMERATOR (y
);
3499 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3501 else if (SCM_BIGP (x
))
3503 if (SCM_I_INUMP (y
))
3505 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3506 scm_remember_upto_here_1 (x
);
3507 return scm_from_bool (sgn
< 0);
3509 else if (SCM_BIGP (y
))
3511 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3512 scm_remember_upto_here_2 (x
, y
);
3513 return scm_from_bool (cmp
< 0);
3515 else if (SCM_REALP (y
))
3518 if (xisnan (SCM_REAL_VALUE (y
)))
3520 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3521 scm_remember_upto_here_1 (x
);
3522 return scm_from_bool (cmp
< 0);
3524 else if (SCM_FRACTIONP (y
))
3527 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3529 else if (SCM_REALP (x
))
3531 if (SCM_I_INUMP (y
))
3532 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3533 else if (SCM_BIGP (y
))
3536 if (xisnan (SCM_REAL_VALUE (x
)))
3538 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3539 scm_remember_upto_here_1 (y
);
3540 return scm_from_bool (cmp
> 0);
3542 else if (SCM_REALP (y
))
3543 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3544 else if (SCM_FRACTIONP (y
))
3546 double xx
= SCM_REAL_VALUE (x
);
3550 return scm_from_bool (xx
< 0.0);
3551 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3555 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3557 else if (SCM_FRACTIONP (x
))
3559 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3561 /* "a/b < y" becomes "a < y*b" */
3562 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3563 x
= SCM_FRACTION_NUMERATOR (x
);
3566 else if (SCM_REALP (y
))
3568 double yy
= SCM_REAL_VALUE (y
);
3572 return scm_from_bool (0.0 < yy
);
3573 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3576 else if (SCM_FRACTIONP (y
))
3578 /* "a/b < c/d" becomes "a*d < c*b" */
3579 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3580 SCM_FRACTION_DENOMINATOR (y
));
3581 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3582 SCM_FRACTION_DENOMINATOR (x
));
3588 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3591 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3595 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3596 /* "Return @code{#t} if the list of parameters is monotonically\n"
3599 #define FUNC_NAME s_scm_gr_p
3601 scm_gr_p (SCM x
, SCM y
)
3603 if (!SCM_NUMBERP (x
))
3604 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3605 else if (!SCM_NUMBERP (y
))
3606 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3608 return scm_less_p (y
, x
);
3613 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3614 /* "Return @code{#t} if the list of parameters is monotonically\n"
3617 #define FUNC_NAME s_scm_leq_p
3619 scm_leq_p (SCM x
, SCM y
)
3621 if (!SCM_NUMBERP (x
))
3622 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3623 else if (!SCM_NUMBERP (y
))
3624 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3625 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3628 return scm_not (scm_less_p (y
, x
));
3633 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3634 /* "Return @code{#t} if the list of parameters is monotonically\n"
3637 #define FUNC_NAME s_scm_geq_p
3639 scm_geq_p (SCM x
, SCM y
)
3641 if (!SCM_NUMBERP (x
))
3642 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3643 else if (!SCM_NUMBERP (y
))
3644 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3645 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3648 return scm_not (scm_less_p (x
, y
));
3653 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3654 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3660 if (SCM_I_INUMP (z
))
3661 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3662 else if (SCM_BIGP (z
))
3664 else if (SCM_REALP (z
))
3665 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3666 else if (SCM_COMPLEXP (z
))
3667 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3668 && SCM_COMPLEX_IMAG (z
) == 0.0);
3669 else if (SCM_FRACTIONP (z
))
3672 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3676 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3677 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3681 scm_positive_p (SCM x
)
3683 if (SCM_I_INUMP (x
))
3684 return scm_from_bool (SCM_I_INUM (x
) > 0);
3685 else if (SCM_BIGP (x
))
3687 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3688 scm_remember_upto_here_1 (x
);
3689 return scm_from_bool (sgn
> 0);
3691 else if (SCM_REALP (x
))
3692 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3693 else if (SCM_FRACTIONP (x
))
3694 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3696 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3700 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3701 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3705 scm_negative_p (SCM x
)
3707 if (SCM_I_INUMP (x
))
3708 return scm_from_bool (SCM_I_INUM (x
) < 0);
3709 else if (SCM_BIGP (x
))
3711 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3712 scm_remember_upto_here_1 (x
);
3713 return scm_from_bool (sgn
< 0);
3715 else if (SCM_REALP (x
))
3716 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3717 else if (SCM_FRACTIONP (x
))
3718 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3720 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3724 /* scm_min and scm_max return an inexact when either argument is inexact, as
3725 required by r5rs. On that basis, for exact/inexact combinations the
3726 exact is converted to inexact to compare and possibly return. This is
3727 unlike scm_less_p above which takes some trouble to preserve all bits in
3728 its test, such trouble is not required for min and max. */
3730 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3731 /* "Return the maximum of all parameter values."
3734 scm_max (SCM x
, SCM y
)
3739 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3740 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3743 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3746 if (SCM_I_INUMP (x
))
3748 long xx
= SCM_I_INUM (x
);
3749 if (SCM_I_INUMP (y
))
3751 long yy
= SCM_I_INUM (y
);
3752 return (xx
< yy
) ? y
: x
;
3754 else if (SCM_BIGP (y
))
3756 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3757 scm_remember_upto_here_1 (y
);
3758 return (sgn
< 0) ? x
: y
;
3760 else if (SCM_REALP (y
))
3763 /* if y==NaN then ">" is false and we return NaN */
3764 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3766 else if (SCM_FRACTIONP (y
))
3769 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3772 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3774 else if (SCM_BIGP (x
))
3776 if (SCM_I_INUMP (y
))
3778 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3779 scm_remember_upto_here_1 (x
);
3780 return (sgn
< 0) ? y
: x
;
3782 else if (SCM_BIGP (y
))
3784 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3785 scm_remember_upto_here_2 (x
, y
);
3786 return (cmp
> 0) ? x
: y
;
3788 else if (SCM_REALP (y
))
3790 /* if y==NaN then xx>yy is false, so we return the NaN y */
3793 xx
= scm_i_big2dbl (x
);
3794 yy
= SCM_REAL_VALUE (y
);
3795 return (xx
> yy
? scm_from_double (xx
) : y
);
3797 else if (SCM_FRACTIONP (y
))
3802 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3804 else if (SCM_REALP (x
))
3806 if (SCM_I_INUMP (y
))
3808 double z
= SCM_I_INUM (y
);
3809 /* if x==NaN then "<" is false and we return NaN */
3810 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3812 else if (SCM_BIGP (y
))
3817 else if (SCM_REALP (y
))
3819 /* if x==NaN then our explicit check means we return NaN
3820 if y==NaN then ">" is false and we return NaN
3821 calling isnan is unavoidable, since it's the only way to know
3822 which of x or y causes any compares to be false */
3823 double xx
= SCM_REAL_VALUE (x
);
3824 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3826 else if (SCM_FRACTIONP (y
))
3828 double yy
= scm_i_fraction2double (y
);
3829 double xx
= SCM_REAL_VALUE (x
);
3830 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3833 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3835 else if (SCM_FRACTIONP (x
))
3837 if (SCM_I_INUMP (y
))
3841 else if (SCM_BIGP (y
))
3845 else if (SCM_REALP (y
))
3847 double xx
= scm_i_fraction2double (x
);
3848 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3850 else if (SCM_FRACTIONP (y
))
3855 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3858 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3862 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3863 /* "Return the minium of all parameter values."
3866 scm_min (SCM x
, SCM y
)
3871 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3872 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3875 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3878 if (SCM_I_INUMP (x
))
3880 long xx
= SCM_I_INUM (x
);
3881 if (SCM_I_INUMP (y
))
3883 long yy
= SCM_I_INUM (y
);
3884 return (xx
< yy
) ? x
: y
;
3886 else if (SCM_BIGP (y
))
3888 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3889 scm_remember_upto_here_1 (y
);
3890 return (sgn
< 0) ? y
: x
;
3892 else if (SCM_REALP (y
))
3895 /* if y==NaN then "<" is false and we return NaN */
3896 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3898 else if (SCM_FRACTIONP (y
))
3901 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3904 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3906 else if (SCM_BIGP (x
))
3908 if (SCM_I_INUMP (y
))
3910 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3911 scm_remember_upto_here_1 (x
);
3912 return (sgn
< 0) ? x
: y
;
3914 else if (SCM_BIGP (y
))
3916 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3917 scm_remember_upto_here_2 (x
, y
);
3918 return (cmp
> 0) ? y
: x
;
3920 else if (SCM_REALP (y
))
3922 /* if y==NaN then xx<yy is false, so we return the NaN y */
3925 xx
= scm_i_big2dbl (x
);
3926 yy
= SCM_REAL_VALUE (y
);
3927 return (xx
< yy
? scm_from_double (xx
) : y
);
3929 else if (SCM_FRACTIONP (y
))
3934 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3936 else if (SCM_REALP (x
))
3938 if (SCM_I_INUMP (y
))
3940 double z
= SCM_I_INUM (y
);
3941 /* if x==NaN then "<" is false and we return NaN */
3942 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3944 else if (SCM_BIGP (y
))
3949 else if (SCM_REALP (y
))
3951 /* if x==NaN then our explicit check means we return NaN
3952 if y==NaN then "<" is false and we return NaN
3953 calling isnan is unavoidable, since it's the only way to know
3954 which of x or y causes any compares to be false */
3955 double xx
= SCM_REAL_VALUE (x
);
3956 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3958 else if (SCM_FRACTIONP (y
))
3960 double yy
= scm_i_fraction2double (y
);
3961 double xx
= SCM_REAL_VALUE (x
);
3962 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3965 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3967 else if (SCM_FRACTIONP (x
))
3969 if (SCM_I_INUMP (y
))
3973 else if (SCM_BIGP (y
))
3977 else if (SCM_REALP (y
))
3979 double xx
= scm_i_fraction2double (x
);
3980 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3982 else if (SCM_FRACTIONP (y
))
3987 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3990 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3994 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3995 /* "Return the sum of all parameter values. Return 0 if called without\n"
3999 scm_sum (SCM x
, SCM y
)
4001 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4003 if (SCM_NUMBERP (x
)) return x
;
4004 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4005 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4008 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4010 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4012 long xx
= SCM_I_INUM (x
);
4013 long yy
= SCM_I_INUM (y
);
4014 long int z
= xx
+ yy
;
4015 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4017 else if (SCM_BIGP (y
))
4022 else if (SCM_REALP (y
))
4024 long int xx
= SCM_I_INUM (x
);
4025 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4027 else if (SCM_COMPLEXP (y
))
4029 long int xx
= SCM_I_INUM (x
);
4030 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4031 SCM_COMPLEX_IMAG (y
));
4033 else if (SCM_FRACTIONP (y
))
4034 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4035 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4036 SCM_FRACTION_DENOMINATOR (y
));
4038 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4039 } else if (SCM_BIGP (x
))
4041 if (SCM_I_INUMP (y
))
4046 inum
= SCM_I_INUM (y
);
4049 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4052 SCM result
= scm_i_mkbig ();
4053 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4054 scm_remember_upto_here_1 (x
);
4055 /* we know the result will have to be a bignum */
4058 return scm_i_normbig (result
);
4062 SCM result
= scm_i_mkbig ();
4063 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4064 scm_remember_upto_here_1 (x
);
4065 /* we know the result will have to be a bignum */
4068 return scm_i_normbig (result
);
4071 else if (SCM_BIGP (y
))
4073 SCM result
= scm_i_mkbig ();
4074 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4075 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4076 mpz_add (SCM_I_BIG_MPZ (result
),
4079 scm_remember_upto_here_2 (x
, y
);
4080 /* we know the result will have to be a bignum */
4083 return scm_i_normbig (result
);
4085 else if (SCM_REALP (y
))
4087 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4088 scm_remember_upto_here_1 (x
);
4089 return scm_from_double (result
);
4091 else if (SCM_COMPLEXP (y
))
4093 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4094 + SCM_COMPLEX_REAL (y
));
4095 scm_remember_upto_here_1 (x
);
4096 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4098 else if (SCM_FRACTIONP (y
))
4099 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4100 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4101 SCM_FRACTION_DENOMINATOR (y
));
4103 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4105 else if (SCM_REALP (x
))
4107 if (SCM_I_INUMP (y
))
4108 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4109 else if (SCM_BIGP (y
))
4111 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4112 scm_remember_upto_here_1 (y
);
4113 return scm_from_double (result
);
4115 else if (SCM_REALP (y
))
4116 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4117 else if (SCM_COMPLEXP (y
))
4118 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4119 SCM_COMPLEX_IMAG (y
));
4120 else if (SCM_FRACTIONP (y
))
4121 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4123 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4125 else if (SCM_COMPLEXP (x
))
4127 if (SCM_I_INUMP (y
))
4128 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4129 SCM_COMPLEX_IMAG (x
));
4130 else if (SCM_BIGP (y
))
4132 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4133 + SCM_COMPLEX_REAL (x
));
4134 scm_remember_upto_here_1 (y
);
4135 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4137 else if (SCM_REALP (y
))
4138 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4139 SCM_COMPLEX_IMAG (x
));
4140 else if (SCM_COMPLEXP (y
))
4141 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4142 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4143 else if (SCM_FRACTIONP (y
))
4144 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4145 SCM_COMPLEX_IMAG (x
));
4147 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4149 else if (SCM_FRACTIONP (x
))
4151 if (SCM_I_INUMP (y
))
4152 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4153 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4154 SCM_FRACTION_DENOMINATOR (x
));
4155 else if (SCM_BIGP (y
))
4156 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4157 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4158 SCM_FRACTION_DENOMINATOR (x
));
4159 else if (SCM_REALP (y
))
4160 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4161 else if (SCM_COMPLEXP (y
))
4162 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4163 SCM_COMPLEX_IMAG (y
));
4164 else if (SCM_FRACTIONP (y
))
4165 /* a/b + c/d = (ad + bc) / bd */
4166 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4167 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4168 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4170 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4173 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4177 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4179 "Return @math{@var{x}+1}.")
4180 #define FUNC_NAME s_scm_oneplus
4182 return scm_sum (x
, SCM_I_MAKINUM (1));
4187 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4188 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4189 * the sum of all but the first argument are subtracted from the first
4191 #define FUNC_NAME s_difference
4193 scm_difference (SCM x
, SCM y
)
4195 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4198 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4200 if (SCM_I_INUMP (x
))
4202 long xx
= -SCM_I_INUM (x
);
4203 if (SCM_FIXABLE (xx
))
4204 return SCM_I_MAKINUM (xx
);
4206 return scm_i_long2big (xx
);
4208 else if (SCM_BIGP (x
))
4209 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4210 bignum, but negating that gives a fixnum. */
4211 return scm_i_normbig (scm_i_clonebig (x
, 0));
4212 else if (SCM_REALP (x
))
4213 return scm_from_double (-SCM_REAL_VALUE (x
));
4214 else if (SCM_COMPLEXP (x
))
4215 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4216 -SCM_COMPLEX_IMAG (x
));
4217 else if (SCM_FRACTIONP (x
))
4218 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4219 SCM_FRACTION_DENOMINATOR (x
));
4221 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4224 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4226 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4228 long int xx
= SCM_I_INUM (x
);
4229 long int yy
= SCM_I_INUM (y
);
4230 long int z
= xx
- yy
;
4231 if (SCM_FIXABLE (z
))
4232 return SCM_I_MAKINUM (z
);
4234 return scm_i_long2big (z
);
4236 else if (SCM_BIGP (y
))
4238 /* inum-x - big-y */
4239 long xx
= SCM_I_INUM (x
);
4242 return scm_i_clonebig (y
, 0);
4245 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4246 SCM result
= scm_i_mkbig ();
4249 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4252 /* x - y == -(y + -x) */
4253 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4254 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4256 scm_remember_upto_here_1 (y
);
4258 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4259 /* we know the result will have to be a bignum */
4262 return scm_i_normbig (result
);
4265 else if (SCM_REALP (y
))
4267 long int xx
= SCM_I_INUM (x
);
4268 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4270 else if (SCM_COMPLEXP (y
))
4272 long int xx
= SCM_I_INUM (x
);
4273 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4274 - SCM_COMPLEX_IMAG (y
));
4276 else if (SCM_FRACTIONP (y
))
4277 /* a - b/c = (ac - b) / c */
4278 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4279 SCM_FRACTION_NUMERATOR (y
)),
4280 SCM_FRACTION_DENOMINATOR (y
));
4282 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4284 else if (SCM_BIGP (x
))
4286 if (SCM_I_INUMP (y
))
4288 /* big-x - inum-y */
4289 long yy
= SCM_I_INUM (y
);
4290 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4292 scm_remember_upto_here_1 (x
);
4294 return (SCM_FIXABLE (-yy
) ?
4295 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4298 SCM result
= scm_i_mkbig ();
4301 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4303 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4304 scm_remember_upto_here_1 (x
);
4306 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4307 /* we know the result will have to be a bignum */
4310 return scm_i_normbig (result
);
4313 else if (SCM_BIGP (y
))
4315 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4316 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4317 SCM result
= scm_i_mkbig ();
4318 mpz_sub (SCM_I_BIG_MPZ (result
),
4321 scm_remember_upto_here_2 (x
, y
);
4322 /* we know the result will have to be a bignum */
4323 if ((sgn_x
== 1) && (sgn_y
== -1))
4325 if ((sgn_x
== -1) && (sgn_y
== 1))
4327 return scm_i_normbig (result
);
4329 else if (SCM_REALP (y
))
4331 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4332 scm_remember_upto_here_1 (x
);
4333 return scm_from_double (result
);
4335 else if (SCM_COMPLEXP (y
))
4337 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4338 - SCM_COMPLEX_REAL (y
));
4339 scm_remember_upto_here_1 (x
);
4340 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4342 else if (SCM_FRACTIONP (y
))
4343 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4344 SCM_FRACTION_NUMERATOR (y
)),
4345 SCM_FRACTION_DENOMINATOR (y
));
4346 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4348 else if (SCM_REALP (x
))
4350 if (SCM_I_INUMP (y
))
4351 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4352 else if (SCM_BIGP (y
))
4354 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4355 scm_remember_upto_here_1 (x
);
4356 return scm_from_double (result
);
4358 else if (SCM_REALP (y
))
4359 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4360 else if (SCM_COMPLEXP (y
))
4361 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4362 -SCM_COMPLEX_IMAG (y
));
4363 else if (SCM_FRACTIONP (y
))
4364 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4366 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4368 else if (SCM_COMPLEXP (x
))
4370 if (SCM_I_INUMP (y
))
4371 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4372 SCM_COMPLEX_IMAG (x
));
4373 else if (SCM_BIGP (y
))
4375 double real_part
= (SCM_COMPLEX_REAL (x
)
4376 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4377 scm_remember_upto_here_1 (x
);
4378 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4380 else if (SCM_REALP (y
))
4381 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4382 SCM_COMPLEX_IMAG (x
));
4383 else if (SCM_COMPLEXP (y
))
4384 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4385 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4386 else if (SCM_FRACTIONP (y
))
4387 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4388 SCM_COMPLEX_IMAG (x
));
4390 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4392 else if (SCM_FRACTIONP (x
))
4394 if (SCM_I_INUMP (y
))
4395 /* a/b - c = (a - cb) / b */
4396 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4397 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4398 SCM_FRACTION_DENOMINATOR (x
));
4399 else if (SCM_BIGP (y
))
4400 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4401 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4402 SCM_FRACTION_DENOMINATOR (x
));
4403 else if (SCM_REALP (y
))
4404 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4405 else if (SCM_COMPLEXP (y
))
4406 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4407 -SCM_COMPLEX_IMAG (y
));
4408 else if (SCM_FRACTIONP (y
))
4409 /* a/b - c/d = (ad - bc) / bd */
4410 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4411 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4412 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4414 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4417 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4422 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4424 "Return @math{@var{x}-1}.")
4425 #define FUNC_NAME s_scm_oneminus
4427 return scm_difference (x
, SCM_I_MAKINUM (1));
4432 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4433 /* "Return the product of all arguments. If called without arguments,\n"
4437 scm_product (SCM x
, SCM y
)
4439 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4442 return SCM_I_MAKINUM (1L);
4443 else if (SCM_NUMBERP (x
))
4446 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4449 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4454 xx
= SCM_I_INUM (x
);
4458 case 0: return x
; break;
4459 case 1: return y
; break;
4462 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4464 long yy
= SCM_I_INUM (y
);
4466 SCM k
= SCM_I_MAKINUM (kk
);
4467 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4471 SCM result
= scm_i_long2big (xx
);
4472 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4473 return scm_i_normbig (result
);
4476 else if (SCM_BIGP (y
))
4478 SCM result
= scm_i_mkbig ();
4479 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4480 scm_remember_upto_here_1 (y
);
4483 else if (SCM_REALP (y
))
4484 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4485 else if (SCM_COMPLEXP (y
))
4486 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4487 xx
* SCM_COMPLEX_IMAG (y
));
4488 else if (SCM_FRACTIONP (y
))
4489 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4490 SCM_FRACTION_DENOMINATOR (y
));
4492 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4494 else if (SCM_BIGP (x
))
4496 if (SCM_I_INUMP (y
))
4501 else if (SCM_BIGP (y
))
4503 SCM result
= scm_i_mkbig ();
4504 mpz_mul (SCM_I_BIG_MPZ (result
),
4507 scm_remember_upto_here_2 (x
, y
);
4510 else if (SCM_REALP (y
))
4512 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4513 scm_remember_upto_here_1 (x
);
4514 return scm_from_double (result
);
4516 else if (SCM_COMPLEXP (y
))
4518 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4519 scm_remember_upto_here_1 (x
);
4520 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4521 z
* SCM_COMPLEX_IMAG (y
));
4523 else if (SCM_FRACTIONP (y
))
4524 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4525 SCM_FRACTION_DENOMINATOR (y
));
4527 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4529 else if (SCM_REALP (x
))
4531 if (SCM_I_INUMP (y
))
4533 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4534 if (scm_is_eq (y
, SCM_INUM0
))
4536 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4538 else if (SCM_BIGP (y
))
4540 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4541 scm_remember_upto_here_1 (y
);
4542 return scm_from_double (result
);
4544 else if (SCM_REALP (y
))
4545 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4546 else if (SCM_COMPLEXP (y
))
4547 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4548 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4549 else if (SCM_FRACTIONP (y
))
4550 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4552 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4554 else if (SCM_COMPLEXP (x
))
4556 if (SCM_I_INUMP (y
))
4558 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4559 if (scm_is_eq (y
, SCM_INUM0
))
4561 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4562 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4564 else if (SCM_BIGP (y
))
4566 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4567 scm_remember_upto_here_1 (y
);
4568 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4569 z
* SCM_COMPLEX_IMAG (x
));
4571 else if (SCM_REALP (y
))
4572 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4573 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4574 else if (SCM_COMPLEXP (y
))
4576 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4577 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4578 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4579 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4581 else if (SCM_FRACTIONP (y
))
4583 double yy
= scm_i_fraction2double (y
);
4584 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4585 yy
* SCM_COMPLEX_IMAG (x
));
4588 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4590 else if (SCM_FRACTIONP (x
))
4592 if (SCM_I_INUMP (y
))
4593 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4594 SCM_FRACTION_DENOMINATOR (x
));
4595 else if (SCM_BIGP (y
))
4596 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4597 SCM_FRACTION_DENOMINATOR (x
));
4598 else if (SCM_REALP (y
))
4599 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4600 else if (SCM_COMPLEXP (y
))
4602 double xx
= scm_i_fraction2double (x
);
4603 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4604 xx
* SCM_COMPLEX_IMAG (y
));
4606 else if (SCM_FRACTIONP (y
))
4607 /* a/b * c/d = ac / bd */
4608 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4609 SCM_FRACTION_NUMERATOR (y
)),
4610 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4611 SCM_FRACTION_DENOMINATOR (y
)));
4613 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4616 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4619 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4620 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4621 #define ALLOW_DIVIDE_BY_ZERO
4622 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4625 /* The code below for complex division is adapted from the GNU
4626 libstdc++, which adapted it from f2c's libF77, and is subject to
4629 /****************************************************************
4630 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4632 Permission to use, copy, modify, and distribute this software
4633 and its documentation for any purpose and without fee is hereby
4634 granted, provided that the above copyright notice appear in all
4635 copies and that both that the copyright notice and this
4636 permission notice and warranty disclaimer appear in supporting
4637 documentation, and that the names of AT&T Bell Laboratories or
4638 Bellcore or any of their entities not be used in advertising or
4639 publicity pertaining to distribution of the software without
4640 specific, written prior permission.
4642 AT&T and Bellcore disclaim all warranties with regard to this
4643 software, including all implied warranties of merchantability
4644 and fitness. In no event shall AT&T or Bellcore be liable for
4645 any special, indirect or consequential damages or any damages
4646 whatsoever resulting from loss of use, data or profits, whether
4647 in an action of contract, negligence or other tortious action,
4648 arising out of or in connection with the use or performance of
4650 ****************************************************************/
4652 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4653 /* Divide the first argument by the product of the remaining
4654 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4656 #define FUNC_NAME s_divide
4658 scm_i_divide (SCM x
, SCM y
, int inexact
)
4662 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4665 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4666 else if (SCM_I_INUMP (x
))
4668 long xx
= SCM_I_INUM (x
);
4669 if (xx
== 1 || xx
== -1)
4671 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4673 scm_num_overflow (s_divide
);
4678 return scm_from_double (1.0 / (double) xx
);
4679 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4682 else if (SCM_BIGP (x
))
4685 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4686 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4688 else if (SCM_REALP (x
))
4690 double xx
= SCM_REAL_VALUE (x
);
4691 #ifndef ALLOW_DIVIDE_BY_ZERO
4693 scm_num_overflow (s_divide
);
4696 return scm_from_double (1.0 / xx
);
4698 else if (SCM_COMPLEXP (x
))
4700 double r
= SCM_COMPLEX_REAL (x
);
4701 double i
= SCM_COMPLEX_IMAG (x
);
4702 if (fabs(r
) <= fabs(i
))
4705 double d
= i
* (1.0 + t
* t
);
4706 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4711 double d
= r
* (1.0 + t
* t
);
4712 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4715 else if (SCM_FRACTIONP (x
))
4716 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4717 SCM_FRACTION_NUMERATOR (x
));
4719 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4722 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4724 long xx
= SCM_I_INUM (x
);
4725 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4727 long yy
= SCM_I_INUM (y
);
4730 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4731 scm_num_overflow (s_divide
);
4733 return scm_from_double ((double) xx
/ (double) yy
);
4736 else if (xx
% yy
!= 0)
4739 return scm_from_double ((double) xx
/ (double) yy
);
4740 else return scm_i_make_ratio (x
, y
);
4745 if (SCM_FIXABLE (z
))
4746 return SCM_I_MAKINUM (z
);
4748 return scm_i_long2big (z
);
4751 else if (SCM_BIGP (y
))
4754 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4755 else return scm_i_make_ratio (x
, y
);
4757 else if (SCM_REALP (y
))
4759 double yy
= SCM_REAL_VALUE (y
);
4760 #ifndef ALLOW_DIVIDE_BY_ZERO
4762 scm_num_overflow (s_divide
);
4765 return scm_from_double ((double) xx
/ yy
);
4767 else if (SCM_COMPLEXP (y
))
4770 complex_div
: /* y _must_ be a complex number */
4772 double r
= SCM_COMPLEX_REAL (y
);
4773 double i
= SCM_COMPLEX_IMAG (y
);
4774 if (fabs(r
) <= fabs(i
))
4777 double d
= i
* (1.0 + t
* t
);
4778 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4783 double d
= r
* (1.0 + t
* t
);
4784 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4788 else if (SCM_FRACTIONP (y
))
4789 /* a / b/c = ac / b */
4790 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4791 SCM_FRACTION_NUMERATOR (y
));
4793 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4795 else if (SCM_BIGP (x
))
4797 if (SCM_I_INUMP (y
))
4799 long int yy
= SCM_I_INUM (y
);
4802 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4803 scm_num_overflow (s_divide
);
4805 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4806 scm_remember_upto_here_1 (x
);
4807 return (sgn
== 0) ? scm_nan () : scm_inf ();
4814 /* FIXME: HMM, what are the relative performance issues here?
4815 We need to test. Is it faster on average to test
4816 divisible_p, then perform whichever operation, or is it
4817 faster to perform the integer div opportunistically and
4818 switch to real if there's a remainder? For now we take the
4819 middle ground: test, then if divisible, use the faster div
4822 long abs_yy
= yy
< 0 ? -yy
: yy
;
4823 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4827 SCM result
= scm_i_mkbig ();
4828 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4829 scm_remember_upto_here_1 (x
);
4831 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4832 return scm_i_normbig (result
);
4837 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4838 else return scm_i_make_ratio (x
, y
);
4842 else if (SCM_BIGP (y
))
4844 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4847 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4848 scm_num_overflow (s_divide
);
4850 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4851 scm_remember_upto_here_1 (x
);
4852 return (sgn
== 0) ? scm_nan () : scm_inf ();
4860 /* It's easily possible for the ratio x/y to fit a double
4861 but one or both x and y be too big to fit a double,
4862 hence the use of mpq_get_d rather than converting and
4865 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4866 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4867 return scm_from_double (mpq_get_d (q
));
4871 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4875 SCM result
= scm_i_mkbig ();
4876 mpz_divexact (SCM_I_BIG_MPZ (result
),
4879 scm_remember_upto_here_2 (x
, y
);
4880 return scm_i_normbig (result
);
4883 return scm_i_make_ratio (x
, y
);
4887 else if (SCM_REALP (y
))
4889 double yy
= SCM_REAL_VALUE (y
);
4890 #ifndef ALLOW_DIVIDE_BY_ZERO
4892 scm_num_overflow (s_divide
);
4895 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4897 else if (SCM_COMPLEXP (y
))
4899 a
= scm_i_big2dbl (x
);
4902 else if (SCM_FRACTIONP (y
))
4903 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4904 SCM_FRACTION_NUMERATOR (y
));
4906 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4908 else if (SCM_REALP (x
))
4910 double rx
= SCM_REAL_VALUE (x
);
4911 if (SCM_I_INUMP (y
))
4913 long int yy
= SCM_I_INUM (y
);
4914 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4916 scm_num_overflow (s_divide
);
4919 return scm_from_double (rx
/ (double) yy
);
4921 else if (SCM_BIGP (y
))
4923 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4924 scm_remember_upto_here_1 (y
);
4925 return scm_from_double (rx
/ dby
);
4927 else if (SCM_REALP (y
))
4929 double yy
= SCM_REAL_VALUE (y
);
4930 #ifndef ALLOW_DIVIDE_BY_ZERO
4932 scm_num_overflow (s_divide
);
4935 return scm_from_double (rx
/ yy
);
4937 else if (SCM_COMPLEXP (y
))
4942 else if (SCM_FRACTIONP (y
))
4943 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4945 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4947 else if (SCM_COMPLEXP (x
))
4949 double rx
= SCM_COMPLEX_REAL (x
);
4950 double ix
= SCM_COMPLEX_IMAG (x
);
4951 if (SCM_I_INUMP (y
))
4953 long int yy
= SCM_I_INUM (y
);
4954 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4956 scm_num_overflow (s_divide
);
4961 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4964 else if (SCM_BIGP (y
))
4966 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4967 scm_remember_upto_here_1 (y
);
4968 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4970 else if (SCM_REALP (y
))
4972 double yy
= SCM_REAL_VALUE (y
);
4973 #ifndef ALLOW_DIVIDE_BY_ZERO
4975 scm_num_overflow (s_divide
);
4978 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4980 else if (SCM_COMPLEXP (y
))
4982 double ry
= SCM_COMPLEX_REAL (y
);
4983 double iy
= SCM_COMPLEX_IMAG (y
);
4984 if (fabs(ry
) <= fabs(iy
))
4987 double d
= iy
* (1.0 + t
* t
);
4988 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4993 double d
= ry
* (1.0 + t
* t
);
4994 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4997 else if (SCM_FRACTIONP (y
))
4999 double yy
= scm_i_fraction2double (y
);
5000 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5003 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5005 else if (SCM_FRACTIONP (x
))
5007 if (SCM_I_INUMP (y
))
5009 long int yy
= SCM_I_INUM (y
);
5010 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5012 scm_num_overflow (s_divide
);
5015 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5016 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5018 else if (SCM_BIGP (y
))
5020 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5021 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5023 else if (SCM_REALP (y
))
5025 double yy
= SCM_REAL_VALUE (y
);
5026 #ifndef ALLOW_DIVIDE_BY_ZERO
5028 scm_num_overflow (s_divide
);
5031 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5033 else if (SCM_COMPLEXP (y
))
5035 a
= scm_i_fraction2double (x
);
5038 else if (SCM_FRACTIONP (y
))
5039 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5040 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5042 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5045 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5049 scm_divide (SCM x
, SCM y
)
5051 return scm_i_divide (x
, y
, 0);
5054 static SCM
scm_divide2real (SCM x
, SCM y
)
5056 return scm_i_divide (x
, y
, 1);
5062 scm_c_truncate (double x
)
5073 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5074 half-way case (ie. when x is an integer plus 0.5) going upwards.
5075 Then half-way cases are identified and adjusted down if the
5076 round-upwards didn't give the desired even integer.
5078 "plus_half == result" identifies a half-way case. If plus_half, which is
5079 x + 0.5, is an integer then x must be an integer plus 0.5.
5081 An odd "result" value is identified with result/2 != floor(result/2).
5082 This is done with plus_half, since that value is ready for use sooner in
5083 a pipelined cpu, and we're already requiring plus_half == result.
5085 Note however that we need to be careful when x is big and already an
5086 integer. In that case "x+0.5" may round to an adjacent integer, causing
5087 us to return such a value, incorrectly. For instance if the hardware is
5088 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5089 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5090 returned. Or if the hardware is in round-upwards mode, then other bigger
5091 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5092 representable value, 2^128+2^76 (or whatever), again incorrect.
5094 These bad roundings of x+0.5 are avoided by testing at the start whether
5095 x is already an integer. If it is then clearly that's the desired result
5096 already. And if it's not then the exponent must be small enough to allow
5097 an 0.5 to be represented, and hence added without a bad rounding. */
5100 scm_c_round (double x
)
5102 double plus_half
, result
;
5107 plus_half
= x
+ 0.5;
5108 result
= floor (plus_half
);
5109 /* Adjust so that the rounding is towards even. */
5110 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5115 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5117 "Round the number @var{x} towards zero.")
5118 #define FUNC_NAME s_scm_truncate_number
5120 if (scm_is_false (scm_negative_p (x
)))
5121 return scm_floor (x
);
5123 return scm_ceiling (x
);
5127 static SCM exactly_one_half
;
5129 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5131 "Round the number @var{x} towards the nearest integer. "
5132 "When it is exactly halfway between two integers, "
5133 "round towards the even one.")
5134 #define FUNC_NAME s_scm_round_number
5136 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5138 else if (SCM_REALP (x
))
5139 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5142 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5143 single quotient+remainder division then examining to see which way
5144 the rounding should go. */
5145 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5146 SCM result
= scm_floor (plus_half
);
5147 /* Adjust so that the rounding is towards even. */
5148 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5149 && scm_is_true (scm_odd_p (result
)))
5150 return scm_difference (result
, SCM_I_MAKINUM (1));
5157 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5159 "Round the number @var{x} towards minus infinity.")
5160 #define FUNC_NAME s_scm_floor
5162 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5164 else if (SCM_REALP (x
))
5165 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5166 else if (SCM_FRACTIONP (x
))
5168 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5169 SCM_FRACTION_DENOMINATOR (x
));
5170 if (scm_is_false (scm_negative_p (x
)))
5172 /* For positive x, rounding towards zero is correct. */
5177 /* For negative x, we need to return q-1 unless x is an
5178 integer. But fractions are never integer, per our
5180 return scm_difference (q
, SCM_I_MAKINUM (1));
5184 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5188 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5190 "Round the number @var{x} towards infinity.")
5191 #define FUNC_NAME s_scm_ceiling
5193 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5195 else if (SCM_REALP (x
))
5196 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5197 else if (SCM_FRACTIONP (x
))
5199 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5200 SCM_FRACTION_DENOMINATOR (x
));
5201 if (scm_is_false (scm_positive_p (x
)))
5203 /* For negative x, rounding towards zero is correct. */
5208 /* For positive x, we need to return q+1 unless x is an
5209 integer. But fractions are never integer, per our
5211 return scm_sum (q
, SCM_I_MAKINUM (1));
5215 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5219 /* sin/cos/tan/asin/acos/atan
5220 sinh/cosh/tanh/asinh/acosh/atanh
5221 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5222 Written by Jerry D. Hedden, (C) FSF.
5223 See the file `COPYING' for terms applying to this program. */
5225 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5227 "Return @var{x} raised to the power of @var{y}.")
5228 #define FUNC_NAME s_scm_expt
5230 if (!SCM_INEXACTP (y
) && scm_is_integer (y
))
5231 return scm_integer_expt (x
, y
);
5232 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5234 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5237 return scm_exp (scm_product (scm_log (x
), y
));
5241 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5243 "Compute the sine of @var{z}.")
5244 #define FUNC_NAME s_scm_sin
5246 if (scm_is_real (z
))
5247 return scm_from_double (sin (scm_to_double (z
)));
5248 else if (SCM_COMPLEXP (z
))
5250 x
= SCM_COMPLEX_REAL (z
);
5251 y
= SCM_COMPLEX_IMAG (z
);
5252 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5253 cos (x
) * sinh (y
));
5256 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5260 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5262 "Compute the cosine of @var{z}.")
5263 #define FUNC_NAME s_scm_cos
5265 if (scm_is_real (z
))
5266 return scm_from_double (cos (scm_to_double (z
)));
5267 else if (SCM_COMPLEXP (z
))
5269 x
= SCM_COMPLEX_REAL (z
);
5270 y
= SCM_COMPLEX_IMAG (z
);
5271 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5272 -sin (x
) * sinh (y
));
5275 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5279 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5281 "Compute the tangent of @var{z}.")
5282 #define FUNC_NAME s_scm_tan
5284 if (scm_is_real (z
))
5285 return scm_from_double (tan (scm_to_double (z
)));
5286 else if (SCM_COMPLEXP (z
))
5288 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5289 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5290 w
= cos (x
) + cosh (y
);
5291 #ifndef ALLOW_DIVIDE_BY_ZERO
5293 scm_num_overflow (s_scm_tan
);
5295 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5298 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5302 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5304 "Compute the hyperbolic sine of @var{z}.")
5305 #define FUNC_NAME s_scm_sinh
5307 if (scm_is_real (z
))
5308 return scm_from_double (sinh (scm_to_double (z
)));
5309 else if (SCM_COMPLEXP (z
))
5311 x
= SCM_COMPLEX_REAL (z
);
5312 y
= SCM_COMPLEX_IMAG (z
);
5313 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5314 cosh (x
) * sin (y
));
5317 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5321 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5323 "Compute the hyperbolic cosine of @var{z}.")
5324 #define FUNC_NAME s_scm_cosh
5326 if (scm_is_real (z
))
5327 return scm_from_double (cosh (scm_to_double (z
)));
5328 else if (SCM_COMPLEXP (z
))
5330 x
= SCM_COMPLEX_REAL (z
);
5331 y
= SCM_COMPLEX_IMAG (z
);
5332 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5333 sinh (x
) * sin (y
));
5336 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5340 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5342 "Compute the hyperbolic tangent of @var{z}.")
5343 #define FUNC_NAME s_scm_tanh
5345 if (scm_is_real (z
))
5346 return scm_from_double (tanh (scm_to_double (z
)));
5347 else if (SCM_COMPLEXP (z
))
5349 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5350 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5351 w
= cosh (x
) + cos (y
);
5352 #ifndef ALLOW_DIVIDE_BY_ZERO
5354 scm_num_overflow (s_scm_tanh
);
5356 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5359 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5363 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5365 "Compute the arc sine of @var{z}.")
5366 #define FUNC_NAME s_scm_asin
5368 if (scm_is_real (z
))
5370 double w
= scm_to_double (z
);
5371 if (w
>= -1.0 && w
<= 1.0)
5372 return scm_from_double (asin (w
));
5374 return scm_product (scm_c_make_rectangular (0, -1),
5375 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5377 else if (SCM_COMPLEXP (z
))
5379 x
= SCM_COMPLEX_REAL (z
);
5380 y
= SCM_COMPLEX_IMAG (z
);
5381 return scm_product (scm_c_make_rectangular (0, -1),
5382 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5385 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5389 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5391 "Compute the arc cosine of @var{z}.")
5392 #define FUNC_NAME s_scm_acos
5394 if (scm_is_real (z
))
5396 double w
= scm_to_double (z
);
5397 if (w
>= -1.0 && w
<= 1.0)
5398 return scm_from_double (acos (w
));
5400 return scm_sum (scm_from_double (acos (0.0)),
5401 scm_product (scm_c_make_rectangular (0, 1),
5402 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5404 else if (SCM_COMPLEXP (z
))
5406 x
= SCM_COMPLEX_REAL (z
);
5407 y
= SCM_COMPLEX_IMAG (z
);
5408 return scm_sum (scm_from_double (acos (0.0)),
5409 scm_product (scm_c_make_rectangular (0, 1),
5410 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5413 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5417 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5419 "With one argument, compute the arc tangent of @var{z}.\n"
5420 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5421 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5422 #define FUNC_NAME s_scm_atan
5426 if (scm_is_real (z
))
5427 return scm_from_double (atan (scm_to_double (z
)));
5428 else if (SCM_COMPLEXP (z
))
5431 v
= SCM_COMPLEX_REAL (z
);
5432 w
= SCM_COMPLEX_IMAG (z
);
5433 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5434 scm_c_make_rectangular (v
, w
+ 1.0))),
5435 scm_c_make_rectangular (0, 2));
5438 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5440 else if (scm_is_real (z
))
5442 if (scm_is_real (y
))
5443 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5445 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5448 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5452 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5454 "Compute the inverse hyperbolic sine of @var{z}.")
5455 #define FUNC_NAME s_scm_sys_asinh
5457 if (scm_is_real (z
))
5458 return scm_from_double (asinh (scm_to_double (z
)));
5459 else if (scm_is_number (z
))
5460 return scm_log (scm_sum (z
,
5461 scm_sqrt (scm_sum (scm_product (z
, z
),
5462 SCM_I_MAKINUM (1)))));
5464 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5468 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5470 "Compute the inverse hyperbolic cosine of @var{z}.")
5471 #define FUNC_NAME s_scm_sys_acosh
5473 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5474 return scm_from_double (acosh (scm_to_double (z
)));
5475 else if (scm_is_number (z
))
5476 return scm_log (scm_sum (z
,
5477 scm_sqrt (scm_difference (scm_product (z
, z
),
5478 SCM_I_MAKINUM (1)))));
5480 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5484 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5486 "Compute the inverse hyperbolic tangent of @var{z}.")
5487 #define FUNC_NAME s_scm_sys_atanh
5489 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5490 return scm_from_double (atanh (scm_to_double (z
)));
5491 else if (scm_is_number (z
))
5492 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5493 scm_difference (SCM_I_MAKINUM (1), z
))),
5496 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5501 scm_c_make_rectangular (double re
, double im
)
5504 return scm_from_double (re
);
5508 SCM_NEWSMOB (z
, scm_tc16_complex
,
5509 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5511 SCM_COMPLEX_REAL (z
) = re
;
5512 SCM_COMPLEX_IMAG (z
) = im
;
5517 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5518 (SCM real_part
, SCM imaginary_part
),
5519 "Return a complex number constructed of the given @var{real-part} "
5520 "and @var{imaginary-part} parts.")
5521 #define FUNC_NAME s_scm_make_rectangular
5523 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5524 SCM_ARG1
, FUNC_NAME
, "real");
5525 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5526 SCM_ARG2
, FUNC_NAME
, "real");
5527 return scm_c_make_rectangular (scm_to_double (real_part
),
5528 scm_to_double (imaginary_part
));
5533 scm_c_make_polar (double mag
, double ang
)
5537 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5538 use it on Glibc-based systems that have it (it's a GNU extension). See
5539 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5541 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5542 sincos (ang
, &s
, &c
);
5547 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5550 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5552 "Return the complex number @var{x} * e^(i * @var{y}).")
5553 #define FUNC_NAME s_scm_make_polar
5555 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5556 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5557 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5562 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5563 /* "Return the real part of the number @var{z}."
5566 scm_real_part (SCM z
)
5568 if (SCM_I_INUMP (z
))
5570 else if (SCM_BIGP (z
))
5572 else if (SCM_REALP (z
))
5574 else if (SCM_COMPLEXP (z
))
5575 return scm_from_double (SCM_COMPLEX_REAL (z
));
5576 else if (SCM_FRACTIONP (z
))
5579 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5583 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5584 /* "Return the imaginary part of the number @var{z}."
5587 scm_imag_part (SCM z
)
5589 if (SCM_I_INUMP (z
))
5591 else if (SCM_BIGP (z
))
5593 else if (SCM_REALP (z
))
5595 else if (SCM_COMPLEXP (z
))
5596 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5597 else if (SCM_FRACTIONP (z
))
5600 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5603 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5604 /* "Return the numerator of the number @var{z}."
5607 scm_numerator (SCM z
)
5609 if (SCM_I_INUMP (z
))
5611 else if (SCM_BIGP (z
))
5613 else if (SCM_FRACTIONP (z
))
5614 return SCM_FRACTION_NUMERATOR (z
);
5615 else if (SCM_REALP (z
))
5616 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5618 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5622 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5623 /* "Return the denominator of the number @var{z}."
5626 scm_denominator (SCM z
)
5628 if (SCM_I_INUMP (z
))
5629 return SCM_I_MAKINUM (1);
5630 else if (SCM_BIGP (z
))
5631 return SCM_I_MAKINUM (1);
5632 else if (SCM_FRACTIONP (z
))
5633 return SCM_FRACTION_DENOMINATOR (z
);
5634 else if (SCM_REALP (z
))
5635 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5637 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5640 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5641 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5642 * "@code{abs} for real arguments, but also allows complex numbers."
5645 scm_magnitude (SCM z
)
5647 if (SCM_I_INUMP (z
))
5649 long int zz
= SCM_I_INUM (z
);
5652 else if (SCM_POSFIXABLE (-zz
))
5653 return SCM_I_MAKINUM (-zz
);
5655 return scm_i_long2big (-zz
);
5657 else if (SCM_BIGP (z
))
5659 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5660 scm_remember_upto_here_1 (z
);
5662 return scm_i_clonebig (z
, 0);
5666 else if (SCM_REALP (z
))
5667 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5668 else if (SCM_COMPLEXP (z
))
5669 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5670 else if (SCM_FRACTIONP (z
))
5672 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5674 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5675 SCM_FRACTION_DENOMINATOR (z
));
5678 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5682 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5683 /* "Return the angle of the complex number @var{z}."
5688 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5689 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5690 But if atan2 follows the floating point rounding mode, then the value
5691 is not a constant. Maybe it'd be close enough though. */
5692 if (SCM_I_INUMP (z
))
5694 if (SCM_I_INUM (z
) >= 0)
5697 return scm_from_double (atan2 (0.0, -1.0));
5699 else if (SCM_BIGP (z
))
5701 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5702 scm_remember_upto_here_1 (z
);
5704 return scm_from_double (atan2 (0.0, -1.0));
5708 else if (SCM_REALP (z
))
5710 if (SCM_REAL_VALUE (z
) >= 0)
5713 return scm_from_double (atan2 (0.0, -1.0));
5715 else if (SCM_COMPLEXP (z
))
5716 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5717 else if (SCM_FRACTIONP (z
))
5719 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5721 else return scm_from_double (atan2 (0.0, -1.0));
5724 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5728 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5729 /* Convert the number @var{x} to its inexact representation.\n"
5732 scm_exact_to_inexact (SCM z
)
5734 if (SCM_I_INUMP (z
))
5735 return scm_from_double ((double) SCM_I_INUM (z
));
5736 else if (SCM_BIGP (z
))
5737 return scm_from_double (scm_i_big2dbl (z
));
5738 else if (SCM_FRACTIONP (z
))
5739 return scm_from_double (scm_i_fraction2double (z
));
5740 else if (SCM_INEXACTP (z
))
5743 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5747 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5749 "Return an exact number that is numerically closest to @var{z}.")
5750 #define FUNC_NAME s_scm_inexact_to_exact
5752 if (SCM_I_INUMP (z
))
5754 else if (SCM_BIGP (z
))
5756 else if (SCM_REALP (z
))
5758 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5759 SCM_OUT_OF_RANGE (1, z
);
5766 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5767 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5768 scm_i_mpz2num (mpq_denref (frac
)));
5770 /* When scm_i_make_ratio throws, we leak the memory allocated
5777 else if (SCM_FRACTIONP (z
))
5780 SCM_WRONG_TYPE_ARG (1, z
);
5784 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5786 "Returns the @emph{simplest} rational number differing\n"
5787 "from @var{x} by no more than @var{eps}.\n"
5789 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5790 "exact result when both its arguments are exact. Thus, you might need\n"
5791 "to use @code{inexact->exact} on the arguments.\n"
5794 "(rationalize (inexact->exact 1.2) 1/100)\n"
5797 #define FUNC_NAME s_scm_rationalize
5799 if (SCM_I_INUMP (x
))
5801 else if (SCM_BIGP (x
))
5803 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5805 /* Use continued fractions to find closest ratio. All
5806 arithmetic is done with exact numbers.
5809 SCM ex
= scm_inexact_to_exact (x
);
5810 SCM int_part
= scm_floor (ex
);
5811 SCM tt
= SCM_I_MAKINUM (1);
5812 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5813 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5817 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5820 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5821 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5823 /* We stop after a million iterations just to be absolutely sure
5824 that we don't go into an infinite loop. The process normally
5825 converges after less than a dozen iterations.
5828 eps
= scm_abs (eps
);
5829 while (++i
< 1000000)
5831 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5832 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5833 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5835 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5836 eps
))) /* abs(x-a/b) <= eps */
5838 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5839 if (scm_is_false (scm_exact_p (x
))
5840 || scm_is_false (scm_exact_p (eps
)))
5841 return scm_exact_to_inexact (res
);
5845 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5847 tt
= scm_floor (rx
); /* tt = floor (rx) */
5853 scm_num_overflow (s_scm_rationalize
);
5856 SCM_WRONG_TYPE_ARG (1, x
);
5860 /* conversion functions */
5863 scm_is_integer (SCM val
)
5865 return scm_is_true (scm_integer_p (val
));
5869 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5871 if (SCM_I_INUMP (val
))
5873 scm_t_signed_bits n
= SCM_I_INUM (val
);
5874 return n
>= min
&& n
<= max
;
5876 else if (SCM_BIGP (val
))
5878 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5880 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5882 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5884 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5885 return n
>= min
&& n
<= max
;
5895 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5896 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5899 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5900 SCM_I_BIG_MPZ (val
));
5902 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5914 return n
>= min
&& n
<= max
;
5922 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5924 if (SCM_I_INUMP (val
))
5926 scm_t_signed_bits n
= SCM_I_INUM (val
);
5927 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5929 else if (SCM_BIGP (val
))
5931 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5933 else if (max
<= ULONG_MAX
)
5935 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5937 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5938 return n
>= min
&& n
<= max
;
5948 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5951 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5952 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5955 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5956 SCM_I_BIG_MPZ (val
));
5958 return n
>= min
&& n
<= max
;
5966 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5968 scm_error (scm_out_of_range_key
,
5970 "Value out of range ~S to ~S: ~S",
5971 scm_list_3 (min
, max
, bad_val
),
5972 scm_list_1 (bad_val
));
5975 #define TYPE scm_t_intmax
5976 #define TYPE_MIN min
5977 #define TYPE_MAX max
5978 #define SIZEOF_TYPE 0
5979 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5980 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5981 #include "libguile/conv-integer.i.c"
5983 #define TYPE scm_t_uintmax
5984 #define TYPE_MIN min
5985 #define TYPE_MAX max
5986 #define SIZEOF_TYPE 0
5987 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5988 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5989 #include "libguile/conv-uinteger.i.c"
5991 #define TYPE scm_t_int8
5992 #define TYPE_MIN SCM_T_INT8_MIN
5993 #define TYPE_MAX SCM_T_INT8_MAX
5994 #define SIZEOF_TYPE 1
5995 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5996 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5997 #include "libguile/conv-integer.i.c"
5999 #define TYPE scm_t_uint8
6001 #define TYPE_MAX SCM_T_UINT8_MAX
6002 #define SIZEOF_TYPE 1
6003 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6004 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6005 #include "libguile/conv-uinteger.i.c"
6007 #define TYPE scm_t_int16
6008 #define TYPE_MIN SCM_T_INT16_MIN
6009 #define TYPE_MAX SCM_T_INT16_MAX
6010 #define SIZEOF_TYPE 2
6011 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6012 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6013 #include "libguile/conv-integer.i.c"
6015 #define TYPE scm_t_uint16
6017 #define TYPE_MAX SCM_T_UINT16_MAX
6018 #define SIZEOF_TYPE 2
6019 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6020 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6021 #include "libguile/conv-uinteger.i.c"
6023 #define TYPE scm_t_int32
6024 #define TYPE_MIN SCM_T_INT32_MIN
6025 #define TYPE_MAX SCM_T_INT32_MAX
6026 #define SIZEOF_TYPE 4
6027 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6028 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6029 #include "libguile/conv-integer.i.c"
6031 #define TYPE scm_t_uint32
6033 #define TYPE_MAX SCM_T_UINT32_MAX
6034 #define SIZEOF_TYPE 4
6035 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6036 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6037 #include "libguile/conv-uinteger.i.c"
6039 #define TYPE scm_t_wchar
6040 #define TYPE_MIN (scm_t_int32)-1
6041 #define TYPE_MAX (scm_t_int32)0x10ffff
6042 #define SIZEOF_TYPE 4
6043 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6044 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6045 #include "libguile/conv-integer.i.c"
6047 #if SCM_HAVE_T_INT64
6049 #define TYPE scm_t_int64
6050 #define TYPE_MIN SCM_T_INT64_MIN
6051 #define TYPE_MAX SCM_T_INT64_MAX
6052 #define SIZEOF_TYPE 8
6053 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6054 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6055 #include "libguile/conv-integer.i.c"
6057 #define TYPE scm_t_uint64
6059 #define TYPE_MAX SCM_T_UINT64_MAX
6060 #define SIZEOF_TYPE 8
6061 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6062 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6063 #include "libguile/conv-uinteger.i.c"
6068 scm_to_mpz (SCM val
, mpz_t rop
)
6070 if (SCM_I_INUMP (val
))
6071 mpz_set_si (rop
, SCM_I_INUM (val
));
6072 else if (SCM_BIGP (val
))
6073 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6075 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6079 scm_from_mpz (mpz_t val
)
6081 return scm_i_mpz2num (val
);
6085 scm_is_real (SCM val
)
6087 return scm_is_true (scm_real_p (val
));
6091 scm_is_rational (SCM val
)
6093 return scm_is_true (scm_rational_p (val
));
6097 scm_to_double (SCM val
)
6099 if (SCM_I_INUMP (val
))
6100 return SCM_I_INUM (val
);
6101 else if (SCM_BIGP (val
))
6102 return scm_i_big2dbl (val
);
6103 else if (SCM_FRACTIONP (val
))
6104 return scm_i_fraction2double (val
);
6105 else if (SCM_REALP (val
))
6106 return SCM_REAL_VALUE (val
);
6108 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6112 scm_from_double (double val
)
6114 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
6115 SCM_REAL_VALUE (z
) = val
;
6119 #if SCM_ENABLE_DISCOURAGED == 1
6122 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6126 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6130 scm_out_of_range (NULL
, num
);
6133 return scm_to_double (num
);
6137 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6141 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6145 scm_out_of_range (NULL
, num
);
6148 return scm_to_double (num
);
6154 scm_is_complex (SCM val
)
6156 return scm_is_true (scm_complex_p (val
));
6160 scm_c_real_part (SCM z
)
6162 if (SCM_COMPLEXP (z
))
6163 return SCM_COMPLEX_REAL (z
);
6166 /* Use the scm_real_part to get proper error checking and
6169 return scm_to_double (scm_real_part (z
));
6174 scm_c_imag_part (SCM z
)
6176 if (SCM_COMPLEXP (z
))
6177 return SCM_COMPLEX_IMAG (z
);
6180 /* Use the scm_imag_part to get proper error checking and
6181 dispatching. The result will almost always be 0.0, but not
6184 return scm_to_double (scm_imag_part (z
));
6189 scm_c_magnitude (SCM z
)
6191 return scm_to_double (scm_magnitude (z
));
6197 return scm_to_double (scm_angle (z
));
6201 scm_is_number (SCM z
)
6203 return scm_is_true (scm_number_p (z
));
6207 /* In the following functions we dispatch to the real-arg funcs like log()
6208 when we know the arg is real, instead of just handing everything to
6209 clog() for instance. This is in case clog() doesn't optimize for a
6210 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6211 well use it to go straight to the applicable C func. */
6213 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6215 "Return the natural logarithm of @var{z}.")
6216 #define FUNC_NAME s_scm_log
6218 if (SCM_COMPLEXP (z
))
6220 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6221 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6223 double re
= SCM_COMPLEX_REAL (z
);
6224 double im
= SCM_COMPLEX_IMAG (z
);
6225 return scm_c_make_rectangular (log (hypot (re
, im
)),
6231 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6232 although the value itself overflows. */
6233 double re
= scm_to_double (z
);
6234 double l
= log (fabs (re
));
6236 return scm_from_double (l
);
6238 return scm_c_make_rectangular (l
, M_PI
);
6244 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6246 "Return the base 10 logarithm of @var{z}.")
6247 #define FUNC_NAME s_scm_log10
6249 if (SCM_COMPLEXP (z
))
6251 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6252 clog() and a multiply by M_LOG10E, rather than the fallback
6253 log10+hypot+atan2.) */
6254 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6255 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6257 double re
= SCM_COMPLEX_REAL (z
);
6258 double im
= SCM_COMPLEX_IMAG (z
);
6259 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6260 M_LOG10E
* atan2 (im
, re
));
6265 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6266 although the value itself overflows. */
6267 double re
= scm_to_double (z
);
6268 double l
= log10 (fabs (re
));
6270 return scm_from_double (l
);
6272 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6278 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6280 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6281 "base of natural logarithms (2.71828@dots{}).")
6282 #define FUNC_NAME s_scm_exp
6284 if (SCM_COMPLEXP (z
))
6286 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6287 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6289 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6290 SCM_COMPLEX_IMAG (z
));
6295 /* When z is a negative bignum the conversion to double overflows,
6296 giving -infinity, but that's ok, the exp is still 0.0. */
6297 return scm_from_double (exp (scm_to_double (z
)));
6303 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6305 "Return the square root of @var{z}. Of the two possible roots\n"
6306 "(positive and negative), the one with the a positive real part\n"
6307 "is returned, or if that's zero then a positive imaginary part.\n"
6311 "(sqrt 9.0) @result{} 3.0\n"
6312 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6313 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6314 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6316 #define FUNC_NAME s_scm_sqrt
6318 if (SCM_COMPLEXP (x
))
6320 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6321 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6323 double re
= SCM_COMPLEX_REAL (x
);
6324 double im
= SCM_COMPLEX_IMAG (x
);
6325 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6326 0.5 * atan2 (im
, re
));
6331 double xx
= scm_to_double (x
);
6333 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6335 return scm_from_double (sqrt (xx
));
6347 mpz_init_set_si (z_negative_one
, -1);
6349 /* It may be possible to tune the performance of some algorithms by using
6350 * the following constants to avoid the creation of bignums. Please, before
6351 * using these values, remember the two rules of program optimization:
6352 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6353 scm_c_define ("most-positive-fixnum",
6354 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6355 scm_c_define ("most-negative-fixnum",
6356 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6358 scm_add_feature ("complex");
6359 scm_add_feature ("inexact");
6360 scm_flo0
= scm_from_double (0.0);
6362 /* determine floating point precision */
6363 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6365 init_dblprec(&scm_dblprec
[i
-2],i
);
6366 init_fx_radix(fx_per_radix
[i
-2],i
);
6369 /* hard code precision for base 10 if the preprocessor tells us to... */
6370 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6373 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6374 SCM_I_MAKINUM (2)));
6375 #include "libguile/numbers.x"