diff options
Diffstat (limited to 'src/floatfns.c')
| -rw-r--r-- | src/floatfns.c | 230 |
1 files changed, 137 insertions, 93 deletions
diff --git a/src/floatfns.c b/src/floatfns.c index 13ecc66fbfa..2d76b97eec7 100644 --- a/src/floatfns.c +++ b/src/floatfns.c @@ -42,18 +42,12 @@ along with GNU Emacs. If not, see <https://www.gnu.org/licenses/>. */ #include <config.h> #include "lisp.h" +#include "bignum.h" #include <math.h> #include <count-leading-zeros.h> -#ifndef isfinite -# define isfinite(x) ((x) - (x) == 0) -#endif -#ifndef isnan -# define isnan(x) ((x) != (x)) -#endif - /* Check that X is a floating point number. */ static void @@ -67,7 +61,7 @@ CHECK_FLOAT (Lisp_Object x) double extract_float (Lisp_Object num) { - CHECK_NUMBER_OR_FLOAT (num); + CHECK_NUMBER (num); return XFLOATINT (num); } @@ -185,7 +179,7 @@ If X is zero, both parts (SGNFCAND and EXP) are zero. */) double f = extract_float (x); int exponent; double sgnfcand = frexp (f, &exponent); - return Fcons (make_float (sgnfcand), make_number (exponent)); + return Fcons (make_float (sgnfcand), make_fixnum (exponent)); } DEFUN ("ldexp", Fldexp, Sldexp, 2, 2, 0, @@ -193,8 +187,8 @@ DEFUN ("ldexp", Fldexp, Sldexp, 2, 2, 0, EXPONENT must be an integer. */) (Lisp_Object sgnfcand, Lisp_Object exponent) { - CHECK_NUMBER (exponent); - int e = min (max (INT_MIN, XINT (exponent)), INT_MAX); + CHECK_FIXNUM (exponent); + int e = min (max (INT_MIN, XFIXNUM (exponent)), INT_MAX); return make_float (ldexp (extract_float (sgnfcand), e)); } @@ -211,29 +205,14 @@ DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0, doc: /* Return the exponential ARG1 ** ARG2. */) (Lisp_Object arg1, Lisp_Object arg2) { - CHECK_NUMBER_OR_FLOAT (arg1); - CHECK_NUMBER_OR_FLOAT (arg2); - if (INTEGERP (arg1) /* common lisp spec */ - && INTEGERP (arg2) /* don't promote, if both are ints, and */ - && XINT (arg2) >= 0) /* we are sure the result is not fractional */ - { /* this can be improved by pre-calculating */ - EMACS_INT y; /* some binary powers of x then accumulating */ - EMACS_UINT acc, x; /* Unsigned so that overflow is well defined. */ - Lisp_Object val; - - x = XINT (arg1); - y = XINT (arg2); - acc = (y & 1 ? x : 1); - - while ((y >>= 1) != 0) - { - x *= x; - if (y & 1) - acc *= x; - } - XSETINT (val, acc); - return val; - } + CHECK_NUMBER (arg1); + CHECK_NUMBER (arg2); + + /* Common Lisp spec: don't promote if both are integers, and if the + result is not fractional. */ + if (INTEGERP (arg1) && !NILP (Fnatnump (arg2))) + return expt_integer (arg1, arg2); + return make_float (pow (XFLOATINT (arg1), XFLOATINT (arg2))); } @@ -273,14 +252,28 @@ DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0, DEFUN ("abs", Fabs, Sabs, 1, 1, 0, doc: /* Return the absolute value of ARG. */) - (register Lisp_Object arg) + (Lisp_Object arg) { - CHECK_NUMBER_OR_FLOAT (arg); + CHECK_NUMBER (arg); - if (FLOATP (arg)) - arg = make_float (fabs (XFLOAT_DATA (arg))); - else if (XINT (arg) < 0) - XSETINT (arg, - XINT (arg)); + if (FIXNUMP (arg)) + { + if (XFIXNUM (arg) < 0) + arg = make_int (-XFIXNUM (arg)); + } + else if (FLOATP (arg)) + { + if (signbit (XFLOAT_DATA (arg))) + arg = make_float (- XFLOAT_DATA (arg)); + } + else + { + if (mpz_sgn (XBIGNUM (arg)->value) < 0) + { + mpz_neg (mpz[0], XBIGNUM (arg)->value); + arg = make_integer_mpz (); + } + } return arg; } @@ -289,12 +282,9 @@ DEFUN ("float", Ffloat, Sfloat, 1, 1, 0, doc: /* Return the floating point number equal to ARG. */) (register Lisp_Object arg) { - CHECK_NUMBER_OR_FLOAT (arg); - - if (INTEGERP (arg)) - return make_float ((double) XINT (arg)); - else /* give 'em the same float back */ - return arg; + CHECK_NUMBER (arg); + /* If ARG is a float, give 'em the same float back. */ + return FLOATP (arg) ? arg : make_float (XFLOATINT (arg)); } static int @@ -311,7 +301,7 @@ This is the same as the exponent of a float. */) (Lisp_Object arg) { EMACS_INT value; - CHECK_NUMBER_OR_FLOAT (arg); + CHECK_NUMBER (arg); if (FLOATP (arg)) { @@ -328,27 +318,42 @@ This is the same as the exponent of a float. */) else value = MOST_POSITIVE_FIXNUM; } + else if (BIGNUMP (arg)) + value = mpz_sizeinbase (XBIGNUM (arg)->value, 2) - 1; else { - EMACS_INT i = eabs (XINT (arg)); + eassert (FIXNUMP (arg)); + EMACS_INT i = eabs (XFIXNUM (arg)); value = (i == 0 ? MOST_NEGATIVE_FIXNUM : EMACS_UINT_WIDTH - 1 - ecount_leading_zeros (i)); } - return make_number (value); + return make_fixnum (value); } +/* True if A is exactly representable as an integer. */ + +static bool +integer_value (Lisp_Object a) +{ + if (FLOATP (a)) + { + double d = XFLOAT_DATA (a); + return d == floor (d) && isfinite (d); + } + return true; +} /* the rounding functions */ static Lisp_Object rounding_driver (Lisp_Object arg, Lisp_Object divisor, double (*double_round) (double), - EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT), - const char *name) + void (*int_divide) (mpz_t, mpz_t const, mpz_t const), + EMACS_INT (*fixnum_divide) (EMACS_INT, EMACS_INT)) { - CHECK_NUMBER_OR_FLOAT (arg); + CHECK_NUMBER (arg); double d; if (NILP (divisor)) @@ -359,18 +364,36 @@ rounding_driver (Lisp_Object arg, Lisp_Object divisor, } else { - CHECK_NUMBER_OR_FLOAT (divisor); - if (!FLOATP (arg) && !FLOATP (divisor)) + CHECK_NUMBER (divisor); + if (integer_value (arg) && integer_value (divisor)) { - if (XINT (divisor) == 0) - xsignal0 (Qarith_error); - return make_number (int_round2 (XINT (arg), XINT (divisor))); + /* Divide as integers. Converting to double might lose + info, even for fixnums; also see the FIXME below. */ + + if (FLOATP (arg)) + arg = double_to_integer (XFLOAT_DATA (arg)); + if (FLOATP (divisor)) + divisor = double_to_integer (XFLOAT_DATA (divisor)); + + if (FIXNUMP (divisor)) + { + if (XFIXNUM (divisor) == 0) + xsignal0 (Qarith_error); + if (FIXNUMP (arg)) + return make_int (fixnum_divide (XFIXNUM (arg), + XFIXNUM (divisor))); + } + int_divide (mpz[0], + *bignum_integer (&mpz[0], arg), + *bignum_integer (&mpz[1], divisor)); + return make_integer_mpz (); } - double f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg); - double f2 = FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor); + double f1 = XFLOATINT (arg); + double f2 = XFLOATINT (divisor); if (! IEEE_FLOATING_POINT && f2 == 0) xsignal0 (Qarith_error); + /* FIXME: This division rounds, so the result is double-rounded. */ d = f1 / f2; } @@ -383,42 +406,61 @@ rounding_driver (Lisp_Object arg, Lisp_Object divisor, { EMACS_INT ir = dr; if (! FIXNUM_OVERFLOW_P (ir)) - return make_number (ir); + return make_fixnum (ir); } - xsignal2 (Qrange_error, build_string (name), arg); + return double_to_integer (dr); } static EMACS_INT -ceiling2 (EMACS_INT i1, EMACS_INT i2) +ceiling2 (EMACS_INT n, EMACS_INT d) { - return i1 / i2 + ((i1 % i2 != 0) & ((i1 < 0) == (i2 < 0))); + return n / d + ((n % d != 0) & ((n < 0) == (d < 0))); } static EMACS_INT -floor2 (EMACS_INT i1, EMACS_INT i2) +floor2 (EMACS_INT n, EMACS_INT d) { - return i1 / i2 - ((i1 % i2 != 0) & ((i1 < 0) != (i2 < 0))); + return n / d - ((n % d != 0) & ((n < 0) != (d < 0))); } static EMACS_INT -truncate2 (EMACS_INT i1, EMACS_INT i2) +truncate2 (EMACS_INT n, EMACS_INT d) { - return i1 / i2; + return n / d; } static EMACS_INT -round2 (EMACS_INT i1, EMACS_INT i2) -{ - /* The C language's division operator gives us one remainder R, but - we want the remainder R1 on the other side of 0 if R1 is closer - to 0 than R is; because we want to round to even, we also want R1 - if R and R1 are the same distance from 0 and if C's quotient is - odd. */ - EMACS_INT q = i1 / i2; - EMACS_INT r = i1 % i2; +round2 (EMACS_INT n, EMACS_INT d) +{ + /* The C language's division operator gives us the remainder R + corresponding to truncated division, but we want the remainder R1 + on the other side of 0 if R1 is closer to 0 than R is; because we + want to round to even, we also want R1 if R and R1 are the same + distance from 0 and if the truncated quotient is odd. */ + EMACS_INT q = n / d; + EMACS_INT r = n % d; + bool neg_d = d < 0; + bool neg_r = r < 0; EMACS_INT abs_r = eabs (r); - EMACS_INT abs_r1 = eabs (i2) - abs_r; - return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1); + EMACS_INT abs_r1 = eabs (d) - abs_r; + if (abs_r1 < abs_r + (q & 1)) + q += neg_d == neg_r ? 1 : -1; + return q; +} + +static void +rounddiv_q (mpz_t q, mpz_t const n, mpz_t const d) +{ + /* Mimic the source code of round2, using mpz_t instead of EMACS_INT. */ + mpz_t *r = &mpz[2], *abs_r = r, *abs_r1 = &mpz[3]; + mpz_tdiv_qr (q, *r, n, d); + bool neg_d = mpz_sgn (d) < 0; + bool neg_r = mpz_sgn (*r) < 0; + mpz_abs (*abs_r, *r); + mpz_abs (*abs_r1, d); + mpz_sub (*abs_r1, *abs_r1, *abs_r); + if (mpz_cmp (*abs_r1, *abs_r) < (mpz_odd_p (q) != 0)) + (neg_d == neg_r ? mpz_add_ui : mpz_sub_ui) (q, q, 1); } /* The code uses emacs_rint, so that it works to undefine HAVE_RINT @@ -435,11 +477,9 @@ emacs_rint (double d) } #endif -#ifdef HAVE_TRUNC -#define emacs_trunc trunc -#else -static double -emacs_trunc (double d) +#ifndef HAVE_TRUNC +double +trunc (double d) { return (d < 0 ? ceil : floor) (d); } @@ -451,7 +491,7 @@ This rounds the value towards +inf. With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR. */) (Lisp_Object arg, Lisp_Object divisor) { - return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling"); + return rounding_driver (arg, divisor, ceil, mpz_cdiv_q, ceiling2); } DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0, @@ -460,7 +500,7 @@ This rounds the value towards -inf. With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR. */) (Lisp_Object arg, Lisp_Object divisor) { - return rounding_driver (arg, divisor, floor, floor2, "floor"); + return rounding_driver (arg, divisor, floor, mpz_fdiv_q, floor2); } DEFUN ("round", Fround, Sround, 1, 2, 0, @@ -473,7 +513,14 @@ your machine. For example, (round 2.5) can return 3 on some systems, but 2 on others. */) (Lisp_Object arg, Lisp_Object divisor) { - return rounding_driver (arg, divisor, emacs_rint, round2, "round"); + return rounding_driver (arg, divisor, emacs_rint, rounddiv_q, round2); +} + +/* Since rounding_driver truncates anyway, no need to call 'trunc'. */ +static double +identity (double x) +{ + return x; } DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0, @@ -482,18 +529,15 @@ Rounds ARG toward zero. With optional DIVISOR, truncate ARG/DIVISOR. */) (Lisp_Object arg, Lisp_Object divisor) { - return rounding_driver (arg, divisor, emacs_trunc, truncate2, - "truncate"); + return rounding_driver (arg, divisor, identity, mpz_tdiv_q, truncate2); } Lisp_Object fmod_float (Lisp_Object x, Lisp_Object y) { - double f1, f2; - - f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x); - f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y); + double f1 = XFLOATINT (x); + double f2 = XFLOATINT (y); f1 = fmod (f1, f2); @@ -543,7 +587,7 @@ DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0, { CHECK_FLOAT (arg); double d = XFLOAT_DATA (arg); - d = emacs_trunc (d); + d = trunc (d); return make_float (d); } |