src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2
   3 Copyright (C) 1988, 1993-1994, 1999, 2001-2017 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
  46 #include <math.h>
  47
  48 #include <count-leading-zeros.h>
  49
  50 #ifndef isfinite
  51 # define isfinite(x) ((x) - (x) == 0)
  52 #endif
  53 #ifndef isnan
  54 # define isnan(x) ((x) != (x))
  55 #endif
  56
  57 /* Check that X is a floating point number.  */
  58
  59 static void
  60 CHECK_FLOAT (Lisp_Object x)
  61 {
  62   CHECK_TYPE (FLOATP (x), Qfloatp, x);
  63 }
  64
  65 /* Extract a Lisp number as a `double', or signal an error.  */
  66
  67 double
  68 extract_float (Lisp_Object num)
  69 {
  70   CHECK_NUMBER_OR_FLOAT (num);
  71   return XFLOATINT (num);
  72 }
  73 \f
  74 /* Trig functions.  */
  75
  76 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  77        doc: /* Return the inverse cosine of ARG.  */)
  78   (Lisp_Object arg)
  79 {
  80   double d = extract_float (arg);
  81   d = acos (d);
  82   return make_float (d);
  83 }
  84
  85 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  86        doc: /* Return the inverse sine of ARG.  */)
  87   (Lisp_Object arg)
  88 {
  89   double d = extract_float (arg);
  90   d = asin (d);
  91   return make_float (d);
  92 }
  93
  94 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
  95        doc: /* Return the inverse tangent of the arguments.
  96 If only one argument Y is given, return the inverse tangent of Y.
  97 If two arguments Y and X are given, return the inverse tangent of Y
  98 divided by X, i.e. the angle in radians between the vector (X, Y)
  99 and the x-axis.  */)
 100   (Lisp_Object y, Lisp_Object x)
 101 {
 102   double d = extract_float (y);
 103
 104   if (NILP (x))
 105     d = atan (d);
 106   else
 107     {
 108       double d2 = extract_float (x);
 109       d = atan2 (d, d2);
 110     }
 111   return make_float (d);
 112 }
 113
 114 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 115        doc: /* Return the cosine of ARG.  */)
 116   (Lisp_Object arg)
 117 {
 118   double d = extract_float (arg);
 119   d = cos (d);
 120   return make_float (d);
 121 }
 122
 123 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 124        doc: /* Return the sine of ARG.  */)
 125   (Lisp_Object arg)
 126 {
 127   double d = extract_float (arg);
 128   d = sin (d);
 129   return make_float (d);
 130 }
 131
 132 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 133        doc: /* Return the tangent of ARG.  */)
 134   (Lisp_Object arg)
 135 {
 136   double d = extract_float (arg);
 137   d = tan (d);
 138   return make_float (d);
 139 }
 140
 141 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
 142        doc: /* Return non nil if argument X is a NaN.  */)
 143   (Lisp_Object x)
 144 {
 145   CHECK_FLOAT (x);
 146   return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
 147 }
 148
 149 /* Although the substitute does not work on NaNs, it is good enough
 150    for platforms lacking the signbit macro.  */
 151 #ifndef signbit
 152 # define signbit(x) ((x) < 0 || (IEEE_FLOATING_POINT && !(x) && 1 / (x) < 0))
 153 #endif
 154
 155 DEFUN ("copysign", Fcopysign, Scopysign, 2, 2, 0,
 156        doc: /* Copy sign of X2 to value of X1, and return the result.
 157 Cause an error if X1 or X2 is not a float.  */)
 158   (Lisp_Object x1, Lisp_Object x2)
 159 {
 160   double f1, f2;
 161
 162   CHECK_FLOAT (x1);
 163   CHECK_FLOAT (x2);
 164
 165   f1 = XFLOAT_DATA (x1);
 166   f2 = XFLOAT_DATA (x2);
 167
 168   /* Use signbit instead of copysign, to avoid calling make_float when
 169      the result is X1.  */
 170   return signbit (f1) != signbit (f2) ? make_float (-f1) : x1;
 171 }
 172
 173 DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
 174        doc: /* Get significand and exponent of a floating point number.
 175 Breaks the floating point number X into its binary significand SGNFCAND
 176 \(a floating point value between 0.5 (included) and 1.0 (excluded))
 177 and an integral exponent EXP for 2, such that:
 178
 179   X = SGNFCAND * 2^EXP
 180
 181 The function returns the cons cell (SGNFCAND . EXP).
 182 If X is zero, both parts (SGNFCAND and EXP) are zero.  */)
 183   (Lisp_Object x)
 184 {
 185   double f = extract_float (x);
 186   int exponent;
 187   double sgnfcand = frexp (f, &exponent);
 188   return Fcons (make_float (sgnfcand), make_number (exponent));
 189 }
 190
 191 DEFUN ("ldexp", Fldexp, Sldexp, 2, 2, 0,
 192        doc: /* Return SGNFCAND * 2**EXPONENT, as a floating point number.
 193 EXPONENT must be an integer.   */)
 194   (Lisp_Object sgnfcand, Lisp_Object exponent)
 195 {
 196   CHECK_NUMBER (exponent);
 197   int e = min (max (INT_MIN, XINT (exponent)), INT_MAX);
 198   return make_float (ldexp (extract_float (sgnfcand), e));
 199 }
 200 \f
 201 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 202        doc: /* Return the exponential base e of ARG.  */)
 203   (Lisp_Object arg)
 204 {
 205   double d = extract_float (arg);
 206   d = exp (d);
 207   return make_float (d);
 208 }
 209
 210 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 211        doc: /* Return the exponential ARG1 ** ARG2.  */)
 212   (Lisp_Object arg1, Lisp_Object arg2)
 213 {
 214   CHECK_NUMBER_OR_FLOAT (arg1);
 215   CHECK_NUMBER_OR_FLOAT (arg2);
 216   if (INTEGERP (arg1)     /* common lisp spec */
 217       && INTEGERP (arg2)   /* don't promote, if both are ints, and */
 218       && XINT (arg2) >= 0) /* we are sure the result is not fractional */
 219     {                           /* this can be improved by pre-calculating */
 220       EMACS_INT y;              /* some binary powers of x then accumulating */
 221       EMACS_UINT acc, x;  /* Unsigned so that overflow is well defined.  */
 222       Lisp_Object val;
 223
 224       x = XINT (arg1);
 225       y = XINT (arg2);
 226       acc = (y & 1 ? x : 1);
 227
 228       while ((y >>= 1) != 0)
 229         {
 230           x *= x;
 231           if (y & 1)
 232             acc *= x;
 233         }
 234       XSETINT (val, acc);
 235       return val;
 236     }
 237   return make_float (pow (XFLOATINT (arg1), XFLOATINT (arg2)));
 238 }
 239
 240 DEFUN ("log", Flog, Slog, 1, 2, 0,
 241        doc: /* Return the natural logarithm of ARG.
 242 If the optional argument BASE is given, return log ARG using that base.  */)
 243   (Lisp_Object arg, Lisp_Object base)
 244 {
 245   double d = extract_float (arg);
 246
 247   if (NILP (base))
 248     d = log (d);
 249   else
 250     {
 251       double b = extract_float (base);
 252
 253       if (b == 10.0)
 254         d = log10 (d);
 255 #if HAVE_LOG2
 256       else if (b == 2.0)
 257         d = log2 (d);
 258 #endif
 259       else
 260         d = log (d) / log (b);
 261     }
 262   return make_float (d);
 263 }
 264
 265 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 266        doc: /* Return the square root of ARG.  */)
 267   (Lisp_Object arg)
 268 {
 269   double d = extract_float (arg);
 270   d = sqrt (d);
 271   return make_float (d);
 272 }
 273 \f
 274 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 275        doc: /* Return the absolute value of ARG.  */)
 276   (register Lisp_Object arg)
 277 {
 278   CHECK_NUMBER_OR_FLOAT (arg);
 279
 280   if (FLOATP (arg))
 281     arg = make_float (fabs (XFLOAT_DATA (arg)));
 282   else if (XINT (arg) < 0)
 283     XSETINT (arg, - XINT (arg));
 284
 285   return arg;
 286 }
 287
 288 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 289        doc: /* Return the floating point number equal to ARG.  */)
 290   (register Lisp_Object arg)
 291 {
 292   CHECK_NUMBER_OR_FLOAT (arg);
 293
 294   if (INTEGERP (arg))
 295     return make_float ((double) XINT (arg));
 296   else                          /* give 'em the same float back */
 297     return arg;
 298 }
 299
 300 static int
 301 ecount_leading_zeros (EMACS_UINT x)
 302 {
 303   return (EMACS_UINT_WIDTH == UINT_WIDTH ? count_leading_zeros (x)
 304           : EMACS_UINT_WIDTH == ULONG_WIDTH ? count_leading_zeros_l (x)
 305           : count_leading_zeros_ll (x));
 306 }
 307
 308 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 309        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 310 This is the same as the exponent of a float.  */)
 311   (Lisp_Object arg)
 312 {
 313   EMACS_INT value;
 314   CHECK_NUMBER_OR_FLOAT (arg);
 315
 316   if (FLOATP (arg))
 317     {
 318       double f = XFLOAT_DATA (arg);
 319
 320       if (f == 0)
 321         value = MOST_NEGATIVE_FIXNUM;
 322       else if (isfinite (f))
 323         {
 324           int ivalue;
 325           frexp (f, &ivalue);
 326           value = ivalue - 1;
 327         }
 328       else
 329         value = MOST_POSITIVE_FIXNUM;
 330     }
 331   else
 332     {
 333       EMACS_INT i = eabs (XINT (arg));
 334       value = (i == 0
 335                ? MOST_NEGATIVE_FIXNUM
 336                : EMACS_UINT_WIDTH - 1 - ecount_leading_zeros (i));
 337     }
 338
 339   return make_number (value);
 340 }
 341
 342
 343 /* the rounding functions  */
 344
 345 static Lisp_Object
 346 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 347                  double (*double_round) (double),
 348                  EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
 349                  const char *name)
 350 {
 351   CHECK_NUMBER_OR_FLOAT (arg);
 352
 353   double d;
 354   if (NILP (divisor))
 355     {
 356       if (! FLOATP (arg))
 357         return arg;
 358       d = XFLOAT_DATA (arg);
 359     }
 360   else
 361     {
 362       CHECK_NUMBER_OR_FLOAT (divisor);
 363       if (!FLOATP (arg) && !FLOATP (divisor))
 364         {
 365           if (XINT (divisor) == 0)
 366             xsignal0 (Qarith_error);
 367           return make_number (int_round2 (XINT (arg), XINT (divisor)));
 368         }
 369
 370       double f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 371       double f2 = FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor);
 372       if (! IEEE_FLOATING_POINT && f2 == 0)
 373         xsignal0 (Qarith_error);
 374       d = f1 / f2;
 375     }
 376
 377   /* Round, coarsely test for fixnum overflow before converting to
 378      EMACS_INT (to avoid undefined C behavior), and then exactly test
 379      for overflow after converting (as FIXNUM_OVERFLOW_P is inaccurate
 380      on floats).  */
 381   double dr = double_round (d);
 382   if (fabs (dr) < 2 * (MOST_POSITIVE_FIXNUM + 1))
 383     {
 384       EMACS_INT ir = dr;
 385       if (! FIXNUM_OVERFLOW_P (ir))
 386         return make_number (ir);
 387     }
 388   xsignal2 (Qrange_error, build_string (name), arg);
 389 }
 390
 391 static EMACS_INT
 392 ceiling2 (EMACS_INT i1, EMACS_INT i2)
 393 {
 394   return i1 / i2 + ((i1 % i2 != 0) & ((i1 < 0) == (i2 < 0)));
 395 }
 396
 397 static EMACS_INT
 398 floor2 (EMACS_INT i1, EMACS_INT i2)
 399 {
 400   return i1 / i2 - ((i1 % i2 != 0) & ((i1 < 0) != (i2 < 0)));
 401 }
 402
 403 static EMACS_INT
 404 truncate2 (EMACS_INT i1, EMACS_INT i2)
 405 {
 406   return i1 / i2;
 407 }
 408
 409 static EMACS_INT
 410 round2 (EMACS_INT i1, EMACS_INT i2)
 411 {
 412   /* The C language's division operator gives us one remainder R, but
 413      we want the remainder R1 on the other side of 0 if R1 is closer
 414      to 0 than R is; because we want to round to even, we also want R1
 415      if R and R1 are the same distance from 0 and if C's quotient is
 416      odd.  */
 417   EMACS_INT q = i1 / i2;
 418   EMACS_INT r = i1 % i2;
 419   EMACS_INT abs_r = eabs (r);
 420   EMACS_INT abs_r1 = eabs (i2) - abs_r;
 421   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 422 }
 423
 424 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 425    if `rint' exists but does not work right.  */
 426 #ifdef HAVE_RINT
 427 #define emacs_rint rint
 428 #else
 429 static double
 430 emacs_rint (double d)
 431 {
 432   double d1 = d + 0.5;
 433   double r = floor (d1);
 434   return r - (r == d1 && fmod (r, 2) != 0);
 435 }
 436 #endif
 437
 438 #ifdef HAVE_TRUNC
 439 #define emacs_trunc trunc
 440 #else
 441 static double
 442 emacs_trunc (double d)
 443 {
 444   return (d < 0 ? ceil : floor) (d);
 445 }
 446 #endif
 447
 448 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 449        doc: /* Return the smallest integer no less than ARG.
 450 This rounds the value towards +inf.
 451 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 452   (Lisp_Object arg, Lisp_Object divisor)
 453 {
 454   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 455 }
 456
 457 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 458        doc: /* Return the largest integer no greater than ARG.
 459 This rounds the value towards -inf.
 460 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 461   (Lisp_Object arg, Lisp_Object divisor)
 462 {
 463   return rounding_driver (arg, divisor, floor, floor2, "floor");
 464 }
 465
 466 DEFUN ("round", Fround, Sround, 1, 2, 0,
 467        doc: /* Return the nearest integer to ARG.
 468 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 469
 470 Rounding a value equidistant between two integers may choose the
 471 integer closer to zero, or it may prefer an even integer, depending on
 472 your machine.  For example, (round 2.5) can return 3 on some
 473 systems, but 2 on others.  */)
 474   (Lisp_Object arg, Lisp_Object divisor)
 475 {
 476   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 477 }
 478
 479 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 480        doc: /* Truncate a floating point number to an int.
 481 Rounds ARG toward zero.
 482 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 483   (Lisp_Object arg, Lisp_Object divisor)
 484 {
 485   return rounding_driver (arg, divisor, emacs_trunc, truncate2,
 486                           "truncate");
 487 }
 488
 489
 490 Lisp_Object
 491 fmod_float (Lisp_Object x, Lisp_Object y)
 492 {
 493   double f1, f2;
 494
 495   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 496   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 497
 498   f1 = fmod (f1, f2);
 499
 500   /* If the "remainder" comes out with the wrong sign, fix it.  */
 501   if (f2 < 0 ? f1 > 0 : f1 < 0)
 502     f1 += f2;
 503
 504   return make_float (f1);
 505 }
 506 \f
 507 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 508        doc: /* Return the smallest integer no less than ARG, as a float.
 509 \(Round toward +inf.)  */)
 510   (Lisp_Object arg)
 511 {
 512   CHECK_FLOAT (arg);
 513   double d = XFLOAT_DATA (arg);
 514   d = ceil (d);
 515   return make_float (d);
 516 }
 517
 518 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 519        doc: /* Return the largest integer no greater than ARG, as a float.
 520 \(Round toward -inf.)  */)
 521   (Lisp_Object arg)
 522 {
 523   CHECK_FLOAT (arg);
 524   double d = XFLOAT_DATA (arg);
 525   d = floor (d);
 526   return make_float (d);
 527 }
 528
 529 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 530        doc: /* Return the nearest integer to ARG, as a float.  */)
 531   (Lisp_Object arg)
 532 {
 533   CHECK_FLOAT (arg);
 534   double d = XFLOAT_DATA (arg);
 535   d = emacs_rint (d);
 536   return make_float (d);
 537 }
 538
 539 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 540        doc: /* Truncate a floating point number to an integral float value.
 541 \(Round toward zero.)  */)
 542   (Lisp_Object arg)
 543 {
 544   CHECK_FLOAT (arg);
 545   double d = XFLOAT_DATA (arg);
 546   d = emacs_trunc (d);
 547   return make_float (d);
 548 }
 549 \f
 550 void
 551 syms_of_floatfns (void)
 552 {
 553   defsubr (&Sacos);
 554   defsubr (&Sasin);
 555   defsubr (&Satan);
 556   defsubr (&Scos);
 557   defsubr (&Ssin);
 558   defsubr (&Stan);
 559   defsubr (&Sisnan);
 560   defsubr (&Scopysign);
 561   defsubr (&Sfrexp);
 562   defsubr (&Sldexp);
 563   defsubr (&Sfceiling);
 564   defsubr (&Sffloor);
 565   defsubr (&Sfround);
 566   defsubr (&Sftruncate);
 567   defsubr (&Sexp);
 568   defsubr (&Sexpt);
 569   defsubr (&Slog);
 570   defsubr (&Ssqrt);
 571
 572   defsubr (&Sabs);
 573   defsubr (&Sfloat);
 574   defsubr (&Slogb);
 575   defsubr (&Sceiling);
 576   defsubr (&Sfloor);
 577   defsubr (&Sround);
 578   defsubr (&Struncate);
 579 }