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