src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2
   3 Copyright (C) 1988, 1993-1994, 1999, 2001-2018 Free Software Foundation,
   4 Inc.
   5
   6 Author: Wolfgang Rupprecht (according to ack.texi)
   7
   8 This file is part of GNU Emacs.
   9
  10 GNU Emacs is free software: you can redistribute it and/or modify
  11 it under the terms of the GNU General Public License as published by
  12 the Free Software Foundation, either version 3 of the License, or (at
  13 your option) any later version.
  14
  15 GNU Emacs is distributed in the hope that it will be useful,
  16 but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18 GNU General Public License for more details.
  19
  20 You should have received a copy of the GNU General Public License
  21 along with GNU Emacs.  If not, see <https://www.gnu.org/licenses/>.  */
  22
  23
  24 /* C89 requires only the following math.h functions, and Emacs omits
  25    the starred functions since we haven't found a use for them:
  26    acos, asin, atan, atan2, ceil, cos, *cosh, exp, fabs, floor, fmod,
  27    frexp, ldexp, log, log10 [via (log X 10)], *modf, pow, sin, *sinh,
  28    sqrt, tan, *tanh.
  29
  30    C99 and C11 require the following math.h functions in addition to
  31    the C89 functions.  Of these, Emacs currently exports only the
  32    starred ones to Lisp, since we haven't found a use for the others:
  33    acosh, atanh, cbrt, *copysign, erf, erfc, exp2, expm1, fdim, fma,
  34    fmax, fmin, fpclassify, hypot, ilogb, isfinite, isgreater,
  35    isgreaterequal, isinf, isless, islessequal, islessgreater, *isnan,
  36    isnormal, isunordered, lgamma, log1p, *log2 [via (log X 2)], *logb
  37    (approximately), lrint/llrint, lround/llround, nan, nearbyint,
  38    nextafter, nexttoward, remainder, remquo, *rint, round, scalbln,
  39    scalbn, signbit, tgamma, *trunc.
  40  */
  41
  42 #include <config.h>
  43
  44 #include "lisp.h"
  45 #include "bignum.h"
  46
  47 #include <math.h>
  48
  49 #include <count-leading-zeros.h>
  50
  51 /* Check that X is a floating point number.  */
  52
  53 static void
  54 CHECK_FLOAT (Lisp_Object x)
  55 {
  56   CHECK_TYPE (FLOATP (x), Qfloatp, x);
  57 }
  58
  59 /* Extract a Lisp number as a `double', or signal an error.  */
  60
  61 double
  62 extract_float (Lisp_Object num)
  63 {
  64   CHECK_NUMBER (num);
  65   return XFLOATINT (num);
  66 }
  67 \f
  68 /* Trig functions.  */
  69
  70 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  71        doc: /* Return the inverse cosine of ARG.  */)
  72   (Lisp_Object arg)
  73 {
  74   double d = extract_float (arg);
  75   d = acos (d);
  76   return make_float (d);
  77 }
  78
  79 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  80        doc: /* Return the inverse sine of ARG.  */)
  81   (Lisp_Object arg)
  82 {
  83   double d = extract_float (arg);
  84   d = asin (d);
  85   return make_float (d);
  86 }
  87
  88 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
  89        doc: /* Return the inverse tangent of the arguments.
  90 If only one argument Y is given, return the inverse tangent of Y.
  91 If two arguments Y and X are given, return the inverse tangent of Y
  92 divided by X, i.e. the angle in radians between the vector (X, Y)
  93 and the x-axis.  */)
  94   (Lisp_Object y, Lisp_Object x)
  95 {
  96   double d = extract_float (y);
  97
  98   if (NILP (x))
  99     d = atan (d);
 100   else
 101     {
 102       double d2 = extract_float (x);
 103       d = atan2 (d, d2);
 104     }
 105   return make_float (d);
 106 }
 107
 108 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 109        doc: /* Return the cosine of ARG.  */)
 110   (Lisp_Object arg)
 111 {
 112   double d = extract_float (arg);
 113   d = cos (d);
 114   return make_float (d);
 115 }
 116
 117 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 118        doc: /* Return the sine of ARG.  */)
 119   (Lisp_Object arg)
 120 {
 121   double d = extract_float (arg);
 122   d = sin (d);
 123   return make_float (d);
 124 }
 125
 126 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 127        doc: /* Return the tangent of ARG.  */)
 128   (Lisp_Object arg)
 129 {
 130   double d = extract_float (arg);
 131   d = tan (d);
 132   return make_float (d);
 133 }
 134
 135 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
 136        doc: /* Return non nil if argument X is a NaN.  */)
 137   (Lisp_Object x)
 138 {
 139   CHECK_FLOAT (x);
 140   return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
 141 }
 142
 143 /* Although the substitute does not work on NaNs, it is good enough
 144    for platforms lacking the signbit macro.  */
 145 #ifndef signbit
 146 # define signbit(x) ((x) < 0 || (IEEE_FLOATING_POINT && !(x) && 1 / (x) < 0))
 147 #endif
 148
 149 DEFUN ("copysign", Fcopysign, Scopysign, 2, 2, 0,
 150        doc: /* Copy sign of X2 to value of X1, and return the result.
 151 Cause an error if X1 or X2 is not a float.  */)
 152   (Lisp_Object x1, Lisp_Object x2)
 153 {
 154   double f1, f2;
 155
 156   CHECK_FLOAT (x1);
 157   CHECK_FLOAT (x2);
 158
 159   f1 = XFLOAT_DATA (x1);
 160   f2 = XFLOAT_DATA (x2);
 161
 162   /* Use signbit instead of copysign, to avoid calling make_float when
 163      the result is X1.  */
 164   return signbit (f1) != signbit (f2) ? make_float (-f1) : x1;
 165 }
 166
 167 DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
 168        doc: /* Get significand and exponent of a floating point number.
 169 Breaks the floating point number X into its binary significand SGNFCAND
 170 \(a floating point value between 0.5 (included) and 1.0 (excluded))
 171 and an integral exponent EXP for 2, such that:
 172
 173   X = SGNFCAND * 2^EXP
 174
 175 The function returns the cons cell (SGNFCAND . EXP).
 176 If X is zero, both parts (SGNFCAND and EXP) are zero.  */)
 177   (Lisp_Object x)
 178 {
 179   double f = extract_float (x);
 180   int exponent;
 181   double sgnfcand = frexp (f, &exponent);
 182   return Fcons (make_float (sgnfcand), make_fixnum (exponent));
 183 }
 184
 185 DEFUN ("ldexp", Fldexp, Sldexp, 2, 2, 0,
 186        doc: /* Return SGNFCAND * 2**EXPONENT, as a floating point number.
 187 EXPONENT must be an integer.   */)
 188   (Lisp_Object sgnfcand, Lisp_Object exponent)
 189 {
 190   CHECK_FIXNUM (exponent);
 191   int e = min (max (INT_MIN, XFIXNUM (exponent)), INT_MAX);
 192   return make_float (ldexp (extract_float (sgnfcand), e));
 193 }
 194 \f
 195 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 196        doc: /* Return the exponential base e of ARG.  */)
 197   (Lisp_Object arg)
 198 {
 199   double d = extract_float (arg);
 200   d = exp (d);
 201   return make_float (d);
 202 }
 203
 204 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 205        doc: /* Return the exponential ARG1 ** ARG2.  */)
 206   (Lisp_Object arg1, Lisp_Object arg2)
 207 {
 208   CHECK_NUMBER (arg1);
 209   CHECK_NUMBER (arg2);
 210
 211   /* Common Lisp spec: don't promote if both are integers, and if the
 212      result is not fractional.  */
 213   if (INTEGERP (arg1) && !NILP (Fnatnump (arg2)))
 214     return expt_integer (arg1, arg2);
 215
 216   return make_float (pow (XFLOATINT (arg1), XFLOATINT (arg2)));
 217 }
 218
 219 DEFUN ("log", Flog, Slog, 1, 2, 0,
 220        doc: /* Return the natural logarithm of ARG.
 221 If the optional argument BASE is given, return log ARG using that base.  */)
 222   (Lisp_Object arg, Lisp_Object base)
 223 {
 224   double d = extract_float (arg);
 225
 226   if (NILP (base))
 227     d = log (d);
 228   else
 229     {
 230       double b = extract_float (base);
 231
 232       if (b == 10.0)
 233         d = log10 (d);
 234 #if HAVE_LOG2
 235       else if (b == 2.0)
 236         d = log2 (d);
 237 #endif
 238       else
 239         d = log (d) / log (b);
 240     }
 241   return make_float (d);
 242 }
 243
 244 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 245        doc: /* Return the square root of ARG.  */)
 246   (Lisp_Object arg)
 247 {
 248   double d = extract_float (arg);
 249   d = sqrt (d);
 250   return make_float (d);
 251 }
 252 \f
 253 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 254        doc: /* Return the absolute value of ARG.  */)
 255   (Lisp_Object arg)
 256 {
 257   CHECK_NUMBER (arg);
 258
 259   if (FIXNUMP (arg))
 260     {
 261       if (XFIXNUM (arg) < 0)
 262         arg = make_int (-XFIXNUM (arg));
 263     }
 264   else if (FLOATP (arg))
 265     {
 266       if (signbit (XFLOAT_DATA (arg)))
 267         arg = make_float (- XFLOAT_DATA (arg));
 268     }
 269   else
 270     {
 271       if (mpz_sgn (XBIGNUM (arg)->value) < 0)
 272         {
 273           mpz_neg (mpz[0], XBIGNUM (arg)->value);
 274           arg = make_integer_mpz ();
 275         }
 276     }
 277
 278   return arg;
 279 }
 280
 281 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 282        doc: /* Return the floating point number equal to ARG.  */)
 283   (register Lisp_Object arg)
 284 {
 285   CHECK_NUMBER (arg);
 286   /* If ARG is a float, give 'em the same float back.  */
 287   return FLOATP (arg) ? arg : make_float (XFLOATINT (arg));
 288 }
 289
 290 static int
 291 ecount_leading_zeros (EMACS_UINT x)
 292 {
 293   return (EMACS_UINT_WIDTH == UINT_WIDTH ? count_leading_zeros (x)
 294           : EMACS_UINT_WIDTH == ULONG_WIDTH ? count_leading_zeros_l (x)
 295           : count_leading_zeros_ll (x));
 296 }
 297
 298 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 299        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 300 This is the same as the exponent of a float.  */)
 301   (Lisp_Object arg)
 302 {
 303   EMACS_INT value;
 304   CHECK_NUMBER (arg);
 305
 306   if (FLOATP (arg))
 307     {
 308       double f = XFLOAT_DATA (arg);
 309
 310       if (f == 0)
 311         value = MOST_NEGATIVE_FIXNUM;
 312       else if (isfinite (f))
 313         {
 314           int ivalue;
 315           frexp (f, &ivalue);
 316           value = ivalue - 1;
 317         }
 318       else
 319         value = MOST_POSITIVE_FIXNUM;
 320     }
 321   else if (BIGNUMP (arg))
 322     value = mpz_sizeinbase (XBIGNUM (arg)->value, 2) - 1;
 323   else
 324     {
 325       eassert (FIXNUMP (arg));
 326       EMACS_INT i = eabs (XFIXNUM (arg));
 327       value = (i == 0
 328                ? MOST_NEGATIVE_FIXNUM
 329                : EMACS_UINT_WIDTH - 1 - ecount_leading_zeros (i));
 330     }
 331
 332   return make_fixnum (value);
 333 }
 334
 335 /* True if A is exactly representable as an integer.  */
 336
 337 static bool
 338 integer_value (Lisp_Object a)
 339 {
 340   if (FLOATP (a))
 341     {
 342       double d = XFLOAT_DATA (a);
 343       return d == floor (d) && isfinite (d);
 344     }
 345   return true;
 346 }
 347
 348 /* the rounding functions  */
 349
 350 static Lisp_Object
 351 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 352                  double (*double_round) (double),
 353                  void (*int_divide) (mpz_t, mpz_t const, mpz_t const),
 354                  EMACS_INT (*fixnum_divide) (EMACS_INT, EMACS_INT))
 355 {
 356   CHECK_NUMBER (arg);
 357
 358   double d;
 359   if (NILP (divisor))
 360     {
 361       if (! FLOATP (arg))
 362         return arg;
 363       d = XFLOAT_DATA (arg);
 364     }
 365   else
 366     {
 367       CHECK_NUMBER (divisor);
 368       if (integer_value (arg) && integer_value (divisor))
 369         {
 370           /* Divide as integers.  Converting to double might lose
 371              info, even for fixnums; also see the FIXME below.  */
 372
 373           if (FLOATP (arg))
 374             arg = double_to_integer (XFLOAT_DATA (arg));
 375           if (FLOATP (divisor))
 376             divisor = double_to_integer (XFLOAT_DATA (divisor));
 377
 378           if (FIXNUMP (divisor))
 379             {
 380               if (XFIXNUM (divisor) == 0)
 381                 xsignal0 (Qarith_error);
 382               if (FIXNUMP (arg))
 383                 return make_int (fixnum_divide (XFIXNUM (arg),
 384                                                 XFIXNUM (divisor)));
 385             }
 386           int_divide (mpz[0],
 387                       *bignum_integer (&mpz[0], arg),
 388                       *bignum_integer (&mpz[1], divisor));
 389           return make_integer_mpz ();
 390         }
 391
 392       double f1 = XFLOATINT (arg);
 393       double f2 = XFLOATINT (divisor);
 394       if (! IEEE_FLOATING_POINT && f2 == 0)
 395         xsignal0 (Qarith_error);
 396       /* FIXME: This division rounds, so the result is double-rounded.  */
 397       d = f1 / f2;
 398     }
 399
 400   /* Round, coarsely test for fixnum overflow before converting to
 401      EMACS_INT (to avoid undefined C behavior), and then exactly test
 402      for overflow after converting (as FIXNUM_OVERFLOW_P is inaccurate
 403      on floats).  */
 404   double dr = double_round (d);
 405   if (fabs (dr) < 2 * (MOST_POSITIVE_FIXNUM + 1))
 406     {
 407       EMACS_INT ir = dr;
 408       if (! FIXNUM_OVERFLOW_P (ir))
 409         return make_fixnum (ir);
 410     }
 411   return double_to_integer (dr);
 412 }
 413
 414 static EMACS_INT
 415 ceiling2 (EMACS_INT n, EMACS_INT d)
 416 {
 417   return n / d + ((n % d != 0) & ((n < 0) == (d < 0)));
 418 }
 419
 420 static EMACS_INT
 421 floor2 (EMACS_INT n, EMACS_INT d)
 422 {
 423   return n / d - ((n % d != 0) & ((n < 0) != (d < 0)));
 424 }
 425
 426 static EMACS_INT
 427 truncate2 (EMACS_INT n, EMACS_INT d)
 428 {
 429   return n / d;
 430 }
 431
 432 static EMACS_INT
 433 round2 (EMACS_INT n, EMACS_INT d)
 434 {
 435   /* The C language's division operator gives us the remainder R
 436      corresponding to truncated division, but we want the remainder R1
 437      on the other side of 0 if R1 is closer to 0 than R is; because we
 438      want to round to even, we also want R1 if R and R1 are the same
 439      distance from 0 and if the truncated quotient is odd.  */
 440   EMACS_INT q = n / d;
 441   EMACS_INT r = n % d;
 442   bool neg_d = d < 0;
 443   bool neg_r = r < 0;
 444   EMACS_INT abs_r = eabs (r);
 445   EMACS_INT abs_r1 = eabs (d) - abs_r;
 446   if (abs_r1 < abs_r + (q & 1))
 447     q += neg_d == neg_r ? 1 : -1;
 448   return q;
 449 }
 450
 451 static void
 452 rounddiv_q (mpz_t q, mpz_t const n, mpz_t const d)
 453 {
 454   /* Mimic the source code of round2, using mpz_t instead of EMACS_INT.  */
 455   mpz_t *r = &mpz[2], *abs_r = r, *abs_r1 = &mpz[3];
 456   mpz_tdiv_qr (q, *r, n, d);
 457   bool neg_d = mpz_sgn (d) < 0;
 458   bool neg_r = mpz_sgn (*r) < 0;
 459   mpz_abs (*abs_r, *r);
 460   mpz_abs (*abs_r1, d);
 461   mpz_sub (*abs_r1, *abs_r1, *abs_r);
 462   if (mpz_cmp (*abs_r1, *abs_r) < (mpz_odd_p (q) != 0))
 463     (neg_d == neg_r ? mpz_add_ui : mpz_sub_ui) (q, q, 1);
 464 }
 465
 466 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 467    if `rint' exists but does not work right.  */
 468 #ifdef HAVE_RINT
 469 #define emacs_rint rint
 470 #else
 471 static double
 472 emacs_rint (double d)
 473 {
 474   double d1 = d + 0.5;
 475   double r = floor (d1);
 476   return r - (r == d1 && fmod (r, 2) != 0);
 477 }
 478 #endif
 479
 480 #ifndef HAVE_TRUNC
 481 double
 482 trunc (double d)
 483 {
 484   return (d < 0 ? ceil : floor) (d);
 485 }
 486 #endif
 487
 488 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 489        doc: /* Return the smallest integer no less than ARG.
 490 This rounds the value towards +inf.
 491 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 492   (Lisp_Object arg, Lisp_Object divisor)
 493 {
 494   return rounding_driver (arg, divisor, ceil, mpz_cdiv_q, ceiling2);
 495 }
 496
 497 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 498        doc: /* Return the largest integer no greater than ARG.
 499 This rounds the value towards -inf.
 500 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 501   (Lisp_Object arg, Lisp_Object divisor)
 502 {
 503   return rounding_driver (arg, divisor, floor, mpz_fdiv_q, floor2);
 504 }
 505
 506 DEFUN ("round", Fround, Sround, 1, 2, 0,
 507        doc: /* Return the nearest integer to ARG.
 508 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 509
 510 Rounding a value equidistant between two integers may choose the
 511 integer closer to zero, or it may prefer an even integer, depending on
 512 your machine.  For example, (round 2.5) can return 3 on some
 513 systems, but 2 on others.  */)
 514   (Lisp_Object arg, Lisp_Object divisor)
 515 {
 516   return rounding_driver (arg, divisor, emacs_rint, rounddiv_q, round2);
 517 }
 518
 519 /* Since rounding_driver truncates anyway, no need to call 'trunc'.  */
 520 static double
 521 identity (double x)
 522 {
 523   return x;
 524 }
 525
 526 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 527        doc: /* Truncate a floating point number to an int.
 528 Rounds ARG toward zero.
 529 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 530   (Lisp_Object arg, Lisp_Object divisor)
 531 {
 532   return rounding_driver (arg, divisor, identity, mpz_tdiv_q, truncate2);
 533 }
 534
 535
 536 Lisp_Object
 537 fmod_float (Lisp_Object x, Lisp_Object y)
 538 {
 539   double f1 = XFLOATINT (x);
 540   double f2 = XFLOATINT (y);
 541
 542   f1 = fmod (f1, f2);
 543
 544   /* If the "remainder" comes out with the wrong sign, fix it.  */
 545   if (f2 < 0 ? f1 > 0 : f1 < 0)
 546     f1 += f2;
 547
 548   return make_float (f1);
 549 }
 550 \f
 551 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 552        doc: /* Return the smallest integer no less than ARG, as a float.
 553 \(Round toward +inf.)  */)
 554   (Lisp_Object arg)
 555 {
 556   CHECK_FLOAT (arg);
 557   double d = XFLOAT_DATA (arg);
 558   d = ceil (d);
 559   return make_float (d);
 560 }
 561
 562 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 563        doc: /* Return the largest integer no greater than ARG, as a float.
 564 \(Round toward -inf.)  */)
 565   (Lisp_Object arg)
 566 {
 567   CHECK_FLOAT (arg);
 568   double d = XFLOAT_DATA (arg);
 569   d = floor (d);
 570   return make_float (d);
 571 }
 572
 573 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 574        doc: /* Return the nearest integer to ARG, as a float.  */)
 575   (Lisp_Object arg)
 576 {
 577   CHECK_FLOAT (arg);
 578   double d = XFLOAT_DATA (arg);
 579   d = emacs_rint (d);
 580   return make_float (d);
 581 }
 582
 583 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 584        doc: /* Truncate a floating point number to an integral float value.
 585 \(Round toward zero.)  */)
 586   (Lisp_Object arg)
 587 {
 588   CHECK_FLOAT (arg);
 589   double d = XFLOAT_DATA (arg);
 590   d = trunc (d);
 591   return make_float (d);
 592 }
 593 \f
 594 void
 595 syms_of_floatfns (void)
 596 {
 597   defsubr (&Sacos);
 598   defsubr (&Sasin);
 599   defsubr (&Satan);
 600   defsubr (&Scos);
 601   defsubr (&Ssin);
 602   defsubr (&Stan);
 603   defsubr (&Sisnan);
 604   defsubr (&Scopysign);
 605   defsubr (&Sfrexp);
 606   defsubr (&Sldexp);
 607   defsubr (&Sfceiling);
 608   defsubr (&Sffloor);
 609   defsubr (&Sfround);
 610   defsubr (&Sftruncate);
 611   defsubr (&Sexp);
 612   defsubr (&Sexpt);
 613   defsubr (&Slog);
 614   defsubr (&Ssqrt);
 615
 616   defsubr (&Sabs);
 617   defsubr (&Sfloat);
 618   defsubr (&Slogb);
 619   defsubr (&Sceiling);
 620   defsubr (&Sfloor);
 621   defsubr (&Sround);
 622   defsubr (&Struncate);
 623 }