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 static EMACS_INT
 381 ceiling2 (EMACS_INT i1, EMACS_INT i2)
 382 {
 383   return i1 / i2 + ((i1 % i2 != 0) & ((i1 < 0) == (i2 < 0)));
 384 }
 385
 386 static EMACS_INT
 387 floor2 (EMACS_INT i1, EMACS_INT i2)
 388 {
 389   return i1 / i2 - ((i1 % i2 != 0) & ((i1 < 0) != (i2 < 0)));
 390 }
 391
 392 static EMACS_INT
 393 truncate2 (EMACS_INT i1, EMACS_INT i2)
 394 {
 395   return i1 / i2;
 396 }
 397
 398 static EMACS_INT
 399 round2 (EMACS_INT i1, EMACS_INT i2)
 400 {
 401   /* The C language's division operator gives us one remainder R, but
 402      we want the remainder R1 on the other side of 0 if R1 is closer
 403      to 0 than R is; because we want to round to even, we also want R1
 404      if R and R1 are the same distance from 0 and if C's quotient is
 405      odd.  */
 406   EMACS_INT q = i1 / i2;
 407   EMACS_INT r = i1 % i2;
 408   EMACS_INT abs_r = eabs (r);
 409   EMACS_INT abs_r1 = eabs (i2) - abs_r;
 410   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 411 }
 412
 413 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 414    if `rint' exists but does not work right.  */
 415 #ifdef HAVE_RINT
 416 #define emacs_rint rint
 417 #else
 418 static double
 419 emacs_rint (double d)
 420 {
 421   double d1 = d + 0.5;
 422   double r = floor (d1);
 423   return r - (r == d1 && fmod (r, 2) != 0);
 424 }
 425 #endif
 426
 427 static double
 428 double_identity (double d)
 429 {
 430   return d;
 431 }
 432
 433 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 434        doc: /* Return the smallest integer no less than ARG.
 435 This rounds the value towards +inf.
 436 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 437   (Lisp_Object arg, Lisp_Object divisor)
 438 {
 439   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 440 }
 441
 442 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 443        doc: /* Return the largest integer no greater than ARG.
 444 This rounds the value towards -inf.
 445 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 446   (Lisp_Object arg, Lisp_Object divisor)
 447 {
 448   return rounding_driver (arg, divisor, floor, floor2, "floor");
 449 }
 450
 451 DEFUN ("round", Fround, Sround, 1, 2, 0,
 452        doc: /* Return the nearest integer to ARG.
 453 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 454
 455 Rounding a value equidistant between two integers may choose the
 456 integer closer to zero, or it may prefer an even integer, depending on
 457 your machine.  For example, (round 2.5) can return 3 on some
 458 systems, but 2 on others.  */)
 459   (Lisp_Object arg, Lisp_Object divisor)
 460 {
 461   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 462 }
 463
 464 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 465        doc: /* Truncate a floating point number to an int.
 466 Rounds ARG toward zero.
 467 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 468   (Lisp_Object arg, Lisp_Object divisor)
 469 {
 470   return rounding_driver (arg, divisor, double_identity, truncate2,
 471                           "truncate");
 472 }
 473
 474
 475 Lisp_Object
 476 fmod_float (Lisp_Object x, Lisp_Object y)
 477 {
 478   double f1, f2;
 479
 480   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 481   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 482
 483   f1 = fmod (f1, f2);
 484
 485   /* If the "remainder" comes out with the wrong sign, fix it.  */
 486   if (f2 < 0 ? f1 > 0 : f1 < 0)
 487     f1 += f2;
 488
 489   return make_float (f1);
 490 }
 491 \f
 492 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 493        doc: /* Return the smallest integer no less than ARG, as a float.
 494 (Round toward +inf.)  */)
 495   (Lisp_Object arg)
 496 {
 497   double d = extract_float (arg);
 498   d = ceil (d);
 499   return make_float (d);
 500 }
 501
 502 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 503        doc: /* Return the largest integer no greater than ARG, as a float.
 504 (Round towards -inf.)  */)
 505   (Lisp_Object arg)
 506 {
 507   double d = extract_float (arg);
 508   d = floor (d);
 509   return make_float (d);
 510 }
 511
 512 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 513        doc: /* Return the nearest integer to ARG, as a float.  */)
 514   (Lisp_Object arg)
 515 {
 516   double d = extract_float (arg);
 517   d = emacs_rint (d);
 518   return make_float (d);
 519 }
 520
 521 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 522        doc: /* Truncate a floating point number to an integral float value.
 523 Rounds the value toward zero.  */)
 524   (Lisp_Object arg)
 525 {
 526   double d = extract_float (arg);
 527   if (d >= 0.0)
 528     d = floor (d);
 529   else
 530     d = ceil (d);
 531   return make_float (d);
 532 }
 533 \f
 534 void
 535 syms_of_floatfns (void)
 536 {
 537   defsubr (&Sacos);
 538   defsubr (&Sasin);
 539   defsubr (&Satan);
 540   defsubr (&Scos);
 541   defsubr (&Ssin);
 542   defsubr (&Stan);
 543   defsubr (&Sisnan);
 544 #ifdef HAVE_COPYSIGN
 545   defsubr (&Scopysign);
 546 #endif
 547   defsubr (&Sfrexp);
 548   defsubr (&Sldexp);
 549   defsubr (&Sfceiling);
 550   defsubr (&Sffloor);
 551   defsubr (&Sfround);
 552   defsubr (&Sftruncate);
 553   defsubr (&Sexp);
 554   defsubr (&Sexpt);
 555   defsubr (&Slog);
 556   defsubr (&Ssqrt);
 557
 558   defsubr (&Sabs);
 559   defsubr (&Sfloat);
 560   defsubr (&Slogb);
 561   defsubr (&Sceiling);
 562   defsubr (&Sfloor);
 563   defsubr (&Sround);
 564   defsubr (&Struncate);
 565 }