src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2
   3 Copyright (C) 1988, 1993-1994, 1999, 2001-2014 Free Software Foundation, Inc.
   4
   5 Author: Wolfgang Rupprecht
   6 (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
  13 (at 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
  71   if (FLOATP (num))
  72     return XFLOAT_DATA (num);
  73   return (double) XINT (num);
  74 }
  75 \f
  76 /* Trig functions.  */
  77
  78 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  79        doc: /* Return the inverse cosine of ARG.  */)
  80   (Lisp_Object arg)
  81 {
  82   double d = extract_float (arg);
  83   d = acos (d);
  84   return make_float (d);
  85 }
  86
  87 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  88        doc: /* Return the inverse sine of ARG.  */)
  89   (Lisp_Object arg)
  90 {
  91   double d = extract_float (arg);
  92   d = asin (d);
  93   return make_float (d);
  94 }
  95
  96 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
  97        doc: /* Return the inverse tangent of the arguments.
  98 If only one argument Y is given, return the inverse tangent of Y.
  99 If two arguments Y and X are given, return the inverse tangent of Y
 100 divided by X, i.e. the angle in radians between the vector (X, Y)
 101 and the x-axis.  */)
 102   (Lisp_Object y, Lisp_Object x)
 103 {
 104   double d = extract_float (y);
 105
 106   if (NILP (x))
 107     d = atan (d);
 108   else
 109     {
 110       double d2 = extract_float (x);
 111       d = atan2 (d, d2);
 112     }
 113   return make_float (d);
 114 }
 115
 116 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 117        doc: /* Return the cosine of ARG.  */)
 118   (Lisp_Object arg)
 119 {
 120   double d = extract_float (arg);
 121   d = cos (d);
 122   return make_float (d);
 123 }
 124
 125 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 126        doc: /* Return the sine of ARG.  */)
 127   (Lisp_Object arg)
 128 {
 129   double d = extract_float (arg);
 130   d = sin (d);
 131   return make_float (d);
 132 }
 133
 134 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 135        doc: /* Return the tangent of ARG.  */)
 136   (Lisp_Object arg)
 137 {
 138   double d = extract_float (arg);
 139   d = tan (d);
 140   return make_float (d);
 141 }
 142
 143 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
 144        doc: /* Return non nil if argument X is a NaN.  */)
 145   (Lisp_Object x)
 146 {
 147   CHECK_FLOAT (x);
 148   return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
 149 }
 150
 151 #ifdef HAVE_COPYSIGN
 152 DEFUN ("copysign", Fcopysign, Scopysign, 2, 2, 0,
 153        doc: /* Copy sign of X2 to value of X1, and return the result.
 154 Cause an error if X1 or X2 is not a float.  */)
 155   (Lisp_Object x1, Lisp_Object x2)
 156 {
 157   double f1, f2;
 158
 159   CHECK_FLOAT (x1);
 160   CHECK_FLOAT (x2);
 161
 162   f1 = XFLOAT_DATA (x1);
 163   f2 = XFLOAT_DATA (x2);
 164
 165   return make_float (copysign (f1, f2));
 166 }
 167 #endif
 168
 169 DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
 170        doc: /* Get significand and exponent of a floating point number.
 171 Breaks the floating point number X into its binary significand SGNFCAND
 172 \(a floating point value between 0.5 (included) and 1.0 (excluded))
 173 and an integral exponent EXP for 2, such that:
 174
 175   X = SGNFCAND * 2^EXP
 176
 177 The function returns the cons cell (SGNFCAND . EXP).
 178 If X is zero, both parts (SGNFCAND and EXP) are zero.  */)
 179   (Lisp_Object x)
 180 {
 181   double f = XFLOATINT (x);
 182   int exponent;
 183   double sgnfcand = frexp (f, &exponent);
 184   return Fcons (make_float (sgnfcand), make_number (exponent));
 185 }
 186
 187 DEFUN ("ldexp", Fldexp, Sldexp, 1, 2, 0,
 188        doc: /* Construct number X from significand SGNFCAND and exponent EXP.
 189 Returns the floating point value resulting from multiplying SGNFCAND
 190 (the significand) by 2 raised to the power of EXP (the exponent).   */)
 191   (Lisp_Object sgnfcand, Lisp_Object exponent)
 192 {
 193   CHECK_NUMBER (exponent);
 194   return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exponent)));
 195 }
 196 \f
 197 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 198        doc: /* Return the exponential base e of ARG.  */)
 199   (Lisp_Object arg)
 200 {
 201   double d = extract_float (arg);
 202   d = exp (d);
 203   return make_float (d);
 204 }
 205
 206 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 207        doc: /* Return the exponential ARG1 ** ARG2.  */)
 208   (Lisp_Object arg1, Lisp_Object arg2)
 209 {
 210   double f1, f2, f3;
 211
 212   CHECK_NUMBER_OR_FLOAT (arg1);
 213   CHECK_NUMBER_OR_FLOAT (arg2);
 214   if (INTEGERP (arg1)     /* common lisp spec */
 215       && INTEGERP (arg2)   /* don't promote, if both are ints, and */
 216       && XINT (arg2) >= 0) /* we are sure the result is not fractional */
 217     {                           /* this can be improved by pre-calculating */
 218       EMACS_INT y;              /* some binary powers of x then accumulating */
 219       EMACS_UINT acc, x;  /* Unsigned so that overflow is well defined.  */
 220       Lisp_Object val;
 221
 222       x = XINT (arg1);
 223       y = XINT (arg2);
 224       acc = (y & 1 ? x : 1);
 225
 226       while ((y >>= 1) != 0)
 227         {
 228           x *= x;
 229           if (y & 1)
 230             acc *= x;
 231         }
 232       XSETINT (val, acc);
 233       return val;
 234     }
 235   f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
 236   f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
 237   f3 = pow (f1, f2);
 238   return make_float (f3);
 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 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 302        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 303 This is the same as the exponent of a float.  */)
 304   (Lisp_Object arg)
 305 {
 306   Lisp_Object val;
 307   EMACS_INT value;
 308   double f = extract_float (arg);
 309
 310   if (f == 0.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   XSETINT (val, value);
 322   return val;
 323 }
 324
 325
 326 /* the rounding functions  */
 327
 328 static Lisp_Object
 329 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 330                  double (*double_round) (double),
 331                  EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
 332                  const char *name)
 333 {
 334   CHECK_NUMBER_OR_FLOAT (arg);
 335
 336   if (! NILP (divisor))
 337     {
 338       EMACS_INT i1, i2;
 339
 340       CHECK_NUMBER_OR_FLOAT (divisor);
 341
 342       if (FLOATP (arg) || FLOATP (divisor))
 343         {
 344           double f1, f2;
 345
 346           f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 347           f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
 348           if (! IEEE_FLOATING_POINT && f2 == 0)
 349             xsignal0 (Qarith_error);
 350
 351           f1 = (*double_round) (f1 / f2);
 352           if (FIXNUM_OVERFLOW_P (f1))
 353             xsignal3 (Qrange_error, build_string (name), arg, divisor);
 354           arg = make_number (f1);
 355           return arg;
 356         }
 357
 358       i1 = XINT (arg);
 359       i2 = XINT (divisor);
 360
 361       if (i2 == 0)
 362         xsignal0 (Qarith_error);
 363
 364       XSETINT (arg, (*int_round2) (i1, i2));
 365       return arg;
 366     }
 367
 368   if (FLOATP (arg))
 369     {
 370       double d = (*double_round) (XFLOAT_DATA (arg));
 371       if (FIXNUM_OVERFLOW_P (d))
 372         xsignal2 (Qrange_error, build_string (name), arg);
 373       arg = make_number (d);
 374     }
 375
 376   return arg;
 377 }
 378
 379 /* With C's /, the result is implementation-defined if either operand
 380    is negative, so take care with negative operands in the following
 381    integer functions.  */
 382
 383 static EMACS_INT
 384 ceiling2 (EMACS_INT i1, EMACS_INT i2)
 385 {
 386   return (i2 < 0
 387           ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
 388           : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
 389 }
 390
 391 static EMACS_INT
 392 floor2 (EMACS_INT i1, EMACS_INT i2)
 393 {
 394   return (i2 < 0
 395           ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 396           : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 397 }
 398
 399 static EMACS_INT
 400 truncate2 (EMACS_INT i1, EMACS_INT i2)
 401 {
 402   return (i2 < 0
 403           ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
 404           : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
 405 }
 406
 407 static EMACS_INT
 408 round2 (EMACS_INT i1, EMACS_INT i2)
 409 {
 410   /* The C language's division operator gives us one remainder R, but
 411      we want the remainder R1 on the other side of 0 if R1 is closer
 412      to 0 than R is; because we want to round to even, we also want R1
 413      if R and R1 are the same distance from 0 and if C's quotient is
 414      odd.  */
 415   EMACS_INT q = i1 / i2;
 416   EMACS_INT r = i1 % i2;
 417   EMACS_INT abs_r = eabs (r);
 418   EMACS_INT abs_r1 = eabs (i2) - abs_r;
 419   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 420 }
 421
 422 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 423    if `rint' exists but does not work right.  */
 424 #ifdef HAVE_RINT
 425 #define emacs_rint rint
 426 #else
 427 static double
 428 emacs_rint (double d)
 429 {
 430   double d1 = d + 0.5;
 431   double r = floor (d1);
 432   return r - (r == d1 && fmod (r, 2) != 0);
 433 }
 434 #endif
 435
 436 static double
 437 double_identity (double d)
 438 {
 439   return d;
 440 }
 441
 442 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 443        doc: /* Return the smallest integer no less than ARG.
 444 This rounds the value towards +inf.
 445 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 446   (Lisp_Object arg, Lisp_Object divisor)
 447 {
 448   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 449 }
 450
 451 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 452        doc: /* Return the largest integer no greater than ARG.
 453 This rounds the value towards -inf.
 454 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 455   (Lisp_Object arg, Lisp_Object divisor)
 456 {
 457   return rounding_driver (arg, divisor, floor, floor2, "floor");
 458 }
 459
 460 DEFUN ("round", Fround, Sround, 1, 2, 0,
 461        doc: /* Return the nearest integer to ARG.
 462 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 463
 464 Rounding a value equidistant between two integers may choose the
 465 integer closer to zero, or it may prefer an even integer, depending on
 466 your machine.  For example, \(round 2.5\) can return 3 on some
 467 systems, but 2 on others.  */)
 468   (Lisp_Object arg, Lisp_Object divisor)
 469 {
 470   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 471 }
 472
 473 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 474        doc: /* Truncate a floating point number to an int.
 475 Rounds ARG toward zero.
 476 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 477   (Lisp_Object arg, Lisp_Object divisor)
 478 {
 479   return rounding_driver (arg, divisor, double_identity, truncate2,
 480                           "truncate");
 481 }
 482
 483
 484 Lisp_Object
 485 fmod_float (Lisp_Object x, Lisp_Object y)
 486 {
 487   double f1, f2;
 488
 489   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 490   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 491
 492   f1 = fmod (f1, f2);
 493
 494   /* If the "remainder" comes out with the wrong sign, fix it.  */
 495   if (f2 < 0 ? f1 > 0 : f1 < 0)
 496     f1 += f2;
 497
 498   return make_float (f1);
 499 }
 500 \f
 501 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 502        doc: /* Return the smallest integer no less than ARG, as a float.
 503 \(Round toward +inf.\)  */)
 504   (Lisp_Object arg)
 505 {
 506   double d = extract_float (arg);
 507   d = ceil (d);
 508   return make_float (d);
 509 }
 510
 511 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 512        doc: /* Return the largest integer no greater than ARG, as a float.
 513 \(Round towards -inf.\)  */)
 514   (Lisp_Object arg)
 515 {
 516   double d = extract_float (arg);
 517   d = floor (d);
 518   return make_float (d);
 519 }
 520
 521 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 522        doc: /* Return the nearest integer to ARG, as a float.  */)
 523   (Lisp_Object arg)
 524 {
 525   double d = extract_float (arg);
 526   d = emacs_rint (d);
 527   return make_float (d);
 528 }
 529
 530 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 531        doc: /* Truncate a floating point number to an integral float value.
 532 Rounds the value toward zero.  */)
 533   (Lisp_Object arg)
 534 {
 535   double d = extract_float (arg);
 536   if (d >= 0.0)
 537     d = floor (d);
 538   else
 539     d = ceil (d);
 540   return make_float (d);
 541 }
 542 \f
 543 void
 544 syms_of_floatfns (void)
 545 {
 546   defsubr (&Sacos);
 547   defsubr (&Sasin);
 548   defsubr (&Satan);
 549   defsubr (&Scos);
 550   defsubr (&Ssin);
 551   defsubr (&Stan);
 552   defsubr (&Sisnan);
 553 #ifdef HAVE_COPYSIGN
 554   defsubr (&Scopysign);
 555 #endif
 556   defsubr (&Sfrexp);
 557   defsubr (&Sldexp);
 558   defsubr (&Sfceiling);
 559   defsubr (&Sffloor);
 560   defsubr (&Sfround);
 561   defsubr (&Sftruncate);
 562   defsubr (&Sexp);
 563   defsubr (&Sexpt);
 564   defsubr (&Slog);
 565   defsubr (&Ssqrt);
 566
 567   defsubr (&Sabs);
 568   defsubr (&Sfloat);
 569   defsubr (&Slogb);
 570   defsubr (&Sceiling);
 571   defsubr (&Sfloor);
 572   defsubr (&Sround);
 573   defsubr (&Struncate);
 574 }