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, *modf, pow, sin, *sinh, sqrt, tan, *tanh.
  29  */
  30
  31 #include <config.h>
  32
  33 #include "lisp.h"
  34
  35 #include <math.h>
  36
  37 #ifndef isfinite
  38 # define isfinite(x) ((x) - (x) == 0)
  39 #endif
  40 #ifndef isnan
  41 # define isnan(x) ((x) != (x))
  42 #endif
  43
  44 /* Extract a Lisp number as a `double', or signal an error.  */
  45
  46 double
  47 extract_float (Lisp_Object num)
  48 {
  49   CHECK_NUMBER_OR_FLOAT (num);
  50
  51   if (FLOATP (num))
  52     return XFLOAT_DATA (num);
  53   return (double) XINT (num);
  54 }
  55 \f
  56 /* Trig functions.  */
  57
  58 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  59        doc: /* Return the inverse cosine of ARG.  */)
  60   (Lisp_Object arg)
  61 {
  62   double d = extract_float (arg);
  63   d = acos (d);
  64   return make_float (d);
  65 }
  66
  67 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  68        doc: /* Return the inverse sine of ARG.  */)
  69   (Lisp_Object arg)
  70 {
  71   double d = extract_float (arg);
  72   d = asin (d);
  73   return make_float (d);
  74 }
  75
  76 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
  77        doc: /* Return the inverse tangent of the arguments.
  78 If only one argument Y is given, return the inverse tangent of Y.
  79 If two arguments Y and X are given, return the inverse tangent of Y
  80 divided by X, i.e. the angle in radians between the vector (X, Y)
  81 and the x-axis.  */)
  82   (Lisp_Object y, Lisp_Object x)
  83 {
  84   double d = extract_float (y);
  85
  86   if (NILP (x))
  87     d = atan (d);
  88   else
  89     {
  90       double d2 = extract_float (x);
  91       d = atan2 (d, d2);
  92     }
  93   return make_float (d);
  94 }
  95
  96 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
  97        doc: /* Return the cosine of ARG.  */)
  98   (Lisp_Object arg)
  99 {
 100   double d = extract_float (arg);
 101   d = cos (d);
 102   return make_float (d);
 103 }
 104
 105 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 106        doc: /* Return the sine of ARG.  */)
 107   (Lisp_Object arg)
 108 {
 109   double d = extract_float (arg);
 110   d = sin (d);
 111   return make_float (d);
 112 }
 113
 114 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 115        doc: /* Return the tangent of ARG.  */)
 116   (Lisp_Object arg)
 117 {
 118   double d = extract_float (arg);
 119   d = tan (d);
 120   return make_float (d);
 121 }
 122
 123 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
 124        doc: /* Return non nil iff argument X is a NaN.  */)
 125   (Lisp_Object x)
 126 {
 127   CHECK_FLOAT (x);
 128   return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
 129 }
 130
 131 #ifdef HAVE_COPYSIGN
 132 DEFUN ("copysign", Fcopysign, Scopysign, 2, 2, 0,
 133        doc: /* Copy sign of X2 to value of X1, and return the result.
 134 Cause an error if X1 or X2 is not a float.  */)
 135   (Lisp_Object x1, Lisp_Object x2)
 136 {
 137   double f1, f2;
 138
 139   CHECK_FLOAT (x1);
 140   CHECK_FLOAT (x2);
 141
 142   f1 = XFLOAT_DATA (x1);
 143   f2 = XFLOAT_DATA (x2);
 144
 145   return make_float (copysign (f1, f2));
 146 }
 147 #endif
 148
 149 DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
 150        doc: /* Get significand and exponent of a floating point number.
 151 Breaks the floating point number X into its binary significand SGNFCAND
 152 \(a floating point value between 0.5 (included) and 1.0 (excluded))
 153 and an integral exponent EXP for 2, such that:
 154
 155   X = SGNFCAND * 2^EXP
 156
 157 The function returns the cons cell (SGNFCAND . EXP).
 158 If X is zero, both parts (SGNFCAND and EXP) are zero.  */)
 159   (Lisp_Object x)
 160 {
 161   double f = XFLOATINT (x);
 162   int exponent;
 163   double sgnfcand = frexp (f, &exponent);
 164   return Fcons (make_float (sgnfcand), make_number (exponent));
 165 }
 166
 167 DEFUN ("ldexp", Fldexp, Sldexp, 1, 2, 0,
 168        doc: /* Construct number X from significand SGNFCAND and exponent EXP.
 169 Returns the floating point value resulting from multiplying SGNFCAND
 170 (the significand) by 2 raised to the power of EXP (the exponent).   */)
 171   (Lisp_Object sgnfcand, Lisp_Object exponent)
 172 {
 173   CHECK_NUMBER (exponent);
 174   return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exponent)));
 175 }
 176 \f
 177 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 178        doc: /* Return the exponential base e of ARG.  */)
 179   (Lisp_Object arg)
 180 {
 181   double d = extract_float (arg);
 182   d = exp (d);
 183   return make_float (d);
 184 }
 185
 186 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 187        doc: /* Return the exponential ARG1 ** ARG2.  */)
 188   (Lisp_Object arg1, Lisp_Object arg2)
 189 {
 190   double f1, f2, f3;
 191
 192   CHECK_NUMBER_OR_FLOAT (arg1);
 193   CHECK_NUMBER_OR_FLOAT (arg2);
 194   if (INTEGERP (arg1)     /* common lisp spec */
 195       && INTEGERP (arg2)   /* don't promote, if both are ints, and */
 196       && 0 <= XINT (arg2)) /* we are sure the result is not fractional */
 197     {                           /* this can be improved by pre-calculating */
 198       EMACS_INT y;              /* some binary powers of x then accumulating */
 199       EMACS_UINT acc, x;  /* Unsigned so that overflow is well defined.  */
 200       Lisp_Object val;
 201
 202       x = XINT (arg1);
 203       y = XINT (arg2);
 204       acc = (y & 1 ? x : 1);
 205
 206       while ((y >>= 1) != 0)
 207         {
 208           x *= x;
 209           if (y & 1)
 210             acc *= x;
 211         }
 212       XSETINT (val, acc);
 213       return val;
 214     }
 215   f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
 216   f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
 217   f3 = pow (f1, f2);
 218   return make_float (f3);
 219 }
 220
 221 DEFUN ("log", Flog, Slog, 1, 2, 0,
 222        doc: /* Return the natural logarithm of ARG.
 223 If the optional argument BASE is given, return log ARG using that base.  */)
 224   (Lisp_Object arg, Lisp_Object base)
 225 {
 226   double d = extract_float (arg);
 227
 228   if (NILP (base))
 229     d = log (d);
 230   else
 231     {
 232       double b = extract_float (base);
 233
 234       if (b == 10.0)
 235         d = log10 (d);
 236       else
 237         d = log (d) / log (b);
 238     }
 239   return make_float (d);
 240 }
 241
 242 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
 243        doc: /* Return the logarithm base 10 of ARG.  */)
 244   (Lisp_Object arg)
 245 {
 246   double d = extract_float (arg);
 247   d = log10 (d);
 248   return make_float (d);
 249 }
 250
 251 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 252        doc: /* Return the square root of ARG.  */)
 253   (Lisp_Object arg)
 254 {
 255   double d = extract_float (arg);
 256   d = sqrt (d);
 257   return make_float (d);
 258 }
 259 \f
 260 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 261        doc: /* Return the absolute value of ARG.  */)
 262   (register Lisp_Object arg)
 263 {
 264   CHECK_NUMBER_OR_FLOAT (arg);
 265
 266   if (FLOATP (arg))
 267     arg = make_float (fabs (XFLOAT_DATA (arg)));
 268   else if (XINT (arg) < 0)
 269     XSETINT (arg, - XINT (arg));
 270
 271   return arg;
 272 }
 273
 274 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 275        doc: /* Return the floating point number equal to ARG.  */)
 276   (register Lisp_Object arg)
 277 {
 278   CHECK_NUMBER_OR_FLOAT (arg);
 279
 280   if (INTEGERP (arg))
 281     return make_float ((double) XINT (arg));
 282   else                          /* give 'em the same float back */
 283     return arg;
 284 }
 285
 286 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 287        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 288 This is the same as the exponent of a float.  */)
 289   (Lisp_Object arg)
 290 {
 291   Lisp_Object val;
 292   EMACS_INT value;
 293   double f = extract_float (arg);
 294
 295   if (f == 0.0)
 296     value = MOST_NEGATIVE_FIXNUM;
 297   else if (isfinite (f))
 298     {
 299       int ivalue;
 300       frexp (f, &ivalue);
 301       value = ivalue - 1;
 302     }
 303   else
 304     value = MOST_POSITIVE_FIXNUM;
 305
 306   XSETINT (val, value);
 307   return val;
 308 }
 309
 310
 311 /* the rounding functions  */
 312
 313 static Lisp_Object
 314 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 315                  double (*double_round) (double),
 316                  EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
 317                  const char *name)
 318 {
 319   CHECK_NUMBER_OR_FLOAT (arg);
 320
 321   if (! NILP (divisor))
 322     {
 323       EMACS_INT i1, i2;
 324
 325       CHECK_NUMBER_OR_FLOAT (divisor);
 326
 327       if (FLOATP (arg) || FLOATP (divisor))
 328         {
 329           double f1, f2;
 330
 331           f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 332           f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
 333           if (! IEEE_FLOATING_POINT && f2 == 0)
 334             xsignal0 (Qarith_error);
 335
 336           f1 = (*double_round) (f1 / f2);
 337           if (FIXNUM_OVERFLOW_P (f1))
 338             xsignal3 (Qrange_error, build_string (name), arg, divisor);
 339           arg = make_number (f1);
 340           return arg;
 341         }
 342
 343       i1 = XINT (arg);
 344       i2 = XINT (divisor);
 345
 346       if (i2 == 0)
 347         xsignal0 (Qarith_error);
 348
 349       XSETINT (arg, (*int_round2) (i1, i2));
 350       return arg;
 351     }
 352
 353   if (FLOATP (arg))
 354     {
 355       double d = (*double_round) (XFLOAT_DATA (arg));
 356       if (FIXNUM_OVERFLOW_P (d))
 357         xsignal2 (Qrange_error, build_string (name), arg);
 358       arg = make_number (d);
 359     }
 360
 361   return arg;
 362 }
 363
 364 /* With C's /, the result is implementation-defined if either operand
 365    is negative, so take care with negative operands in the following
 366    integer functions.  */
 367
 368 static EMACS_INT
 369 ceiling2 (EMACS_INT i1, EMACS_INT i2)
 370 {
 371   return (i2 < 0
 372           ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
 373           : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
 374 }
 375
 376 static EMACS_INT
 377 floor2 (EMACS_INT i1, EMACS_INT i2)
 378 {
 379   return (i2 < 0
 380           ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 381           : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 382 }
 383
 384 static EMACS_INT
 385 truncate2 (EMACS_INT i1, EMACS_INT i2)
 386 {
 387   return (i2 < 0
 388           ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
 389           : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
 390 }
 391
 392 static EMACS_INT
 393 round2 (EMACS_INT i1, EMACS_INT i2)
 394 {
 395   /* The C language's division operator gives us one remainder R, but
 396      we want the remainder R1 on the other side of 0 if R1 is closer
 397      to 0 than R is; because we want to round to even, we also want R1
 398      if R and R1 are the same distance from 0 and if C's quotient is
 399      odd.  */
 400   EMACS_INT q = i1 / i2;
 401   EMACS_INT r = i1 % i2;
 402   EMACS_INT abs_r = eabs (r);
 403   EMACS_INT abs_r1 = eabs (i2) - abs_r;
 404   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 405 }
 406
 407 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 408    if `rint' exists but does not work right.  */
 409 #ifdef HAVE_RINT
 410 #define emacs_rint rint
 411 #else
 412 static double
 413 emacs_rint (double d)
 414 {
 415   return floor (d + 0.5);
 416 }
 417 #endif
 418
 419 static double
 420 double_identity (double d)
 421 {
 422   return d;
 423 }
 424
 425 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 426        doc: /* Return the smallest integer no less than ARG.
 427 This rounds the value towards +inf.
 428 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 429   (Lisp_Object arg, Lisp_Object divisor)
 430 {
 431   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 432 }
 433
 434 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 435        doc: /* Return the largest integer no greater than ARG.
 436 This rounds the value towards -inf.
 437 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 438   (Lisp_Object arg, Lisp_Object divisor)
 439 {
 440   return rounding_driver (arg, divisor, floor, floor2, "floor");
 441 }
 442
 443 DEFUN ("round", Fround, Sround, 1, 2, 0,
 444        doc: /* Return the nearest integer to ARG.
 445 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 446
 447 Rounding a value equidistant between two integers may choose the
 448 integer closer to zero, or it may prefer an even integer, depending on
 449 your machine.  For example, \(round 2.5\) can return 3 on some
 450 systems, but 2 on others.  */)
 451   (Lisp_Object arg, Lisp_Object divisor)
 452 {
 453   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 454 }
 455
 456 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 457        doc: /* Truncate a floating point number to an int.
 458 Rounds ARG toward zero.
 459 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 460   (Lisp_Object arg, Lisp_Object divisor)
 461 {
 462   return rounding_driver (arg, divisor, double_identity, truncate2,
 463                           "truncate");
 464 }
 465
 466
 467 Lisp_Object
 468 fmod_float (Lisp_Object x, Lisp_Object y)
 469 {
 470   double f1, f2;
 471
 472   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 473   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 474
 475   f1 = fmod (f1, f2);
 476
 477   /* If the "remainder" comes out with the wrong sign, fix it.  */
 478   if (f2 < 0 ? 0 < f1 : f1 < 0)
 479     f1 += f2;
 480
 481   return make_float (f1);
 482 }
 483 \f
 484 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 485        doc: /* Return the smallest integer no less than ARG, as a float.
 486 \(Round toward +inf.\)  */)
 487   (Lisp_Object arg)
 488 {
 489   double d = extract_float (arg);
 490   d = ceil (d);
 491   return make_float (d);
 492 }
 493
 494 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 495        doc: /* Return the largest integer no greater than ARG, as a float.
 496 \(Round towards -inf.\)  */)
 497   (Lisp_Object arg)
 498 {
 499   double d = extract_float (arg);
 500   d = floor (d);
 501   return make_float (d);
 502 }
 503
 504 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 505        doc: /* Return the nearest integer to ARG, as a float.  */)
 506   (Lisp_Object arg)
 507 {
 508   double d = extract_float (arg);
 509   d = emacs_rint (d);
 510   return make_float (d);
 511 }
 512
 513 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 514        doc: /* Truncate a floating point number to an integral float value.
 515 Rounds the value toward zero.  */)
 516   (Lisp_Object arg)
 517 {
 518   double d = extract_float (arg);
 519   if (d >= 0.0)
 520     d = floor (d);
 521   else
 522     d = ceil (d);
 523   return make_float (d);
 524 }
 525 \f
 526 void
 527 syms_of_floatfns (void)
 528 {
 529   defsubr (&Sacos);
 530   defsubr (&Sasin);
 531   defsubr (&Satan);
 532   defsubr (&Scos);
 533   defsubr (&Ssin);
 534   defsubr (&Stan);
 535   defsubr (&Sisnan);
 536 #ifdef HAVE_COPYSIGN
 537   defsubr (&Scopysign);
 538 #endif
 539   defsubr (&Sfrexp);
 540   defsubr (&Sldexp);
 541   defsubr (&Sfceiling);
 542   defsubr (&Sffloor);
 543   defsubr (&Sfround);
 544   defsubr (&Sftruncate);
 545   defsubr (&Sexp);
 546   defsubr (&Sexpt);
 547   defsubr (&Slog);
 548   defsubr (&Slog10);
 549   defsubr (&Ssqrt);
 550
 551   defsubr (&Sabs);
 552   defsubr (&Sfloat);
 553   defsubr (&Slogb);
 554   defsubr (&Sceiling);
 555   defsubr (&Sfloor);
 556   defsubr (&Sround);
 557   defsubr (&Struncate);
 558 }