6fc4c2abf9677cd560c98553b576bf93370e4ad3
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_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
3995 (SCM x
, SCM y
, SCM rest
),
3996 "Return the sum of all parameter values. Return 0 if called without\n"
3998 #define FUNC_NAME s_scm_i_sum
4000 while (!scm_is_null (rest
))
4001 { x
= scm_sum (x
, y
);
4003 rest
= scm_cdr (rest
);
4005 return scm_sum (x
, y
);
4009 #define s_sum s_scm_i_sum
4010 #define g_sum g_scm_i_sum
4013 scm_sum (SCM x
, SCM y
)
4015 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4017 if (SCM_NUMBERP (x
)) return x
;
4018 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4019 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4022 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4024 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4026 long xx
= SCM_I_INUM (x
);
4027 long yy
= SCM_I_INUM (y
);
4028 long int z
= xx
+ yy
;
4029 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4031 else if (SCM_BIGP (y
))
4036 else if (SCM_REALP (y
))
4038 long int xx
= SCM_I_INUM (x
);
4039 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4041 else if (SCM_COMPLEXP (y
))
4043 long int xx
= SCM_I_INUM (x
);
4044 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4045 SCM_COMPLEX_IMAG (y
));
4047 else if (SCM_FRACTIONP (y
))
4048 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4049 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4050 SCM_FRACTION_DENOMINATOR (y
));
4052 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4053 } else if (SCM_BIGP (x
))
4055 if (SCM_I_INUMP (y
))
4060 inum
= SCM_I_INUM (y
);
4063 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4066 SCM result
= scm_i_mkbig ();
4067 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4068 scm_remember_upto_here_1 (x
);
4069 /* we know the result will have to be a bignum */
4072 return scm_i_normbig (result
);
4076 SCM result
= scm_i_mkbig ();
4077 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4078 scm_remember_upto_here_1 (x
);
4079 /* we know the result will have to be a bignum */
4082 return scm_i_normbig (result
);
4085 else if (SCM_BIGP (y
))
4087 SCM result
= scm_i_mkbig ();
4088 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4089 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4090 mpz_add (SCM_I_BIG_MPZ (result
),
4093 scm_remember_upto_here_2 (x
, y
);
4094 /* we know the result will have to be a bignum */
4097 return scm_i_normbig (result
);
4099 else if (SCM_REALP (y
))
4101 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4102 scm_remember_upto_here_1 (x
);
4103 return scm_from_double (result
);
4105 else if (SCM_COMPLEXP (y
))
4107 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4108 + SCM_COMPLEX_REAL (y
));
4109 scm_remember_upto_here_1 (x
);
4110 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4112 else if (SCM_FRACTIONP (y
))
4113 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4114 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4115 SCM_FRACTION_DENOMINATOR (y
));
4117 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4119 else if (SCM_REALP (x
))
4121 if (SCM_I_INUMP (y
))
4122 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4123 else if (SCM_BIGP (y
))
4125 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4126 scm_remember_upto_here_1 (y
);
4127 return scm_from_double (result
);
4129 else if (SCM_REALP (y
))
4130 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4131 else if (SCM_COMPLEXP (y
))
4132 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4133 SCM_COMPLEX_IMAG (y
));
4134 else if (SCM_FRACTIONP (y
))
4135 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4137 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4139 else if (SCM_COMPLEXP (x
))
4141 if (SCM_I_INUMP (y
))
4142 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4143 SCM_COMPLEX_IMAG (x
));
4144 else if (SCM_BIGP (y
))
4146 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4147 + SCM_COMPLEX_REAL (x
));
4148 scm_remember_upto_here_1 (y
);
4149 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4151 else if (SCM_REALP (y
))
4152 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4153 SCM_COMPLEX_IMAG (x
));
4154 else if (SCM_COMPLEXP (y
))
4155 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4156 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4157 else if (SCM_FRACTIONP (y
))
4158 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4159 SCM_COMPLEX_IMAG (x
));
4161 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4163 else if (SCM_FRACTIONP (x
))
4165 if (SCM_I_INUMP (y
))
4166 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4167 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4168 SCM_FRACTION_DENOMINATOR (x
));
4169 else if (SCM_BIGP (y
))
4170 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4171 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4172 SCM_FRACTION_DENOMINATOR (x
));
4173 else if (SCM_REALP (y
))
4174 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4175 else if (SCM_COMPLEXP (y
))
4176 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4177 SCM_COMPLEX_IMAG (y
));
4178 else if (SCM_FRACTIONP (y
))
4179 /* a/b + c/d = (ad + bc) / bd */
4180 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4181 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4182 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4184 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4187 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4191 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4193 "Return @math{@var{x}+1}.")
4194 #define FUNC_NAME s_scm_oneplus
4196 return scm_sum (x
, SCM_I_MAKINUM (1));
4201 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4202 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4203 * the sum of all but the first argument are subtracted from the first
4205 #define FUNC_NAME s_difference
4207 scm_difference (SCM x
, SCM y
)
4209 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4212 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4214 if (SCM_I_INUMP (x
))
4216 long xx
= -SCM_I_INUM (x
);
4217 if (SCM_FIXABLE (xx
))
4218 return SCM_I_MAKINUM (xx
);
4220 return scm_i_long2big (xx
);
4222 else if (SCM_BIGP (x
))
4223 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4224 bignum, but negating that gives a fixnum. */
4225 return scm_i_normbig (scm_i_clonebig (x
, 0));
4226 else if (SCM_REALP (x
))
4227 return scm_from_double (-SCM_REAL_VALUE (x
));
4228 else if (SCM_COMPLEXP (x
))
4229 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4230 -SCM_COMPLEX_IMAG (x
));
4231 else if (SCM_FRACTIONP (x
))
4232 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4233 SCM_FRACTION_DENOMINATOR (x
));
4235 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4238 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4240 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4242 long int xx
= SCM_I_INUM (x
);
4243 long int yy
= SCM_I_INUM (y
);
4244 long int z
= xx
- yy
;
4245 if (SCM_FIXABLE (z
))
4246 return SCM_I_MAKINUM (z
);
4248 return scm_i_long2big (z
);
4250 else if (SCM_BIGP (y
))
4252 /* inum-x - big-y */
4253 long xx
= SCM_I_INUM (x
);
4256 return scm_i_clonebig (y
, 0);
4259 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4260 SCM result
= scm_i_mkbig ();
4263 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4266 /* x - y == -(y + -x) */
4267 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4268 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4270 scm_remember_upto_here_1 (y
);
4272 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4273 /* we know the result will have to be a bignum */
4276 return scm_i_normbig (result
);
4279 else if (SCM_REALP (y
))
4281 long int xx
= SCM_I_INUM (x
);
4282 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4284 else if (SCM_COMPLEXP (y
))
4286 long int xx
= SCM_I_INUM (x
);
4287 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4288 - SCM_COMPLEX_IMAG (y
));
4290 else if (SCM_FRACTIONP (y
))
4291 /* a - b/c = (ac - b) / c */
4292 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4293 SCM_FRACTION_NUMERATOR (y
)),
4294 SCM_FRACTION_DENOMINATOR (y
));
4296 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4298 else if (SCM_BIGP (x
))
4300 if (SCM_I_INUMP (y
))
4302 /* big-x - inum-y */
4303 long yy
= SCM_I_INUM (y
);
4304 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4306 scm_remember_upto_here_1 (x
);
4308 return (SCM_FIXABLE (-yy
) ?
4309 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4312 SCM result
= scm_i_mkbig ();
4315 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4317 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4318 scm_remember_upto_here_1 (x
);
4320 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4321 /* we know the result will have to be a bignum */
4324 return scm_i_normbig (result
);
4327 else if (SCM_BIGP (y
))
4329 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4330 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4331 SCM result
= scm_i_mkbig ();
4332 mpz_sub (SCM_I_BIG_MPZ (result
),
4335 scm_remember_upto_here_2 (x
, y
);
4336 /* we know the result will have to be a bignum */
4337 if ((sgn_x
== 1) && (sgn_y
== -1))
4339 if ((sgn_x
== -1) && (sgn_y
== 1))
4341 return scm_i_normbig (result
);
4343 else if (SCM_REALP (y
))
4345 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4346 scm_remember_upto_here_1 (x
);
4347 return scm_from_double (result
);
4349 else if (SCM_COMPLEXP (y
))
4351 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4352 - SCM_COMPLEX_REAL (y
));
4353 scm_remember_upto_here_1 (x
);
4354 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4356 else if (SCM_FRACTIONP (y
))
4357 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4358 SCM_FRACTION_NUMERATOR (y
)),
4359 SCM_FRACTION_DENOMINATOR (y
));
4360 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4362 else if (SCM_REALP (x
))
4364 if (SCM_I_INUMP (y
))
4365 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4366 else if (SCM_BIGP (y
))
4368 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4369 scm_remember_upto_here_1 (x
);
4370 return scm_from_double (result
);
4372 else if (SCM_REALP (y
))
4373 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4374 else if (SCM_COMPLEXP (y
))
4375 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4376 -SCM_COMPLEX_IMAG (y
));
4377 else if (SCM_FRACTIONP (y
))
4378 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4380 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4382 else if (SCM_COMPLEXP (x
))
4384 if (SCM_I_INUMP (y
))
4385 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4386 SCM_COMPLEX_IMAG (x
));
4387 else if (SCM_BIGP (y
))
4389 double real_part
= (SCM_COMPLEX_REAL (x
)
4390 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4391 scm_remember_upto_here_1 (x
);
4392 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4394 else if (SCM_REALP (y
))
4395 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4396 SCM_COMPLEX_IMAG (x
));
4397 else if (SCM_COMPLEXP (y
))
4398 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4399 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4400 else if (SCM_FRACTIONP (y
))
4401 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4402 SCM_COMPLEX_IMAG (x
));
4404 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4406 else if (SCM_FRACTIONP (x
))
4408 if (SCM_I_INUMP (y
))
4409 /* a/b - c = (a - cb) / b */
4410 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4411 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4412 SCM_FRACTION_DENOMINATOR (x
));
4413 else if (SCM_BIGP (y
))
4414 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4415 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4416 SCM_FRACTION_DENOMINATOR (x
));
4417 else if (SCM_REALP (y
))
4418 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4419 else if (SCM_COMPLEXP (y
))
4420 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4421 -SCM_COMPLEX_IMAG (y
));
4422 else if (SCM_FRACTIONP (y
))
4423 /* a/b - c/d = (ad - bc) / bd */
4424 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4425 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4426 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4428 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4431 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4436 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4438 "Return @math{@var{x}-1}.")
4439 #define FUNC_NAME s_scm_oneminus
4441 return scm_difference (x
, SCM_I_MAKINUM (1));
4446 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4447 /* "Return the product of all arguments. If called without arguments,\n"
4451 scm_product (SCM x
, SCM y
)
4453 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4456 return SCM_I_MAKINUM (1L);
4457 else if (SCM_NUMBERP (x
))
4460 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4463 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4468 xx
= SCM_I_INUM (x
);
4472 case 0: return x
; break;
4473 case 1: return y
; break;
4476 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4478 long yy
= SCM_I_INUM (y
);
4480 SCM k
= SCM_I_MAKINUM (kk
);
4481 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4485 SCM result
= scm_i_long2big (xx
);
4486 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4487 return scm_i_normbig (result
);
4490 else if (SCM_BIGP (y
))
4492 SCM result
= scm_i_mkbig ();
4493 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4494 scm_remember_upto_here_1 (y
);
4497 else if (SCM_REALP (y
))
4498 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4499 else if (SCM_COMPLEXP (y
))
4500 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4501 xx
* SCM_COMPLEX_IMAG (y
));
4502 else if (SCM_FRACTIONP (y
))
4503 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4504 SCM_FRACTION_DENOMINATOR (y
));
4506 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4508 else if (SCM_BIGP (x
))
4510 if (SCM_I_INUMP (y
))
4515 else if (SCM_BIGP (y
))
4517 SCM result
= scm_i_mkbig ();
4518 mpz_mul (SCM_I_BIG_MPZ (result
),
4521 scm_remember_upto_here_2 (x
, y
);
4524 else if (SCM_REALP (y
))
4526 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4527 scm_remember_upto_here_1 (x
);
4528 return scm_from_double (result
);
4530 else if (SCM_COMPLEXP (y
))
4532 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4533 scm_remember_upto_here_1 (x
);
4534 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4535 z
* SCM_COMPLEX_IMAG (y
));
4537 else if (SCM_FRACTIONP (y
))
4538 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4539 SCM_FRACTION_DENOMINATOR (y
));
4541 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4543 else if (SCM_REALP (x
))
4545 if (SCM_I_INUMP (y
))
4547 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4548 if (scm_is_eq (y
, SCM_INUM0
))
4550 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4552 else if (SCM_BIGP (y
))
4554 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4555 scm_remember_upto_here_1 (y
);
4556 return scm_from_double (result
);
4558 else if (SCM_REALP (y
))
4559 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4560 else if (SCM_COMPLEXP (y
))
4561 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4562 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4563 else if (SCM_FRACTIONP (y
))
4564 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4566 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4568 else if (SCM_COMPLEXP (x
))
4570 if (SCM_I_INUMP (y
))
4572 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4573 if (scm_is_eq (y
, SCM_INUM0
))
4575 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4576 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4578 else if (SCM_BIGP (y
))
4580 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4581 scm_remember_upto_here_1 (y
);
4582 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4583 z
* SCM_COMPLEX_IMAG (x
));
4585 else if (SCM_REALP (y
))
4586 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4587 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4588 else if (SCM_COMPLEXP (y
))
4590 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4591 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4592 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4593 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4595 else if (SCM_FRACTIONP (y
))
4597 double yy
= scm_i_fraction2double (y
);
4598 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4599 yy
* SCM_COMPLEX_IMAG (x
));
4602 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4604 else if (SCM_FRACTIONP (x
))
4606 if (SCM_I_INUMP (y
))
4607 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4608 SCM_FRACTION_DENOMINATOR (x
));
4609 else if (SCM_BIGP (y
))
4610 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4611 SCM_FRACTION_DENOMINATOR (x
));
4612 else if (SCM_REALP (y
))
4613 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4614 else if (SCM_COMPLEXP (y
))
4616 double xx
= scm_i_fraction2double (x
);
4617 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4618 xx
* SCM_COMPLEX_IMAG (y
));
4620 else if (SCM_FRACTIONP (y
))
4621 /* a/b * c/d = ac / bd */
4622 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4623 SCM_FRACTION_NUMERATOR (y
)),
4624 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4625 SCM_FRACTION_DENOMINATOR (y
)));
4627 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4630 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4633 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4634 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4635 #define ALLOW_DIVIDE_BY_ZERO
4636 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4639 /* The code below for complex division is adapted from the GNU
4640 libstdc++, which adapted it from f2c's libF77, and is subject to
4643 /****************************************************************
4644 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4646 Permission to use, copy, modify, and distribute this software
4647 and its documentation for any purpose and without fee is hereby
4648 granted, provided that the above copyright notice appear in all
4649 copies and that both that the copyright notice and this
4650 permission notice and warranty disclaimer appear in supporting
4651 documentation, and that the names of AT&T Bell Laboratories or
4652 Bellcore or any of their entities not be used in advertising or
4653 publicity pertaining to distribution of the software without
4654 specific, written prior permission.
4656 AT&T and Bellcore disclaim all warranties with regard to this
4657 software, including all implied warranties of merchantability
4658 and fitness. In no event shall AT&T or Bellcore be liable for
4659 any special, indirect or consequential damages or any damages
4660 whatsoever resulting from loss of use, data or profits, whether
4661 in an action of contract, negligence or other tortious action,
4662 arising out of or in connection with the use or performance of
4664 ****************************************************************/
4666 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4667 /* Divide the first argument by the product of the remaining
4668 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4670 #define FUNC_NAME s_divide
4672 scm_i_divide (SCM x
, SCM y
, int inexact
)
4676 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4679 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4680 else if (SCM_I_INUMP (x
))
4682 long xx
= SCM_I_INUM (x
);
4683 if (xx
== 1 || xx
== -1)
4685 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4687 scm_num_overflow (s_divide
);
4692 return scm_from_double (1.0 / (double) xx
);
4693 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4696 else if (SCM_BIGP (x
))
4699 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4700 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4702 else if (SCM_REALP (x
))
4704 double xx
= SCM_REAL_VALUE (x
);
4705 #ifndef ALLOW_DIVIDE_BY_ZERO
4707 scm_num_overflow (s_divide
);
4710 return scm_from_double (1.0 / xx
);
4712 else if (SCM_COMPLEXP (x
))
4714 double r
= SCM_COMPLEX_REAL (x
);
4715 double i
= SCM_COMPLEX_IMAG (x
);
4716 if (fabs(r
) <= fabs(i
))
4719 double d
= i
* (1.0 + t
* t
);
4720 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4725 double d
= r
* (1.0 + t
* t
);
4726 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4729 else if (SCM_FRACTIONP (x
))
4730 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4731 SCM_FRACTION_NUMERATOR (x
));
4733 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4736 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4738 long xx
= SCM_I_INUM (x
);
4739 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4741 long yy
= SCM_I_INUM (y
);
4744 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4745 scm_num_overflow (s_divide
);
4747 return scm_from_double ((double) xx
/ (double) yy
);
4750 else if (xx
% yy
!= 0)
4753 return scm_from_double ((double) xx
/ (double) yy
);
4754 else return scm_i_make_ratio (x
, y
);
4759 if (SCM_FIXABLE (z
))
4760 return SCM_I_MAKINUM (z
);
4762 return scm_i_long2big (z
);
4765 else if (SCM_BIGP (y
))
4768 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4769 else return scm_i_make_ratio (x
, y
);
4771 else if (SCM_REALP (y
))
4773 double yy
= SCM_REAL_VALUE (y
);
4774 #ifndef ALLOW_DIVIDE_BY_ZERO
4776 scm_num_overflow (s_divide
);
4779 return scm_from_double ((double) xx
/ yy
);
4781 else if (SCM_COMPLEXP (y
))
4784 complex_div
: /* y _must_ be a complex number */
4786 double r
= SCM_COMPLEX_REAL (y
);
4787 double i
= SCM_COMPLEX_IMAG (y
);
4788 if (fabs(r
) <= fabs(i
))
4791 double d
= i
* (1.0 + t
* t
);
4792 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4797 double d
= r
* (1.0 + t
* t
);
4798 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4802 else if (SCM_FRACTIONP (y
))
4803 /* a / b/c = ac / b */
4804 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4805 SCM_FRACTION_NUMERATOR (y
));
4807 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4809 else if (SCM_BIGP (x
))
4811 if (SCM_I_INUMP (y
))
4813 long int yy
= SCM_I_INUM (y
);
4816 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4817 scm_num_overflow (s_divide
);
4819 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4820 scm_remember_upto_here_1 (x
);
4821 return (sgn
== 0) ? scm_nan () : scm_inf ();
4828 /* FIXME: HMM, what are the relative performance issues here?
4829 We need to test. Is it faster on average to test
4830 divisible_p, then perform whichever operation, or is it
4831 faster to perform the integer div opportunistically and
4832 switch to real if there's a remainder? For now we take the
4833 middle ground: test, then if divisible, use the faster div
4836 long abs_yy
= yy
< 0 ? -yy
: yy
;
4837 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4841 SCM result
= scm_i_mkbig ();
4842 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4843 scm_remember_upto_here_1 (x
);
4845 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4846 return scm_i_normbig (result
);
4851 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4852 else return scm_i_make_ratio (x
, y
);
4856 else if (SCM_BIGP (y
))
4858 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4861 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4862 scm_num_overflow (s_divide
);
4864 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4865 scm_remember_upto_here_1 (x
);
4866 return (sgn
== 0) ? scm_nan () : scm_inf ();
4874 /* It's easily possible for the ratio x/y to fit a double
4875 but one or both x and y be too big to fit a double,
4876 hence the use of mpq_get_d rather than converting and
4879 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4880 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4881 return scm_from_double (mpq_get_d (q
));
4885 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4889 SCM result
= scm_i_mkbig ();
4890 mpz_divexact (SCM_I_BIG_MPZ (result
),
4893 scm_remember_upto_here_2 (x
, y
);
4894 return scm_i_normbig (result
);
4897 return scm_i_make_ratio (x
, y
);
4901 else if (SCM_REALP (y
))
4903 double yy
= SCM_REAL_VALUE (y
);
4904 #ifndef ALLOW_DIVIDE_BY_ZERO
4906 scm_num_overflow (s_divide
);
4909 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4911 else if (SCM_COMPLEXP (y
))
4913 a
= scm_i_big2dbl (x
);
4916 else if (SCM_FRACTIONP (y
))
4917 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4918 SCM_FRACTION_NUMERATOR (y
));
4920 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4922 else if (SCM_REALP (x
))
4924 double rx
= SCM_REAL_VALUE (x
);
4925 if (SCM_I_INUMP (y
))
4927 long int yy
= SCM_I_INUM (y
);
4928 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4930 scm_num_overflow (s_divide
);
4933 return scm_from_double (rx
/ (double) yy
);
4935 else if (SCM_BIGP (y
))
4937 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4938 scm_remember_upto_here_1 (y
);
4939 return scm_from_double (rx
/ dby
);
4941 else if (SCM_REALP (y
))
4943 double yy
= SCM_REAL_VALUE (y
);
4944 #ifndef ALLOW_DIVIDE_BY_ZERO
4946 scm_num_overflow (s_divide
);
4949 return scm_from_double (rx
/ yy
);
4951 else if (SCM_COMPLEXP (y
))
4956 else if (SCM_FRACTIONP (y
))
4957 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4959 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4961 else if (SCM_COMPLEXP (x
))
4963 double rx
= SCM_COMPLEX_REAL (x
);
4964 double ix
= SCM_COMPLEX_IMAG (x
);
4965 if (SCM_I_INUMP (y
))
4967 long int yy
= SCM_I_INUM (y
);
4968 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4970 scm_num_overflow (s_divide
);
4975 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4978 else if (SCM_BIGP (y
))
4980 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4981 scm_remember_upto_here_1 (y
);
4982 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4984 else if (SCM_REALP (y
))
4986 double yy
= SCM_REAL_VALUE (y
);
4987 #ifndef ALLOW_DIVIDE_BY_ZERO
4989 scm_num_overflow (s_divide
);
4992 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4994 else if (SCM_COMPLEXP (y
))
4996 double ry
= SCM_COMPLEX_REAL (y
);
4997 double iy
= SCM_COMPLEX_IMAG (y
);
4998 if (fabs(ry
) <= fabs(iy
))
5001 double d
= iy
* (1.0 + t
* t
);
5002 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5007 double d
= ry
* (1.0 + t
* t
);
5008 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5011 else if (SCM_FRACTIONP (y
))
5013 double yy
= scm_i_fraction2double (y
);
5014 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5017 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5019 else if (SCM_FRACTIONP (x
))
5021 if (SCM_I_INUMP (y
))
5023 long int yy
= SCM_I_INUM (y
);
5024 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5026 scm_num_overflow (s_divide
);
5029 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5030 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5032 else if (SCM_BIGP (y
))
5034 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5035 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5037 else if (SCM_REALP (y
))
5039 double yy
= SCM_REAL_VALUE (y
);
5040 #ifndef ALLOW_DIVIDE_BY_ZERO
5042 scm_num_overflow (s_divide
);
5045 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5047 else if (SCM_COMPLEXP (y
))
5049 a
= scm_i_fraction2double (x
);
5052 else if (SCM_FRACTIONP (y
))
5053 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5054 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5056 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5059 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5063 scm_divide (SCM x
, SCM y
)
5065 return scm_i_divide (x
, y
, 0);
5068 static SCM
scm_divide2real (SCM x
, SCM y
)
5070 return scm_i_divide (x
, y
, 1);
5076 scm_c_truncate (double x
)
5087 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5088 half-way case (ie. when x is an integer plus 0.5) going upwards.
5089 Then half-way cases are identified and adjusted down if the
5090 round-upwards didn't give the desired even integer.
5092 "plus_half == result" identifies a half-way case. If plus_half, which is
5093 x + 0.5, is an integer then x must be an integer plus 0.5.
5095 An odd "result" value is identified with result/2 != floor(result/2).
5096 This is done with plus_half, since that value is ready for use sooner in
5097 a pipelined cpu, and we're already requiring plus_half == result.
5099 Note however that we need to be careful when x is big and already an
5100 integer. In that case "x+0.5" may round to an adjacent integer, causing
5101 us to return such a value, incorrectly. For instance if the hardware is
5102 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5103 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5104 returned. Or if the hardware is in round-upwards mode, then other bigger
5105 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5106 representable value, 2^128+2^76 (or whatever), again incorrect.
5108 These bad roundings of x+0.5 are avoided by testing at the start whether
5109 x is already an integer. If it is then clearly that's the desired result
5110 already. And if it's not then the exponent must be small enough to allow
5111 an 0.5 to be represented, and hence added without a bad rounding. */
5114 scm_c_round (double x
)
5116 double plus_half
, result
;
5121 plus_half
= x
+ 0.5;
5122 result
= floor (plus_half
);
5123 /* Adjust so that the rounding is towards even. */
5124 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5129 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5131 "Round the number @var{x} towards zero.")
5132 #define FUNC_NAME s_scm_truncate_number
5134 if (scm_is_false (scm_negative_p (x
)))
5135 return scm_floor (x
);
5137 return scm_ceiling (x
);
5141 static SCM exactly_one_half
;
5143 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5145 "Round the number @var{x} towards the nearest integer. "
5146 "When it is exactly halfway between two integers, "
5147 "round towards the even one.")
5148 #define FUNC_NAME s_scm_round_number
5150 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5152 else if (SCM_REALP (x
))
5153 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5156 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5157 single quotient+remainder division then examining to see which way
5158 the rounding should go. */
5159 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5160 SCM result
= scm_floor (plus_half
);
5161 /* Adjust so that the rounding is towards even. */
5162 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5163 && scm_is_true (scm_odd_p (result
)))
5164 return scm_difference (result
, SCM_I_MAKINUM (1));
5171 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5173 "Round the number @var{x} towards minus infinity.")
5174 #define FUNC_NAME s_scm_floor
5176 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5178 else if (SCM_REALP (x
))
5179 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5180 else if (SCM_FRACTIONP (x
))
5182 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5183 SCM_FRACTION_DENOMINATOR (x
));
5184 if (scm_is_false (scm_negative_p (x
)))
5186 /* For positive x, rounding towards zero is correct. */
5191 /* For negative x, we need to return q-1 unless x is an
5192 integer. But fractions are never integer, per our
5194 return scm_difference (q
, SCM_I_MAKINUM (1));
5198 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5202 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5204 "Round the number @var{x} towards infinity.")
5205 #define FUNC_NAME s_scm_ceiling
5207 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5209 else if (SCM_REALP (x
))
5210 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5211 else if (SCM_FRACTIONP (x
))
5213 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5214 SCM_FRACTION_DENOMINATOR (x
));
5215 if (scm_is_false (scm_positive_p (x
)))
5217 /* For negative x, rounding towards zero is correct. */
5222 /* For positive x, we need to return q+1 unless x is an
5223 integer. But fractions are never integer, per our
5225 return scm_sum (q
, SCM_I_MAKINUM (1));
5229 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5233 /* sin/cos/tan/asin/acos/atan
5234 sinh/cosh/tanh/asinh/acosh/atanh
5235 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5236 Written by Jerry D. Hedden, (C) FSF.
5237 See the file `COPYING' for terms applying to this program. */
5239 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5241 "Return @var{x} raised to the power of @var{y}.")
5242 #define FUNC_NAME s_scm_expt
5244 if (!SCM_INEXACTP (y
) && scm_is_integer (y
))
5245 return scm_integer_expt (x
, y
);
5246 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5248 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5251 return scm_exp (scm_product (scm_log (x
), y
));
5255 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5257 "Compute the sine of @var{z}.")
5258 #define FUNC_NAME s_scm_sin
5260 if (scm_is_real (z
))
5261 return scm_from_double (sin (scm_to_double (z
)));
5262 else if (SCM_COMPLEXP (z
))
5264 x
= SCM_COMPLEX_REAL (z
);
5265 y
= SCM_COMPLEX_IMAG (z
);
5266 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5267 cos (x
) * sinh (y
));
5270 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5274 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5276 "Compute the cosine of @var{z}.")
5277 #define FUNC_NAME s_scm_cos
5279 if (scm_is_real (z
))
5280 return scm_from_double (cos (scm_to_double (z
)));
5281 else if (SCM_COMPLEXP (z
))
5283 x
= SCM_COMPLEX_REAL (z
);
5284 y
= SCM_COMPLEX_IMAG (z
);
5285 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5286 -sin (x
) * sinh (y
));
5289 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5293 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5295 "Compute the tangent of @var{z}.")
5296 #define FUNC_NAME s_scm_tan
5298 if (scm_is_real (z
))
5299 return scm_from_double (tan (scm_to_double (z
)));
5300 else if (SCM_COMPLEXP (z
))
5302 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5303 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5304 w
= cos (x
) + cosh (y
);
5305 #ifndef ALLOW_DIVIDE_BY_ZERO
5307 scm_num_overflow (s_scm_tan
);
5309 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5312 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5316 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5318 "Compute the hyperbolic sine of @var{z}.")
5319 #define FUNC_NAME s_scm_sinh
5321 if (scm_is_real (z
))
5322 return scm_from_double (sinh (scm_to_double (z
)));
5323 else if (SCM_COMPLEXP (z
))
5325 x
= SCM_COMPLEX_REAL (z
);
5326 y
= SCM_COMPLEX_IMAG (z
);
5327 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5328 cosh (x
) * sin (y
));
5331 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5335 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5337 "Compute the hyperbolic cosine of @var{z}.")
5338 #define FUNC_NAME s_scm_cosh
5340 if (scm_is_real (z
))
5341 return scm_from_double (cosh (scm_to_double (z
)));
5342 else if (SCM_COMPLEXP (z
))
5344 x
= SCM_COMPLEX_REAL (z
);
5345 y
= SCM_COMPLEX_IMAG (z
);
5346 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5347 sinh (x
) * sin (y
));
5350 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5354 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5356 "Compute the hyperbolic tangent of @var{z}.")
5357 #define FUNC_NAME s_scm_tanh
5359 if (scm_is_real (z
))
5360 return scm_from_double (tanh (scm_to_double (z
)));
5361 else if (SCM_COMPLEXP (z
))
5363 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5364 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5365 w
= cosh (x
) + cos (y
);
5366 #ifndef ALLOW_DIVIDE_BY_ZERO
5368 scm_num_overflow (s_scm_tanh
);
5370 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5373 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5377 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5379 "Compute the arc sine of @var{z}.")
5380 #define FUNC_NAME s_scm_asin
5382 if (scm_is_real (z
))
5384 double w
= scm_to_double (z
);
5385 if (w
>= -1.0 && w
<= 1.0)
5386 return scm_from_double (asin (w
));
5388 return scm_product (scm_c_make_rectangular (0, -1),
5389 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5391 else if (SCM_COMPLEXP (z
))
5393 x
= SCM_COMPLEX_REAL (z
);
5394 y
= SCM_COMPLEX_IMAG (z
);
5395 return scm_product (scm_c_make_rectangular (0, -1),
5396 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5399 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5403 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5405 "Compute the arc cosine of @var{z}.")
5406 #define FUNC_NAME s_scm_acos
5408 if (scm_is_real (z
))
5410 double w
= scm_to_double (z
);
5411 if (w
>= -1.0 && w
<= 1.0)
5412 return scm_from_double (acos (w
));
5414 return scm_sum (scm_from_double (acos (0.0)),
5415 scm_product (scm_c_make_rectangular (0, 1),
5416 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5418 else if (SCM_COMPLEXP (z
))
5420 x
= SCM_COMPLEX_REAL (z
);
5421 y
= SCM_COMPLEX_IMAG (z
);
5422 return scm_sum (scm_from_double (acos (0.0)),
5423 scm_product (scm_c_make_rectangular (0, 1),
5424 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5427 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5431 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5433 "With one argument, compute the arc tangent of @var{z}.\n"
5434 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5435 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5436 #define FUNC_NAME s_scm_atan
5440 if (scm_is_real (z
))
5441 return scm_from_double (atan (scm_to_double (z
)));
5442 else if (SCM_COMPLEXP (z
))
5445 v
= SCM_COMPLEX_REAL (z
);
5446 w
= SCM_COMPLEX_IMAG (z
);
5447 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5448 scm_c_make_rectangular (v
, w
+ 1.0))),
5449 scm_c_make_rectangular (0, 2));
5452 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5454 else if (scm_is_real (z
))
5456 if (scm_is_real (y
))
5457 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5459 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5462 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5466 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5468 "Compute the inverse hyperbolic sine of @var{z}.")
5469 #define FUNC_NAME s_scm_sys_asinh
5471 if (scm_is_real (z
))
5472 return scm_from_double (asinh (scm_to_double (z
)));
5473 else if (scm_is_number (z
))
5474 return scm_log (scm_sum (z
,
5475 scm_sqrt (scm_sum (scm_product (z
, z
),
5476 SCM_I_MAKINUM (1)))));
5478 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5482 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5484 "Compute the inverse hyperbolic cosine of @var{z}.")
5485 #define FUNC_NAME s_scm_sys_acosh
5487 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5488 return scm_from_double (acosh (scm_to_double (z
)));
5489 else if (scm_is_number (z
))
5490 return scm_log (scm_sum (z
,
5491 scm_sqrt (scm_difference (scm_product (z
, z
),
5492 SCM_I_MAKINUM (1)))));
5494 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5498 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5500 "Compute the inverse hyperbolic tangent of @var{z}.")
5501 #define FUNC_NAME s_scm_sys_atanh
5503 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5504 return scm_from_double (atanh (scm_to_double (z
)));
5505 else if (scm_is_number (z
))
5506 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5507 scm_difference (SCM_I_MAKINUM (1), z
))),
5510 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5515 scm_c_make_rectangular (double re
, double im
)
5518 return scm_from_double (re
);
5522 SCM_NEWSMOB (z
, scm_tc16_complex
,
5523 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5525 SCM_COMPLEX_REAL (z
) = re
;
5526 SCM_COMPLEX_IMAG (z
) = im
;
5531 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5532 (SCM real_part
, SCM imaginary_part
),
5533 "Return a complex number constructed of the given @var{real-part} "
5534 "and @var{imaginary-part} parts.")
5535 #define FUNC_NAME s_scm_make_rectangular
5537 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5538 SCM_ARG1
, FUNC_NAME
, "real");
5539 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5540 SCM_ARG2
, FUNC_NAME
, "real");
5541 return scm_c_make_rectangular (scm_to_double (real_part
),
5542 scm_to_double (imaginary_part
));
5547 scm_c_make_polar (double mag
, double ang
)
5551 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5552 use it on Glibc-based systems that have it (it's a GNU extension). See
5553 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5555 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5556 sincos (ang
, &s
, &c
);
5561 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5564 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5566 "Return the complex number @var{x} * e^(i * @var{y}).")
5567 #define FUNC_NAME s_scm_make_polar
5569 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5570 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5571 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5576 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5577 /* "Return the real part of the number @var{z}."
5580 scm_real_part (SCM z
)
5582 if (SCM_I_INUMP (z
))
5584 else if (SCM_BIGP (z
))
5586 else if (SCM_REALP (z
))
5588 else if (SCM_COMPLEXP (z
))
5589 return scm_from_double (SCM_COMPLEX_REAL (z
));
5590 else if (SCM_FRACTIONP (z
))
5593 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5597 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5598 /* "Return the imaginary part of the number @var{z}."
5601 scm_imag_part (SCM z
)
5603 if (SCM_I_INUMP (z
))
5605 else if (SCM_BIGP (z
))
5607 else if (SCM_REALP (z
))
5609 else if (SCM_COMPLEXP (z
))
5610 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5611 else if (SCM_FRACTIONP (z
))
5614 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5617 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5618 /* "Return the numerator of the number @var{z}."
5621 scm_numerator (SCM z
)
5623 if (SCM_I_INUMP (z
))
5625 else if (SCM_BIGP (z
))
5627 else if (SCM_FRACTIONP (z
))
5628 return SCM_FRACTION_NUMERATOR (z
);
5629 else if (SCM_REALP (z
))
5630 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5632 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5636 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5637 /* "Return the denominator of the number @var{z}."
5640 scm_denominator (SCM z
)
5642 if (SCM_I_INUMP (z
))
5643 return SCM_I_MAKINUM (1);
5644 else if (SCM_BIGP (z
))
5645 return SCM_I_MAKINUM (1);
5646 else if (SCM_FRACTIONP (z
))
5647 return SCM_FRACTION_DENOMINATOR (z
);
5648 else if (SCM_REALP (z
))
5649 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5651 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5654 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5655 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5656 * "@code{abs} for real arguments, but also allows complex numbers."
5659 scm_magnitude (SCM z
)
5661 if (SCM_I_INUMP (z
))
5663 long int zz
= SCM_I_INUM (z
);
5666 else if (SCM_POSFIXABLE (-zz
))
5667 return SCM_I_MAKINUM (-zz
);
5669 return scm_i_long2big (-zz
);
5671 else if (SCM_BIGP (z
))
5673 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5674 scm_remember_upto_here_1 (z
);
5676 return scm_i_clonebig (z
, 0);
5680 else if (SCM_REALP (z
))
5681 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5682 else if (SCM_COMPLEXP (z
))
5683 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5684 else if (SCM_FRACTIONP (z
))
5686 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5688 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5689 SCM_FRACTION_DENOMINATOR (z
));
5692 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5696 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5697 /* "Return the angle of the complex number @var{z}."
5702 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5703 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5704 But if atan2 follows the floating point rounding mode, then the value
5705 is not a constant. Maybe it'd be close enough though. */
5706 if (SCM_I_INUMP (z
))
5708 if (SCM_I_INUM (z
) >= 0)
5711 return scm_from_double (atan2 (0.0, -1.0));
5713 else if (SCM_BIGP (z
))
5715 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5716 scm_remember_upto_here_1 (z
);
5718 return scm_from_double (atan2 (0.0, -1.0));
5722 else if (SCM_REALP (z
))
5724 if (SCM_REAL_VALUE (z
) >= 0)
5727 return scm_from_double (atan2 (0.0, -1.0));
5729 else if (SCM_COMPLEXP (z
))
5730 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5731 else if (SCM_FRACTIONP (z
))
5733 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5735 else return scm_from_double (atan2 (0.0, -1.0));
5738 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5742 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5743 /* Convert the number @var{x} to its inexact representation.\n"
5746 scm_exact_to_inexact (SCM z
)
5748 if (SCM_I_INUMP (z
))
5749 return scm_from_double ((double) SCM_I_INUM (z
));
5750 else if (SCM_BIGP (z
))
5751 return scm_from_double (scm_i_big2dbl (z
));
5752 else if (SCM_FRACTIONP (z
))
5753 return scm_from_double (scm_i_fraction2double (z
));
5754 else if (SCM_INEXACTP (z
))
5757 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5761 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5763 "Return an exact number that is numerically closest to @var{z}.")
5764 #define FUNC_NAME s_scm_inexact_to_exact
5766 if (SCM_I_INUMP (z
))
5768 else if (SCM_BIGP (z
))
5770 else if (SCM_REALP (z
))
5772 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5773 SCM_OUT_OF_RANGE (1, z
);
5780 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5781 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5782 scm_i_mpz2num (mpq_denref (frac
)));
5784 /* When scm_i_make_ratio throws, we leak the memory allocated
5791 else if (SCM_FRACTIONP (z
))
5794 SCM_WRONG_TYPE_ARG (1, z
);
5798 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5800 "Returns the @emph{simplest} rational number differing\n"
5801 "from @var{x} by no more than @var{eps}.\n"
5803 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5804 "exact result when both its arguments are exact. Thus, you might need\n"
5805 "to use @code{inexact->exact} on the arguments.\n"
5808 "(rationalize (inexact->exact 1.2) 1/100)\n"
5811 #define FUNC_NAME s_scm_rationalize
5813 if (SCM_I_INUMP (x
))
5815 else if (SCM_BIGP (x
))
5817 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5819 /* Use continued fractions to find closest ratio. All
5820 arithmetic is done with exact numbers.
5823 SCM ex
= scm_inexact_to_exact (x
);
5824 SCM int_part
= scm_floor (ex
);
5825 SCM tt
= SCM_I_MAKINUM (1);
5826 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5827 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5831 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5834 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5835 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5837 /* We stop after a million iterations just to be absolutely sure
5838 that we don't go into an infinite loop. The process normally
5839 converges after less than a dozen iterations.
5842 eps
= scm_abs (eps
);
5843 while (++i
< 1000000)
5845 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5846 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5847 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5849 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5850 eps
))) /* abs(x-a/b) <= eps */
5852 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5853 if (scm_is_false (scm_exact_p (x
))
5854 || scm_is_false (scm_exact_p (eps
)))
5855 return scm_exact_to_inexact (res
);
5859 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5861 tt
= scm_floor (rx
); /* tt = floor (rx) */
5867 scm_num_overflow (s_scm_rationalize
);
5870 SCM_WRONG_TYPE_ARG (1, x
);
5874 /* conversion functions */
5877 scm_is_integer (SCM val
)
5879 return scm_is_true (scm_integer_p (val
));
5883 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5885 if (SCM_I_INUMP (val
))
5887 scm_t_signed_bits n
= SCM_I_INUM (val
);
5888 return n
>= min
&& n
<= max
;
5890 else if (SCM_BIGP (val
))
5892 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5894 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5896 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5898 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5899 return n
>= min
&& n
<= max
;
5909 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5910 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5913 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5914 SCM_I_BIG_MPZ (val
));
5916 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5928 return n
>= min
&& n
<= max
;
5936 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5938 if (SCM_I_INUMP (val
))
5940 scm_t_signed_bits n
= SCM_I_INUM (val
);
5941 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5943 else if (SCM_BIGP (val
))
5945 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5947 else if (max
<= ULONG_MAX
)
5949 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5951 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5952 return n
>= min
&& n
<= max
;
5962 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5965 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5966 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5969 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5970 SCM_I_BIG_MPZ (val
));
5972 return n
>= min
&& n
<= max
;
5980 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5982 scm_error (scm_out_of_range_key
,
5984 "Value out of range ~S to ~S: ~S",
5985 scm_list_3 (min
, max
, bad_val
),
5986 scm_list_1 (bad_val
));
5989 #define TYPE scm_t_intmax
5990 #define TYPE_MIN min
5991 #define TYPE_MAX max
5992 #define SIZEOF_TYPE 0
5993 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5994 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5995 #include "libguile/conv-integer.i.c"
5997 #define TYPE scm_t_uintmax
5998 #define TYPE_MIN min
5999 #define TYPE_MAX max
6000 #define SIZEOF_TYPE 0
6001 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6002 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6003 #include "libguile/conv-uinteger.i.c"
6005 #define TYPE scm_t_int8
6006 #define TYPE_MIN SCM_T_INT8_MIN
6007 #define TYPE_MAX SCM_T_INT8_MAX
6008 #define SIZEOF_TYPE 1
6009 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6010 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6011 #include "libguile/conv-integer.i.c"
6013 #define TYPE scm_t_uint8
6015 #define TYPE_MAX SCM_T_UINT8_MAX
6016 #define SIZEOF_TYPE 1
6017 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6018 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6019 #include "libguile/conv-uinteger.i.c"
6021 #define TYPE scm_t_int16
6022 #define TYPE_MIN SCM_T_INT16_MIN
6023 #define TYPE_MAX SCM_T_INT16_MAX
6024 #define SIZEOF_TYPE 2
6025 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6026 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6027 #include "libguile/conv-integer.i.c"
6029 #define TYPE scm_t_uint16
6031 #define TYPE_MAX SCM_T_UINT16_MAX
6032 #define SIZEOF_TYPE 2
6033 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6034 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6035 #include "libguile/conv-uinteger.i.c"
6037 #define TYPE scm_t_int32
6038 #define TYPE_MIN SCM_T_INT32_MIN
6039 #define TYPE_MAX SCM_T_INT32_MAX
6040 #define SIZEOF_TYPE 4
6041 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6042 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6043 #include "libguile/conv-integer.i.c"
6045 #define TYPE scm_t_uint32
6047 #define TYPE_MAX SCM_T_UINT32_MAX
6048 #define SIZEOF_TYPE 4
6049 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6050 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6051 #include "libguile/conv-uinteger.i.c"
6053 #define TYPE scm_t_wchar
6054 #define TYPE_MIN (scm_t_int32)-1
6055 #define TYPE_MAX (scm_t_int32)0x10ffff
6056 #define SIZEOF_TYPE 4
6057 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6058 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6059 #include "libguile/conv-integer.i.c"
6061 #if SCM_HAVE_T_INT64
6063 #define TYPE scm_t_int64
6064 #define TYPE_MIN SCM_T_INT64_MIN
6065 #define TYPE_MAX SCM_T_INT64_MAX
6066 #define SIZEOF_TYPE 8
6067 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6068 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6069 #include "libguile/conv-integer.i.c"
6071 #define TYPE scm_t_uint64
6073 #define TYPE_MAX SCM_T_UINT64_MAX
6074 #define SIZEOF_TYPE 8
6075 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6076 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6077 #include "libguile/conv-uinteger.i.c"
6082 scm_to_mpz (SCM val
, mpz_t rop
)
6084 if (SCM_I_INUMP (val
))
6085 mpz_set_si (rop
, SCM_I_INUM (val
));
6086 else if (SCM_BIGP (val
))
6087 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6089 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6093 scm_from_mpz (mpz_t val
)
6095 return scm_i_mpz2num (val
);
6099 scm_is_real (SCM val
)
6101 return scm_is_true (scm_real_p (val
));
6105 scm_is_rational (SCM val
)
6107 return scm_is_true (scm_rational_p (val
));
6111 scm_to_double (SCM val
)
6113 if (SCM_I_INUMP (val
))
6114 return SCM_I_INUM (val
);
6115 else if (SCM_BIGP (val
))
6116 return scm_i_big2dbl (val
);
6117 else if (SCM_FRACTIONP (val
))
6118 return scm_i_fraction2double (val
);
6119 else if (SCM_REALP (val
))
6120 return SCM_REAL_VALUE (val
);
6122 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6126 scm_from_double (double val
)
6128 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
6129 SCM_REAL_VALUE (z
) = val
;
6133 #if SCM_ENABLE_DISCOURAGED == 1
6136 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6140 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6144 scm_out_of_range (NULL
, num
);
6147 return scm_to_double (num
);
6151 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6155 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6159 scm_out_of_range (NULL
, num
);
6162 return scm_to_double (num
);
6168 scm_is_complex (SCM val
)
6170 return scm_is_true (scm_complex_p (val
));
6174 scm_c_real_part (SCM z
)
6176 if (SCM_COMPLEXP (z
))
6177 return SCM_COMPLEX_REAL (z
);
6180 /* Use the scm_real_part to get proper error checking and
6183 return scm_to_double (scm_real_part (z
));
6188 scm_c_imag_part (SCM z
)
6190 if (SCM_COMPLEXP (z
))
6191 return SCM_COMPLEX_IMAG (z
);
6194 /* Use the scm_imag_part to get proper error checking and
6195 dispatching. The result will almost always be 0.0, but not
6198 return scm_to_double (scm_imag_part (z
));
6203 scm_c_magnitude (SCM z
)
6205 return scm_to_double (scm_magnitude (z
));
6211 return scm_to_double (scm_angle (z
));
6215 scm_is_number (SCM z
)
6217 return scm_is_true (scm_number_p (z
));
6221 /* In the following functions we dispatch to the real-arg funcs like log()
6222 when we know the arg is real, instead of just handing everything to
6223 clog() for instance. This is in case clog() doesn't optimize for a
6224 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6225 well use it to go straight to the applicable C func. */
6227 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6229 "Return the natural logarithm of @var{z}.")
6230 #define FUNC_NAME s_scm_log
6232 if (SCM_COMPLEXP (z
))
6234 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6235 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6237 double re
= SCM_COMPLEX_REAL (z
);
6238 double im
= SCM_COMPLEX_IMAG (z
);
6239 return scm_c_make_rectangular (log (hypot (re
, im
)),
6245 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6246 although the value itself overflows. */
6247 double re
= scm_to_double (z
);
6248 double l
= log (fabs (re
));
6250 return scm_from_double (l
);
6252 return scm_c_make_rectangular (l
, M_PI
);
6258 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6260 "Return the base 10 logarithm of @var{z}.")
6261 #define FUNC_NAME s_scm_log10
6263 if (SCM_COMPLEXP (z
))
6265 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6266 clog() and a multiply by M_LOG10E, rather than the fallback
6267 log10+hypot+atan2.) */
6268 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6269 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6271 double re
= SCM_COMPLEX_REAL (z
);
6272 double im
= SCM_COMPLEX_IMAG (z
);
6273 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6274 M_LOG10E
* atan2 (im
, re
));
6279 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6280 although the value itself overflows. */
6281 double re
= scm_to_double (z
);
6282 double l
= log10 (fabs (re
));
6284 return scm_from_double (l
);
6286 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6292 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6294 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6295 "base of natural logarithms (2.71828@dots{}).")
6296 #define FUNC_NAME s_scm_exp
6298 if (SCM_COMPLEXP (z
))
6300 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6301 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6303 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6304 SCM_COMPLEX_IMAG (z
));
6309 /* When z is a negative bignum the conversion to double overflows,
6310 giving -infinity, but that's ok, the exp is still 0.0. */
6311 return scm_from_double (exp (scm_to_double (z
)));
6317 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6319 "Return the square root of @var{z}. Of the two possible roots\n"
6320 "(positive and negative), the one with the a positive real part\n"
6321 "is returned, or if that's zero then a positive imaginary part.\n"
6325 "(sqrt 9.0) @result{} 3.0\n"
6326 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6327 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6328 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6330 #define FUNC_NAME s_scm_sqrt
6332 if (SCM_COMPLEXP (x
))
6334 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6335 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6337 double re
= SCM_COMPLEX_REAL (x
);
6338 double im
= SCM_COMPLEX_IMAG (x
);
6339 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6340 0.5 * atan2 (im
, re
));
6345 double xx
= scm_to_double (x
);
6347 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6349 return scm_from_double (sqrt (xx
));
6361 mpz_init_set_si (z_negative_one
, -1);
6363 /* It may be possible to tune the performance of some algorithms by using
6364 * the following constants to avoid the creation of bignums. Please, before
6365 * using these values, remember the two rules of program optimization:
6366 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6367 scm_c_define ("most-positive-fixnum",
6368 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6369 scm_c_define ("most-negative-fixnum",
6370 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6372 scm_add_feature ("complex");
6373 scm_add_feature ("inexact");
6374 scm_flo0
= scm_from_double (0.0);
6376 /* determine floating point precision */
6377 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6379 init_dblprec(&scm_dblprec
[i
-2],i
);
6380 init_fx_radix(fx_per_radix
[i
-2],i
);
6383 /* hard code precision for base 10 if the preprocessor tells us to... */
6384 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6387 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6388 SCM_I_MAKINUM (2)));
6389 #include "libguile/numbers.x"