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