src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2
   3 Copyright (C) 1988, 1993-1994, 1999, 2001-2017 Free Software Foundation,
   4 Inc.
   5
   6 Author: Wolfgang Rupprecht (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 (at
  13 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 #include <count-leading-zeros.h>
  49
  50 /* 'isfinite' and 'isnan' cause build failures on Solaris 10 with the
  51    bundled GCC in c99 mode.  Work around the bugs with simple
  52    implementations that are good enough.  */
  53 #undef isfinite
  54 #define isfinite(x) ((x) - (x) == 0)
  55 #undef isnan
  56 #define isnan(x) ((x) != (x))
  57
  58 /* Check that X is a floating point number.  */
  59
  60 static void
  61 CHECK_FLOAT (Lisp_Object x)
  62 {
  63   CHECK_TYPE (FLOATP (x), Qfloatp, x);
  64 }
  65
  66 /* Extract a Lisp number as a `double', or signal an error.  */
  67
  68 double
  69 extract_float (Lisp_Object num)
  70 {
  71   CHECK_NUMBER_OR_FLOAT (num);
  72   return XFLOATINT (num);
  73 }
  74 \f
  75 /* Trig functions.  */
  76
  77 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  78        doc: /* Return the inverse cosine of ARG.  */)
  79   (Lisp_Object arg)
  80 {
  81   double d = extract_float (arg);
  82   d = acos (d);
  83   return make_float (d);
  84 }
  85
  86 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  87        doc: /* Return the inverse sine of ARG.  */)
  88   (Lisp_Object arg)
  89 {
  90   double d = extract_float (arg);
  91   d = asin (d);
  92   return make_float (d);
  93 }
  94
  95 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
  96        doc: /* Return the inverse tangent of the arguments.
  97 If only one argument Y is given, return the inverse tangent of Y.
  98 If two arguments Y and X are given, return the inverse tangent of Y
  99 divided by X, i.e. the angle in radians between the vector (X, Y)
 100 and the x-axis.  */)
 101   (Lisp_Object y, Lisp_Object x)
 102 {
 103   double d = extract_float (y);
 104
 105   if (NILP (x))
 106     d = atan (d);
 107   else
 108     {
 109       double d2 = extract_float (x);
 110       d = atan2 (d, d2);
 111     }
 112   return make_float (d);
 113 }
 114
 115 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 116        doc: /* Return the cosine of ARG.  */)
 117   (Lisp_Object arg)
 118 {
 119   double d = extract_float (arg);
 120   d = cos (d);
 121   return make_float (d);
 122 }
 123
 124 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 125        doc: /* Return the sine of ARG.  */)
 126   (Lisp_Object arg)
 127 {
 128   double d = extract_float (arg);
 129   d = sin (d);
 130   return make_float (d);
 131 }
 132
 133 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 134        doc: /* Return the tangent of ARG.  */)
 135   (Lisp_Object arg)
 136 {
 137   double d = extract_float (arg);
 138   d = tan (d);
 139   return make_float (d);
 140 }
 141
 142 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
 143        doc: /* Return non nil if argument X is a NaN.  */)
 144   (Lisp_Object x)
 145 {
 146   CHECK_FLOAT (x);
 147   return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
 148 }
 149
 150 #ifdef HAVE_COPYSIGN
 151 DEFUN ("copysign", Fcopysign, Scopysign, 2, 2, 0,
 152        doc: /* Copy sign of X2 to value of X1, and return the result.
 153 Cause an error if X1 or X2 is not a float.  */)
 154   (Lisp_Object x1, Lisp_Object x2)
 155 {
 156   double f1, f2;
 157
 158   CHECK_FLOAT (x1);
 159   CHECK_FLOAT (x2);
 160
 161   f1 = XFLOAT_DATA (x1);
 162   f2 = XFLOAT_DATA (x2);
 163
 164   return make_float (copysign (f1, f2));
 165 }
 166 #endif
 167
 168 DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
 169        doc: /* Get significand and exponent of a floating point number.
 170 Breaks the floating point number X into its binary significand SGNFCAND
 171 \(a floating point value between 0.5 (included) and 1.0 (excluded))
 172 and an integral exponent EXP for 2, such that:
 173
 174   X = SGNFCAND * 2^EXP
 175
 176 The function returns the cons cell (SGNFCAND . EXP).
 177 If X is zero, both parts (SGNFCAND and EXP) are zero.  */)
 178   (Lisp_Object x)
 179 {
 180   double f = extract_float (x);
 181   int exponent;
 182   double sgnfcand = frexp (f, &exponent);
 183   return Fcons (make_float (sgnfcand), make_number (exponent));
 184 }
 185
 186 DEFUN ("ldexp", Fldexp, Sldexp, 2, 2, 0,
 187        doc: /* Return SGNFCAND * 2**EXPONENT, as a floating point number.
 188 EXPONENT must be an integer.   */)
 189   (Lisp_Object sgnfcand, Lisp_Object exponent)
 190 {
 191   CHECK_NUMBER (exponent);
 192   int e = min (max (INT_MIN, XINT (exponent)), INT_MAX);
 193   return make_float (ldexp (extract_float (sgnfcand), e));
 194 }
 195 \f
 196 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 197        doc: /* Return the exponential base e of ARG.  */)
 198   (Lisp_Object arg)
 199 {
 200   double d = extract_float (arg);
 201   d = exp (d);
 202   return make_float (d);
 203 }
 204
 205 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 206        doc: /* Return the exponential ARG1 ** ARG2.  */)
 207   (Lisp_Object arg1, Lisp_Object arg2)
 208 {
 209   CHECK_NUMBER_OR_FLOAT (arg1);
 210   CHECK_NUMBER_OR_FLOAT (arg2);
 211   if (INTEGERP (arg1)     /* common lisp spec */
 212       && INTEGERP (arg2)   /* don't promote, if both are ints, and */
 213       && XINT (arg2) >= 0) /* we are sure the result is not fractional */
 214     {                           /* this can be improved by pre-calculating */
 215       EMACS_INT y;              /* some binary powers of x then accumulating */
 216       EMACS_UINT acc, x;  /* Unsigned so that overflow is well defined.  */
 217       Lisp_Object val;
 218
 219       x = XINT (arg1);
 220       y = XINT (arg2);
 221       acc = (y & 1 ? x : 1);
 222
 223       while ((y >>= 1) != 0)
 224         {
 225           x *= x;
 226           if (y & 1)
 227             acc *= x;
 228         }
 229       XSETINT (val, acc);
 230       return val;
 231     }
 232   return make_float (pow (XFLOATINT (arg1), XFLOATINT (arg2)));
 233 }
 234
 235 DEFUN ("log", Flog, Slog, 1, 2, 0,
 236        doc: /* Return the natural logarithm of ARG.
 237 If the optional argument BASE is given, return log ARG using that base.  */)
 238   (Lisp_Object arg, Lisp_Object base)
 239 {
 240   double d = extract_float (arg);
 241
 242   if (NILP (base))
 243     d = log (d);
 244   else
 245     {
 246       double b = extract_float (base);
 247
 248       if (b == 10.0)
 249         d = log10 (d);
 250 #if HAVE_LOG2
 251       else if (b == 2.0)
 252         d = log2 (d);
 253 #endif
 254       else
 255         d = log (d) / log (b);
 256     }
 257   return make_float (d);
 258 }
 259
 260 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 261        doc: /* Return the square root of ARG.  */)
 262   (Lisp_Object arg)
 263 {
 264   double d = extract_float (arg);
 265   d = sqrt (d);
 266   return make_float (d);
 267 }
 268 \f
 269 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 270        doc: /* Return the absolute value of ARG.  */)
 271   (register Lisp_Object arg)
 272 {
 273   CHECK_NUMBER_OR_FLOAT (arg);
 274
 275   if (FLOATP (arg))
 276     arg = make_float (fabs (XFLOAT_DATA (arg)));
 277   else if (XINT (arg) < 0)
 278     XSETINT (arg, - XINT (arg));
 279
 280   return arg;
 281 }
 282
 283 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 284        doc: /* Return the floating point number equal to ARG.  */)
 285   (register Lisp_Object arg)
 286 {
 287   CHECK_NUMBER_OR_FLOAT (arg);
 288
 289   if (INTEGERP (arg))
 290     return make_float ((double) XINT (arg));
 291   else                          /* give 'em the same float back */
 292     return arg;
 293 }
 294
 295 static int
 296 ecount_leading_zeros (EMACS_UINT x)
 297 {
 298   return (EMACS_UINT_WIDTH == UINT_WIDTH ? count_leading_zeros (x)
 299           : EMACS_UINT_WIDTH == ULONG_WIDTH ? count_leading_zeros_l (x)
 300           : count_leading_zeros_ll (x));
 301 }
 302
 303 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 304        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 305 This is the same as the exponent of a float.  */)
 306   (Lisp_Object arg)
 307 {
 308   EMACS_INT value;
 309   CHECK_NUMBER_OR_FLOAT (arg);
 310
 311   if (FLOATP (arg))
 312     {
 313       double f = XFLOAT_DATA (arg);
 314
 315       if (f == 0)
 316         value = MOST_NEGATIVE_FIXNUM;
 317       else if (isfinite (f))
 318         {
 319           int ivalue;
 320           frexp (f, &ivalue);
 321           value = ivalue - 1;
 322         }
 323       else
 324         value = MOST_POSITIVE_FIXNUM;
 325     }
 326   else
 327     {
 328       EMACS_INT i = eabs (XINT (arg));
 329       value = (i == 0
 330                ? MOST_NEGATIVE_FIXNUM
 331                : EMACS_UINT_WIDTH - 1 - ecount_leading_zeros (i));
 332     }
 333
 334   return make_number (value);
 335 }
 336
 337
 338 /* the rounding functions  */
 339
 340 static Lisp_Object
 341 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 342                  double (*double_round) (double),
 343                  EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
 344                  const char *name)
 345 {
 346   CHECK_NUMBER_OR_FLOAT (arg);
 347
 348   double d;
 349   if (NILP (divisor))
 350     {
 351       if (! FLOATP (arg))
 352         return arg;
 353       d = XFLOAT_DATA (arg);
 354     }
 355   else
 356     {
 357       CHECK_NUMBER_OR_FLOAT (divisor);
 358       if (!FLOATP (arg) && !FLOATP (divisor))
 359         {
 360           if (XINT (divisor) == 0)
 361             xsignal0 (Qarith_error);
 362           return make_number (int_round2 (XINT (arg), XINT (divisor)));
 363         }
 364
 365       double f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 366       double f2 = FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor);
 367       if (! IEEE_FLOATING_POINT && f2 == 0)
 368         xsignal0 (Qarith_error);
 369       d = f1 / f2;
 370     }
 371
 372   /* Round, coarsely test for fixnum overflow before converting to
 373      EMACS_INT (to avoid undefined C behavior), and then exactly test
 374      for overflow after converting (as FIXNUM_OVERFLOW_P is inaccurate
 375      on floats).  */
 376   double dr = double_round (d);
 377   if (fabs (dr) < 2 * (MOST_POSITIVE_FIXNUM + 1))
 378     {
 379       EMACS_INT ir = dr;
 380       if (! FIXNUM_OVERFLOW_P (ir))
 381         return make_number (ir);
 382     }
 383   xsignal2 (Qrange_error, build_string (name), arg);
 384 }
 385
 386 static EMACS_INT
 387 ceiling2 (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 floor2 (EMACS_INT i1, EMACS_INT i2)
 394 {
 395   return i1 / i2 - ((i1 % i2 != 0) & ((i1 < 0) != (i2 < 0)));
 396 }
 397
 398 static EMACS_INT
 399 truncate2 (EMACS_INT i1, EMACS_INT i2)
 400 {
 401   return i1 / i2;
 402 }
 403
 404 static EMACS_INT
 405 round2 (EMACS_INT i1, EMACS_INT i2)
 406 {
 407   /* The C language's division operator gives us one remainder R, but
 408      we want the remainder R1 on the other side of 0 if R1 is closer
 409      to 0 than R is; because we want to round to even, we also want R1
 410      if R and R1 are the same distance from 0 and if C's quotient is
 411      odd.  */
 412   EMACS_INT q = i1 / i2;
 413   EMACS_INT r = i1 % i2;
 414   EMACS_INT abs_r = eabs (r);
 415   EMACS_INT abs_r1 = eabs (i2) - abs_r;
 416   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 417 }
 418
 419 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 420    if `rint' exists but does not work right.  */
 421 #ifdef HAVE_RINT
 422 #define emacs_rint rint
 423 #else
 424 static double
 425 emacs_rint (double d)
 426 {
 427   double d1 = d + 0.5;
 428   double r = floor (d1);
 429   return r - (r == d1 && fmod (r, 2) != 0);
 430 }
 431 #endif
 432
 433 #ifdef HAVE_TRUNC
 434 #define emacs_trunc trunc
 435 #else
 436 static double
 437 emacs_trunc (double d)
 438 {
 439   return (d < 0 ? ceil : floor) (d);
 440 }
 441 #endif
 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, emacs_trunc, 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 }