src/floatfns.c

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