src/floatfns.c

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