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