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