src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2
   3 Copyright (C) 1988, 1993-1994, 1999, 2001-2013 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 #ifndef isfinite
  50 # define isfinite(x) ((x) - (x) == 0)
  51 #endif
  52 #ifndef isnan
  53 # define isnan(x) ((x) != (x))
  54 #endif
  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 iff 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   return floor (d + 0.5);
 431 }
 432 #endif
 433
 434 static double
 435 double_identity (double d)
 436 {
 437   return d;
 438 }
 439
 440 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 441        doc: /* Return the smallest integer no less than ARG.
 442 This rounds the value towards +inf.
 443 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 444   (Lisp_Object arg, Lisp_Object divisor)
 445 {
 446   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 447 }
 448
 449 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 450        doc: /* Return the largest integer no greater than ARG.
 451 This rounds the value towards -inf.
 452 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 453   (Lisp_Object arg, Lisp_Object divisor)
 454 {
 455   return rounding_driver (arg, divisor, floor, floor2, "floor");
 456 }
 457
 458 DEFUN ("round", Fround, Sround, 1, 2, 0,
 459        doc: /* Return the nearest integer to ARG.
 460 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 461
 462 Rounding a value equidistant between two integers may choose the
 463 integer closer to zero, or it may prefer an even integer, depending on
 464 your machine.  For example, \(round 2.5\) can return 3 on some
 465 systems, but 2 on others.  */)
 466   (Lisp_Object arg, Lisp_Object divisor)
 467 {
 468   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 469 }
 470
 471 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 472        doc: /* Truncate a floating point number to an int.
 473 Rounds ARG toward zero.
 474 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 475   (Lisp_Object arg, Lisp_Object divisor)
 476 {
 477   return rounding_driver (arg, divisor, double_identity, truncate2,
 478                           "truncate");
 479 }
 480
 481
 482 Lisp_Object
 483 fmod_float (Lisp_Object x, Lisp_Object y)
 484 {
 485   double f1, f2;
 486
 487   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 488   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 489
 490   f1 = fmod (f1, f2);
 491
 492   /* If the "remainder" comes out with the wrong sign, fix it.  */
 493   if (f2 < 0 ? f1 > 0 : f1 < 0)
 494     f1 += f2;
 495
 496   return make_float (f1);
 497 }
 498 \f
 499 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 500        doc: /* Return the smallest integer no less than ARG, as a float.
 501 \(Round toward +inf.\)  */)
 502   (Lisp_Object arg)
 503 {
 504   double d = extract_float (arg);
 505   d = ceil (d);
 506   return make_float (d);
 507 }
 508
 509 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 510        doc: /* Return the largest integer no greater than ARG, as a float.
 511 \(Round towards -inf.\)  */)
 512   (Lisp_Object arg)
 513 {
 514   double d = extract_float (arg);
 515   d = floor (d);
 516   return make_float (d);
 517 }
 518
 519 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 520        doc: /* Return the nearest integer to ARG, as a float.  */)
 521   (Lisp_Object arg)
 522 {
 523   double d = extract_float (arg);
 524   d = emacs_rint (d);
 525   return make_float (d);
 526 }
 527
 528 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 529        doc: /* Truncate a floating point number to an integral float value.
 530 Rounds the value toward zero.  */)
 531   (Lisp_Object arg)
 532 {
 533   double d = extract_float (arg);
 534   if (d >= 0.0)
 535     d = floor (d);
 536   else
 537     d = ceil (d);
 538   return make_float (d);
 539 }
 540 \f
 541 void
 542 syms_of_floatfns (void)
 543 {
 544   defsubr (&Sacos);
 545   defsubr (&Sasin);
 546   defsubr (&Satan);
 547   defsubr (&Scos);
 548   defsubr (&Ssin);
 549   defsubr (&Stan);
 550   defsubr (&Sisnan);
 551 #ifdef HAVE_COPYSIGN
 552   defsubr (&Scopysign);
 553 #endif
 554   defsubr (&Sfrexp);
 555   defsubr (&Sldexp);
 556   defsubr (&Sfceiling);
 557   defsubr (&Sffloor);
 558   defsubr (&Sfround);
 559   defsubr (&Sftruncate);
 560   defsubr (&Sexp);
 561   defsubr (&Sexpt);
 562   defsubr (&Slog);
 563   defsubr (&Ssqrt);
 564
 565   defsubr (&Sabs);
 566   defsubr (&Sfloat);
 567   defsubr (&Slogb);
 568   defsubr (&Sceiling);
 569   defsubr (&Sfloor);
 570   defsubr (&Sround);
 571   defsubr (&Struncate);
 572 }