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