code.delx.au - gnu-emacs/blob - 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    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, *logb (approximately),
  37    lrint/llrint, lround/llround, nan, nearbyint, nextafter,
  38    nexttoward, remainder, remquo, *rint, round, scalbln, scalbn,
  39    signbit, tgamma, trunc.
  40  */
  41
  42 #include <config.h>
  43
  44 #include "lisp.h"
  45
  46 #include <math.h>
  47
  48 #ifndef isfinite
  49 # define isfinite(x) ((x) - (x) == 0)
  50 #endif
  51 #ifndef isnan
  52 # define isnan(x) ((x) != (x))
  53 #endif
  54
  55 /* Check that X is a floating point number.  */
  56
  57 static void
  58 CHECK_FLOAT (Lisp_Object x)
  59 {
  60   CHECK_TYPE (FLOATP (x), Qfloatp, x);
  61 }
  62
  63 /* Extract a Lisp number as a `double', or signal an error.  */
  64
  65 double
  66 extract_float (Lisp_Object num)
  67 {
  68   CHECK_NUMBER_OR_FLOAT (num);
  69
  70   if (FLOATP (num))
  71     return XFLOAT_DATA (num);
  72   return (double) XINT (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 iff 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 = XFLOATINT (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, 1, 2, 0,
 187        doc: /* Construct number X from significand SGNFCAND and exponent EXP.
 188 Returns the floating point value resulting from multiplying SGNFCAND
 189 (the significand) by 2 raised to the power of EXP (the exponent).   */)
 190   (Lisp_Object sgnfcand, Lisp_Object exponent)
 191 {
 192   CHECK_NUMBER (exponent);
 193   return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exponent)));
 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   double f1, f2, f3;
 210
 211   CHECK_NUMBER_OR_FLOAT (arg1);
 212   CHECK_NUMBER_OR_FLOAT (arg2);
 213   if (INTEGERP (arg1)     /* common lisp spec */
 214       && INTEGERP (arg2)   /* don't promote, if both are ints, and */
 215       && XINT (arg2) >= 0) /* we are sure the result is not fractional */
 216     {                           /* this can be improved by pre-calculating */
 217       EMACS_INT y;              /* some binary powers of x then accumulating */
 218       EMACS_UINT acc, x;  /* Unsigned so that overflow is well defined.  */
 219       Lisp_Object val;
 220
 221       x = XINT (arg1);
 222       y = XINT (arg2);
 223       acc = (y & 1 ? x : 1);
 224
 225       while ((y >>= 1) != 0)
 226         {
 227           x *= x;
 228           if (y & 1)
 229             acc *= x;
 230         }
 231       XSETINT (val, acc);
 232       return val;
 233     }
 234   f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
 235   f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
 236   f3 = pow (f1, f2);
 237   return make_float (f3);
 238 }
 239
 240 DEFUN ("log", Flog, Slog, 1, 2, 0,
 241        doc: /* Return the natural logarithm of ARG.
 242 If the optional argument BASE is given, return log ARG using that base.  */)
 243   (Lisp_Object arg, Lisp_Object base)
 244 {
 245   double d = extract_float (arg);
 246
 247   if (NILP (base))
 248     d = log (d);
 249   else
 250     {
 251       double b = extract_float (base);
 252
 253       if (b == 10.0)
 254         d = log10 (d);
 255       else
 256         d = log (d) / log (b);
 257     }
 258   return make_float (d);
 259 }
 260
 261 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
 262        doc: /* Return the logarithm base 10 of ARG.  */)
 263   (Lisp_Object arg)
 264 {
 265   double d = extract_float (arg);
 266   d = log10 (d);
 267   return make_float (d);
 268 }
 269
 270 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 271        doc: /* Return the square root of ARG.  */)
 272   (Lisp_Object arg)
 273 {
 274   double d = extract_float (arg);
 275   d = sqrt (d);
 276   return make_float (d);
 277 }
 278 \f
 279 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 280        doc: /* Return the absolute value of ARG.  */)
 281   (register Lisp_Object arg)
 282 {
 283   CHECK_NUMBER_OR_FLOAT (arg);
 284
 285   if (FLOATP (arg))
 286     arg = make_float (fabs (XFLOAT_DATA (arg)));
 287   else if (XINT (arg) < 0)
 288     XSETINT (arg, - XINT (arg));
 289
 290   return arg;
 291 }
 292
 293 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 294        doc: /* Return the floating point number equal to ARG.  */)
 295   (register Lisp_Object arg)
 296 {
 297   CHECK_NUMBER_OR_FLOAT (arg);
 298
 299   if (INTEGERP (arg))
 300     return make_float ((double) XINT (arg));
 301   else                          /* give 'em the same float back */
 302     return arg;
 303 }
 304
 305 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 306        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 307 This is the same as the exponent of a float.  */)
 308   (Lisp_Object arg)
 309 {
 310   Lisp_Object val;
 311   EMACS_INT value;
 312   double f = extract_float (arg);
 313
 314   if (f == 0.0)
 315     value = MOST_NEGATIVE_FIXNUM;
 316   else if (isfinite (f))
 317     {
 318       int ivalue;
 319       frexp (f, &ivalue);
 320       value = ivalue - 1;
 321     }
 322   else
 323     value = MOST_POSITIVE_FIXNUM;
 324
 325   XSETINT (val, value);
 326   return val;
 327 }
 328
 329
 330 /* the rounding functions  */
 331
 332 static Lisp_Object
 333 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 334                  double (*double_round) (double),
 335                  EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
 336                  const char *name)
 337 {
 338   CHECK_NUMBER_OR_FLOAT (arg);
 339
 340   if (! NILP (divisor))
 341     {
 342       EMACS_INT i1, i2;
 343
 344       CHECK_NUMBER_OR_FLOAT (divisor);
 345
 346       if (FLOATP (arg) || FLOATP (divisor))
 347         {
 348           double f1, f2;
 349
 350           f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 351           f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
 352           if (! IEEE_FLOATING_POINT && f2 == 0)
 353             xsignal0 (Qarith_error);
 354
 355           f1 = (*double_round) (f1 / f2);
 356           if (FIXNUM_OVERFLOW_P (f1))
 357             xsignal3 (Qrange_error, build_string (name), arg, divisor);
 358           arg = make_number (f1);
 359           return arg;
 360         }
 361
 362       i1 = XINT (arg);
 363       i2 = XINT (divisor);
 364
 365       if (i2 == 0)
 366         xsignal0 (Qarith_error);
 367
 368       XSETINT (arg, (*int_round2) (i1, i2));
 369       return arg;
 370     }
 371
 372   if (FLOATP (arg))
 373     {
 374       double d = (*double_round) (XFLOAT_DATA (arg));
 375       if (FIXNUM_OVERFLOW_P (d))
 376         xsignal2 (Qrange_error, build_string (name), arg);
 377       arg = make_number (d);
 378     }
 379
 380   return arg;
 381 }
 382
 383 /* With C's /, the result is implementation-defined if either operand
 384    is negative, so take care with negative operands in the following
 385    integer functions.  */
 386
 387 static EMACS_INT
 388 ceiling2 (EMACS_INT i1, EMACS_INT i2)
 389 {
 390   return (i2 < 0
 391           ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
 392           : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
 393 }
 394
 395 static EMACS_INT
 396 floor2 (EMACS_INT i1, EMACS_INT i2)
 397 {
 398   return (i2 < 0
 399           ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 400           : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 401 }
 402
 403 static EMACS_INT
 404 truncate2 (EMACS_INT i1, EMACS_INT i2)
 405 {
 406   return (i2 < 0
 407           ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
 408           : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
 409 }
 410
 411 static EMACS_INT
 412 round2 (EMACS_INT i1, EMACS_INT i2)
 413 {
 414   /* The C language's division operator gives us one remainder R, but
 415      we want the remainder R1 on the other side of 0 if R1 is closer
 416      to 0 than R is; because we want to round to even, we also want R1
 417      if R and R1 are the same distance from 0 and if C's quotient is
 418      odd.  */
 419   EMACS_INT q = i1 / i2;
 420   EMACS_INT r = i1 % i2;
 421   EMACS_INT abs_r = eabs (r);
 422   EMACS_INT abs_r1 = eabs (i2) - abs_r;
 423   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 424 }
 425
 426 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 427    if `rint' exists but does not work right.  */
 428 #ifdef HAVE_RINT
 429 #define emacs_rint rint
 430 #else
 431 static double
 432 emacs_rint (double d)
 433 {
 434   return floor (d + 0.5);
 435 }
 436 #endif
 437
 438 static double
 439 double_identity (double d)
 440 {
 441   return d;
 442 }
 443
 444 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 445        doc: /* Return the smallest integer no less than ARG.
 446 This rounds the value towards +inf.
 447 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 448   (Lisp_Object arg, Lisp_Object divisor)
 449 {
 450   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 451 }
 452
 453 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 454        doc: /* Return the largest integer no greater than ARG.
 455 This rounds the value towards -inf.
 456 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 457   (Lisp_Object arg, Lisp_Object divisor)
 458 {
 459   return rounding_driver (arg, divisor, floor, floor2, "floor");
 460 }
 461
 462 DEFUN ("round", Fround, Sround, 1, 2, 0,
 463        doc: /* Return the nearest integer to ARG.
 464 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 465
 466 Rounding a value equidistant between two integers may choose the
 467 integer closer to zero, or it may prefer an even integer, depending on
 468 your machine.  For example, \(round 2.5\) can return 3 on some
 469 systems, but 2 on others.  */)
 470   (Lisp_Object arg, Lisp_Object divisor)
 471 {
 472   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 473 }
 474
 475 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 476        doc: /* Truncate a floating point number to an int.
 477 Rounds ARG toward zero.
 478 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 479   (Lisp_Object arg, Lisp_Object divisor)
 480 {
 481   return rounding_driver (arg, divisor, double_identity, truncate2,
 482                           "truncate");
 483 }
 484
 485
 486 Lisp_Object
 487 fmod_float (Lisp_Object x, Lisp_Object y)
 488 {
 489   double f1, f2;
 490
 491   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 492   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 493
 494   f1 = fmod (f1, f2);
 495
 496   /* If the "remainder" comes out with the wrong sign, fix it.  */
 497   if (f2 < 0 ? f1 > 0 : f1 < 0)
 498     f1 += f2;
 499
 500   return make_float (f1);
 501 }
 502 \f
 503 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 504        doc: /* Return the smallest integer no less than ARG, as a float.
 505 \(Round toward +inf.\)  */)
 506   (Lisp_Object arg)
 507 {
 508   double d = extract_float (arg);
 509   d = ceil (d);
 510   return make_float (d);
 511 }
 512
 513 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 514        doc: /* Return the largest integer no greater than ARG, as a float.
 515 \(Round towards -inf.\)  */)
 516   (Lisp_Object arg)
 517 {
 518   double d = extract_float (arg);
 519   d = floor (d);
 520   return make_float (d);
 521 }
 522
 523 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 524        doc: /* Return the nearest integer to ARG, as a float.  */)
 525   (Lisp_Object arg)
 526 {
 527   double d = extract_float (arg);
 528   d = emacs_rint (d);
 529   return make_float (d);
 530 }
 531
 532 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 533        doc: /* Truncate a floating point number to an integral float value.
 534 Rounds the value toward zero.  */)
 535   (Lisp_Object arg)
 536 {
 537   double d = extract_float (arg);
 538   if (d >= 0.0)
 539     d = floor (d);
 540   else
 541     d = ceil (d);
 542   return make_float (d);
 543 }
 544 \f
 545 void
 546 syms_of_floatfns (void)
 547 {
 548   defsubr (&Sacos);
 549   defsubr (&Sasin);
 550   defsubr (&Satan);
 551   defsubr (&Scos);
 552   defsubr (&Ssin);
 553   defsubr (&Stan);
 554   defsubr (&Sisnan);
 555 #ifdef HAVE_COPYSIGN
 556   defsubr (&Scopysign);
 557 #endif
 558   defsubr (&Sfrexp);
 559   defsubr (&Sldexp);
 560   defsubr (&Sfceiling);
 561   defsubr (&Sffloor);
 562   defsubr (&Sfround);
 563   defsubr (&Sftruncate);
 564   defsubr (&Sexp);
 565   defsubr (&Sexpt);
 566   defsubr (&Slog);
 567   defsubr (&Slog10);
 568   defsubr (&Ssqrt);
 569
 570   defsubr (&Sabs);
 571   defsubr (&Sfloat);
 572   defsubr (&Slogb);
 573   defsubr (&Sceiling);
 574   defsubr (&Sfloor);
 575   defsubr (&Sround);
 576   defsubr (&Struncate);
 577 }