src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2
   3 Copyright (C) 1988, 1993-1994, 1999, 2001-2011
   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 /* ANSI C requires only these float functions:
  26    acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
  27    frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
  28
  29    Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
  30    Define HAVE_CBRT if you have cbrt.
  31    Define HAVE_RINT if you have a working rint.
  32    If you don't define these, then the appropriate routines will be simulated.
  33
  34    Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
  35    (This should happen automatically.)
  36
  37    Define FLOAT_CHECK_ERRNO if the float library routines set errno.
  38    This has no effect if HAVE_MATHERR is defined.
  39
  40    Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
  41    (What systems actually do this?  Please let us know.)
  42
  43    Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
  44    either setting errno, or signaling SIGFPE/SIGILL.  Otherwise, domain and
  45    range checking will happen before calling the float routines.  This has
  46    no effect if HAVE_MATHERR is defined (since matherr will be called when
  47    a domain error occurs.)
  48  */
  49
  50 #include <config.h>
  51 #include <signal.h>
  52 #include <setjmp.h>
  53 #include "lisp.h"
  54 #include "syssignal.h"
  55
  56 #if STDC_HEADERS
  57 #include <float.h>
  58 #endif
  59
  60 /* If IEEE_FLOATING_POINT isn't defined, default it from FLT_*. */
  61 #ifndef IEEE_FLOATING_POINT
  62 #if (FLT_RADIX == 2 && FLT_MANT_DIG == 24 \
  63      && FLT_MIN_EXP == -125 && FLT_MAX_EXP == 128)
  64 #define IEEE_FLOATING_POINT 1
  65 #else
  66 #define IEEE_FLOATING_POINT 0
  67 #endif
  68 #endif
  69
  70 #include <math.h>
  71
  72 /* This declaration is omitted on some systems, like Ultrix.  */
  73 #if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
  74 extern double logb (double);
  75 #endif /* not HPUX and HAVE_LOGB and no logb macro */
  76
  77 #if defined(DOMAIN) && defined(SING) && defined(OVERFLOW)
  78     /* If those are defined, then this is probably a `matherr' machine. */
  79 # ifndef HAVE_MATHERR
  80 #  define HAVE_MATHERR
  81 # endif
  82 #endif
  83
  84 #ifdef NO_MATHERR
  85 #undef HAVE_MATHERR
  86 #endif
  87
  88 #ifdef HAVE_MATHERR
  89 # ifdef FLOAT_CHECK_ERRNO
  90 #  undef FLOAT_CHECK_ERRNO
  91 # endif
  92 # ifdef FLOAT_CHECK_DOMAIN
  93 #  undef FLOAT_CHECK_DOMAIN
  94 # endif
  95 #endif
  96
  97 #ifndef NO_FLOAT_CHECK_ERRNO
  98 #define FLOAT_CHECK_ERRNO
  99 #endif
 100
 101 #ifdef FLOAT_CHECK_ERRNO
 102 # include <errno.h>
 103 #endif
 104
 105 #ifdef FLOAT_CATCH_SIGILL
 106 static void float_error ();
 107 #endif
 108
 109 /* Nonzero while executing in floating point.
 110    This tells float_error what to do.  */
 111
 112 static int in_float;
 113
 114 /* If an argument is out of range for a mathematical function,
 115    here is the actual argument value to use in the error message.
 116    These variables are used only across the floating point library call
 117    so there is no need to staticpro them.  */
 118
 119 static Lisp_Object float_error_arg, float_error_arg2;
 120
 121 static const char *float_error_fn_name;
 122
 123 /* Evaluate the floating point expression D, recording NUM
 124    as the original argument for error messages.
 125    D is normally an assignment expression.
 126    Handle errors which may result in signals or may set errno.
 127
 128    Note that float_error may be declared to return void, so you can't
 129    just cast the zero after the colon to (void) to make the types
 130    check properly.  */
 131
 132 #ifdef FLOAT_CHECK_ERRNO
 133 #define IN_FLOAT(d, name, num)                          \
 134   do {                                                  \
 135     float_error_arg = num;                              \
 136     float_error_fn_name = name;                         \
 137     in_float = 1; errno = 0; (d); in_float = 0;         \
 138     switch (errno) {                                    \
 139     case 0: break;                                      \
 140     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 141     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 142     default:     arith_error (float_error_fn_name, float_error_arg);    \
 143     }                                                   \
 144   } while (0)
 145 #define IN_FLOAT2(d, name, num, num2)                   \
 146   do {                                                  \
 147     float_error_arg = num;                              \
 148     float_error_arg2 = num2;                            \
 149     float_error_fn_name = name;                         \
 150     in_float = 1; errno = 0; (d); in_float = 0;         \
 151     switch (errno) {                                    \
 152     case 0: break;                                      \
 153     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 154     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 155     default:     arith_error (float_error_fn_name, float_error_arg);    \
 156     }                                                   \
 157   } while (0)
 158 #else
 159 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
 160 #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
 161 #endif
 162
 163 /* Convert float to Lisp_Int if it fits, else signal a range error
 164    using the given arguments.  */
 165 #define FLOAT_TO_INT(x, i, name, num)                                   \
 166   do                                                                    \
 167     {                                                                   \
 168       if (FIXNUM_OVERFLOW_P (x))                                        \
 169         range_error (name, num);                                        \
 170       XSETINT (i,  (EMACS_INT)(x));                                     \
 171     }                                                                   \
 172   while (0)
 173 #define FLOAT_TO_INT2(x, i, name, num1, num2)                           \
 174   do                                                                    \
 175     {                                                                   \
 176       if (FIXNUM_OVERFLOW_P (x))                                        \
 177         range_error2 (name, num1, num2);                                \
 178       XSETINT (i,  (EMACS_INT)(x));                                     \
 179     }                                                                   \
 180   while (0)
 181
 182 #define arith_error(op,arg) \
 183   xsignal2 (Qarith_error, build_string ((op)), (arg))
 184 #define range_error(op,arg) \
 185   xsignal2 (Qrange_error, build_string ((op)), (arg))
 186 #define range_error2(op,a1,a2) \
 187   xsignal3 (Qrange_error, build_string ((op)), (a1), (a2))
 188 #define domain_error(op,arg) \
 189   xsignal2 (Qdomain_error, build_string ((op)), (arg))
 190 #ifdef FLOAT_CHECK_DOMAIN
 191 #define domain_error2(op,a1,a2) \
 192   xsignal3 (Qdomain_error, build_string ((op)), (a1), (a2))
 193 #endif
 194
 195 /* Extract a Lisp number as a `double', or signal an error.  */
 196
 197 double
 198 extract_float (Lisp_Object num)
 199 {
 200   CHECK_NUMBER_OR_FLOAT (num);
 201
 202   if (FLOATP (num))
 203     return XFLOAT_DATA (num);
 204   return (double) XINT (num);
 205 }
 206 \f
 207 /* Trig functions.  */
 208
 209 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
 210        doc: /* Return the inverse cosine of ARG.  */)
 211   (register Lisp_Object arg)
 212 {
 213   double d = extract_float (arg);
 214 #ifdef FLOAT_CHECK_DOMAIN
 215   if (d > 1.0 || d < -1.0)
 216     domain_error ("acos", arg);
 217 #endif
 218   IN_FLOAT (d = acos (d), "acos", arg);
 219   return make_float (d);
 220 }
 221
 222 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
 223        doc: /* Return the inverse sine of ARG.  */)
 224   (register Lisp_Object arg)
 225 {
 226   double d = extract_float (arg);
 227 #ifdef FLOAT_CHECK_DOMAIN
 228   if (d > 1.0 || d < -1.0)
 229     domain_error ("asin", arg);
 230 #endif
 231   IN_FLOAT (d = asin (d), "asin", arg);
 232   return make_float (d);
 233 }
 234
 235 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
 236        doc: /* Return the inverse tangent of the arguments.
 237 If only one argument Y is given, return the inverse tangent of Y.
 238 If two arguments Y and X are given, return the inverse tangent of Y
 239 divided by X, i.e. the angle in radians between the vector (X, Y)
 240 and the x-axis.  */)
 241   (register Lisp_Object y, Lisp_Object x)
 242 {
 243   double d = extract_float (y);
 244
 245   if (NILP (x))
 246     IN_FLOAT (d = atan (d), "atan", y);
 247   else
 248     {
 249       double d2 = extract_float (x);
 250
 251       IN_FLOAT2 (d = atan2 (d, d2), "atan", y, x);
 252     }
 253   return make_float (d);
 254 }
 255
 256 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 257        doc: /* Return the cosine of ARG.  */)
 258   (register Lisp_Object arg)
 259 {
 260   double d = extract_float (arg);
 261   IN_FLOAT (d = cos (d), "cos", arg);
 262   return make_float (d);
 263 }
 264
 265 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 266        doc: /* Return the sine of ARG.  */)
 267   (register Lisp_Object arg)
 268 {
 269   double d = extract_float (arg);
 270   IN_FLOAT (d = sin (d), "sin", arg);
 271   return make_float (d);
 272 }
 273
 274 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 275        doc: /* Return the tangent of ARG.  */)
 276   (register Lisp_Object arg)
 277 {
 278   double d = extract_float (arg);
 279   double c = cos (d);
 280 #ifdef FLOAT_CHECK_DOMAIN
 281   if (c == 0.0)
 282     domain_error ("tan", arg);
 283 #endif
 284   IN_FLOAT (d = sin (d) / c, "tan", arg);
 285   return make_float (d);
 286 }
 287
 288 #if defined HAVE_ISNAN && defined HAVE_COPYSIGN
 289 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
 290        doc: /* Return non nil iff argument X is a NaN.  */)
 291   (Lisp_Object x)
 292 {
 293   CHECK_FLOAT (x);
 294   return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
 295 }
 296
 297 DEFUN ("copysign", Fcopysign, Scopysign, 1, 2, 0,
 298        doc: /* Copy sign of X2 to value of X1, and return the result.
 299 Cause an error if X1 or X2 is not a float.  */)
 300   (Lisp_Object x1, Lisp_Object x2)
 301 {
 302   double f1, f2;
 303
 304   CHECK_FLOAT (x1);
 305   CHECK_FLOAT (x2);
 306
 307   f1 = XFLOAT_DATA (x1);
 308   f2 = XFLOAT_DATA (x2);
 309
 310   return make_float (copysign (f1, f2));
 311 }
 312
 313 DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
 314        doc: /* Get significand and exponent of a floating point number.
 315 Breaks the floating point number X into its binary significand SGNFCAND
 316 \(a floating point value between 0.5 (included) and 1.0 (excluded))
 317 and an integral exponent EXP for 2, such that:
 318
 319   X = SGNFCAND * 2^EXP
 320
 321 The function returns the cons cell (SGNFCAND . EXP).
 322 If X is zero, both parts (SGNFCAND and EXP) are zero.  */)
 323   (Lisp_Object x)
 324 {
 325   double f = XFLOATINT (x);
 326
 327   if (f == 0.0)
 328     return Fcons (make_float (0.0), make_number (0));
 329   else
 330     {
 331       int exponent;
 332       double sgnfcand = frexp (f, &exponent);
 333       return Fcons (make_float (sgnfcand), make_number (exponent));
 334     }
 335 }
 336
 337 DEFUN ("ldexp", Fldexp, Sldexp, 1, 2, 0,
 338        doc: /* Construct number X from significand SGNFCAND and exponent EXP.
 339 Returns the floating point value resulting from multiplying SGNFCAND
 340 (the significand) by 2 raised to the power of EXP (the exponent).   */)
 341   (Lisp_Object sgnfcand, Lisp_Object exponent)
 342 {
 343   CHECK_NUMBER (exponent);
 344   return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exponent)));
 345 }
 346 #endif
 347 \f
 348 #if 0 /* Leave these out unless we find there's a reason for them.  */
 349
 350 DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
 351        doc: /* Return the bessel function j0 of ARG.  */)
 352   (register Lisp_Object arg)
 353 {
 354   double d = extract_float (arg);
 355   IN_FLOAT (d = j0 (d), "bessel-j0", arg);
 356   return make_float (d);
 357 }
 358
 359 DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
 360        doc: /* Return the bessel function j1 of ARG.  */)
 361   (register Lisp_Object arg)
 362 {
 363   double d = extract_float (arg);
 364   IN_FLOAT (d = j1 (d), "bessel-j1", arg);
 365   return make_float (d);
 366 }
 367
 368 DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
 369        doc: /* Return the order N bessel function output jn of ARG.
 370 The first arg (the order) is truncated to an integer.  */)
 371   (register Lisp_Object n, Lisp_Object arg)
 372 {
 373   int i1 = extract_float (n);
 374   double f2 = extract_float (arg);
 375
 376   IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
 377   return make_float (f2);
 378 }
 379
 380 DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
 381        doc: /* Return the bessel function y0 of ARG.  */)
 382   (register Lisp_Object arg)
 383 {
 384   double d = extract_float (arg);
 385   IN_FLOAT (d = y0 (d), "bessel-y0", arg);
 386   return make_float (d);
 387 }
 388
 389 DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
 390        doc: /* Return the bessel function y1 of ARG.  */)
 391   (register Lisp_Object arg)
 392 {
 393   double d = extract_float (arg);
 394   IN_FLOAT (d = y1 (d), "bessel-y0", arg);
 395   return make_float (d);
 396 }
 397
 398 DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
 399        doc: /* Return the order N bessel function output yn of ARG.
 400 The first arg (the order) is truncated to an integer.  */)
 401   (register Lisp_Object n, Lisp_Object arg)
 402 {
 403   int i1 = extract_float (n);
 404   double f2 = extract_float (arg);
 405
 406   IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
 407   return make_float (f2);
 408 }
 409
 410 #endif
 411 \f
 412 #if 0 /* Leave these out unless we see they are worth having.  */
 413
 414 DEFUN ("erf", Ferf, Serf, 1, 1, 0,
 415        doc: /* Return the mathematical error function of ARG.  */)
 416   (register Lisp_Object arg)
 417 {
 418   double d = extract_float (arg);
 419   IN_FLOAT (d = erf (d), "erf", arg);
 420   return make_float (d);
 421 }
 422
 423 DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
 424        doc: /* Return the complementary error function of ARG.  */)
 425   (register Lisp_Object arg)
 426 {
 427   double d = extract_float (arg);
 428   IN_FLOAT (d = erfc (d), "erfc", arg);
 429   return make_float (d);
 430 }
 431
 432 DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
 433        doc: /* Return the log gamma of ARG.  */)
 434   (register Lisp_Object arg)
 435 {
 436   double d = extract_float (arg);
 437   IN_FLOAT (d = lgamma (d), "log-gamma", arg);
 438   return make_float (d);
 439 }
 440
 441 DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
 442        doc: /* Return the cube root of ARG.  */)
 443   (register Lisp_Object arg)
 444 {
 445   double d = extract_float (arg);
 446 #ifdef HAVE_CBRT
 447   IN_FLOAT (d = cbrt (d), "cube-root", arg);
 448 #else
 449   if (d >= 0.0)
 450     IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
 451   else
 452     IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
 453 #endif
 454   return make_float (d);
 455 }
 456
 457 #endif
 458 \f
 459 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 460        doc: /* Return the exponential base e of ARG.  */)
 461   (register Lisp_Object arg)
 462 {
 463   double d = extract_float (arg);
 464 #ifdef FLOAT_CHECK_DOMAIN
 465   if (d > 709.7827)   /* Assume IEEE doubles here */
 466     range_error ("exp", arg);
 467   else if (d < -709.0)
 468     return make_float (0.0);
 469   else
 470 #endif
 471     IN_FLOAT (d = exp (d), "exp", arg);
 472   return make_float (d);
 473 }
 474
 475 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 476        doc: /* Return the exponential ARG1 ** ARG2.  */)
 477   (register Lisp_Object arg1, Lisp_Object arg2)
 478 {
 479   double f1, f2, f3;
 480
 481   CHECK_NUMBER_OR_FLOAT (arg1);
 482   CHECK_NUMBER_OR_FLOAT (arg2);
 483   if (INTEGERP (arg1)     /* common lisp spec */
 484       && INTEGERP (arg2)   /* don't promote, if both are ints, and */
 485       && 0 <= XINT (arg2)) /* we are sure the result is not fractional */
 486     {                           /* this can be improved by pre-calculating */
 487       EMACS_INT acc, x, y;      /* some binary powers of x then accumulating */
 488       Lisp_Object val;
 489
 490       x = XINT (arg1);
 491       y = XINT (arg2);
 492       acc = 1;
 493
 494       if (y < 0)
 495         {
 496           if (x == 1)
 497             acc = 1;
 498           else if (x == -1)
 499             acc = (y & 1) ? -1 : 1;
 500           else
 501             acc = 0;
 502         }
 503       else
 504         {
 505           while (y > 0)
 506             {
 507               if (y & 1)
 508                 acc *= x;
 509               x *= x;
 510               y = (unsigned)y >> 1;
 511             }
 512         }
 513       XSETINT (val, acc);
 514       return val;
 515     }
 516   f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
 517   f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
 518   /* Really should check for overflow, too */
 519   if (f1 == 0.0 && f2 == 0.0)
 520     f1 = 1.0;
 521 #ifdef FLOAT_CHECK_DOMAIN
 522   else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor(f2)))
 523     domain_error2 ("expt", arg1, arg2);
 524 #endif
 525   IN_FLOAT2 (f3 = pow (f1, f2), "expt", arg1, arg2);
 526   /* Check for overflow in the result.  */
 527   if (f1 != 0.0 && f3 == 0.0)
 528     range_error ("expt", arg1);
 529   return make_float (f3);
 530 }
 531
 532 DEFUN ("log", Flog, Slog, 1, 2, 0,
 533        doc: /* Return the natural logarithm of ARG.
 534 If the optional argument BASE is given, return log ARG using that base.  */)
 535   (register Lisp_Object arg, Lisp_Object base)
 536 {
 537   double d = extract_float (arg);
 538
 539 #ifdef FLOAT_CHECK_DOMAIN
 540   if (d <= 0.0)
 541     domain_error2 ("log", arg, base);
 542 #endif
 543   if (NILP (base))
 544     IN_FLOAT (d = log (d), "log", arg);
 545   else
 546     {
 547       double b = extract_float (base);
 548
 549 #ifdef FLOAT_CHECK_DOMAIN
 550       if (b <= 0.0 || b == 1.0)
 551         domain_error2 ("log", arg, base);
 552 #endif
 553       if (b == 10.0)
 554         IN_FLOAT2 (d = log10 (d), "log", arg, base);
 555       else
 556         IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
 557     }
 558   return make_float (d);
 559 }
 560
 561 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
 562        doc: /* Return the logarithm base 10 of ARG.  */)
 563   (register Lisp_Object arg)
 564 {
 565   double d = extract_float (arg);
 566 #ifdef FLOAT_CHECK_DOMAIN
 567   if (d <= 0.0)
 568     domain_error ("log10", arg);
 569 #endif
 570   IN_FLOAT (d = log10 (d), "log10", arg);
 571   return make_float (d);
 572 }
 573
 574 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 575        doc: /* Return the square root of ARG.  */)
 576   (register Lisp_Object arg)
 577 {
 578   double d = extract_float (arg);
 579 #ifdef FLOAT_CHECK_DOMAIN
 580   if (d < 0.0)
 581     domain_error ("sqrt", arg);
 582 #endif
 583   IN_FLOAT (d = sqrt (d), "sqrt", arg);
 584   return make_float (d);
 585 }
 586 \f
 587 #if 0 /* Not clearly worth adding.  */
 588
 589 DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
 590        doc: /* Return the inverse hyperbolic cosine of ARG.  */)
 591   (register Lisp_Object arg)
 592 {
 593   double d = extract_float (arg);
 594 #ifdef FLOAT_CHECK_DOMAIN
 595   if (d < 1.0)
 596     domain_error ("acosh", arg);
 597 #endif
 598 #ifdef HAVE_INVERSE_HYPERBOLIC
 599   IN_FLOAT (d = acosh (d), "acosh", arg);
 600 #else
 601   IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
 602 #endif
 603   return make_float (d);
 604 }
 605
 606 DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
 607        doc: /* Return the inverse hyperbolic sine of ARG.  */)
 608   (register Lisp_Object arg)
 609 {
 610   double d = extract_float (arg);
 611 #ifdef HAVE_INVERSE_HYPERBOLIC
 612   IN_FLOAT (d = asinh (d), "asinh", arg);
 613 #else
 614   IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
 615 #endif
 616   return make_float (d);
 617 }
 618
 619 DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
 620        doc: /* Return the inverse hyperbolic tangent of ARG.  */)
 621   (register Lisp_Object arg)
 622 {
 623   double d = extract_float (arg);
 624 #ifdef FLOAT_CHECK_DOMAIN
 625   if (d >= 1.0 || d <= -1.0)
 626     domain_error ("atanh", arg);
 627 #endif
 628 #ifdef HAVE_INVERSE_HYPERBOLIC
 629   IN_FLOAT (d = atanh (d), "atanh", arg);
 630 #else
 631   IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
 632 #endif
 633   return make_float (d);
 634 }
 635
 636 DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
 637        doc: /* Return the hyperbolic cosine of ARG.  */)
 638   (register Lisp_Object arg)
 639 {
 640   double d = extract_float (arg);
 641 #ifdef FLOAT_CHECK_DOMAIN
 642   if (d > 710.0 || d < -710.0)
 643     range_error ("cosh", arg);
 644 #endif
 645   IN_FLOAT (d = cosh (d), "cosh", arg);
 646   return make_float (d);
 647 }
 648
 649 DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
 650        doc: /* Return the hyperbolic sine of ARG.  */)
 651   (register Lisp_Object arg)
 652 {
 653   double d = extract_float (arg);
 654 #ifdef FLOAT_CHECK_DOMAIN
 655   if (d > 710.0 || d < -710.0)
 656     range_error ("sinh", arg);
 657 #endif
 658   IN_FLOAT (d = sinh (d), "sinh", arg);
 659   return make_float (d);
 660 }
 661
 662 DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
 663        doc: /* Return the hyperbolic tangent of ARG.  */)
 664   (register Lisp_Object arg)
 665 {
 666   double d = extract_float (arg);
 667   IN_FLOAT (d = tanh (d), "tanh", arg);
 668   return make_float (d);
 669 }
 670 #endif
 671 \f
 672 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 673        doc: /* Return the absolute value of ARG.  */)
 674   (register Lisp_Object arg)
 675 {
 676   CHECK_NUMBER_OR_FLOAT (arg);
 677
 678   if (FLOATP (arg))
 679     IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
 680   else if (XINT (arg) < 0)
 681     XSETINT (arg, - XINT (arg));
 682
 683   return arg;
 684 }
 685
 686 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 687        doc: /* Return the floating point number equal to ARG.  */)
 688   (register Lisp_Object arg)
 689 {
 690   CHECK_NUMBER_OR_FLOAT (arg);
 691
 692   if (INTEGERP (arg))
 693     return make_float ((double) XINT (arg));
 694   else                          /* give 'em the same float back */
 695     return arg;
 696 }
 697
 698 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 699        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 700 This is the same as the exponent of a float.  */)
 701   (Lisp_Object arg)
 702 {
 703   Lisp_Object val;
 704   EMACS_INT value;
 705   double f = extract_float (arg);
 706
 707   if (f == 0.0)
 708     value = MOST_NEGATIVE_FIXNUM;
 709   else
 710     {
 711 #ifdef HAVE_LOGB
 712       IN_FLOAT (value = logb (f), "logb", arg);
 713 #else
 714 #ifdef HAVE_FREXP
 715       int ivalue;
 716       IN_FLOAT (frexp (f, &ivalue), "logb", arg);
 717       value = ivalue - 1;
 718 #else
 719       int i;
 720       double d;
 721       if (f < 0.0)
 722         f = -f;
 723       value = -1;
 724       while (f < 0.5)
 725         {
 726           for (i = 1, d = 0.5; d * d >= f; i += i)
 727             d *= d;
 728           f /= d;
 729           value -= i;
 730         }
 731       while (f >= 1.0)
 732         {
 733           for (i = 1, d = 2.0; d * d <= f; i += i)
 734             d *= d;
 735           f /= d;
 736           value += i;
 737         }
 738 #endif
 739 #endif
 740     }
 741   XSETINT (val, value);
 742   return val;
 743 }
 744
 745
 746 /* the rounding functions  */
 747
 748 static Lisp_Object
 749 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 750                  double (*double_round) (double),
 751                  EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
 752                  const char *name)
 753 {
 754   CHECK_NUMBER_OR_FLOAT (arg);
 755
 756   if (! NILP (divisor))
 757     {
 758       EMACS_INT i1, i2;
 759
 760       CHECK_NUMBER_OR_FLOAT (divisor);
 761
 762       if (FLOATP (arg) || FLOATP (divisor))
 763         {
 764           double f1, f2;
 765
 766           f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 767           f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
 768           if (! IEEE_FLOATING_POINT && f2 == 0)
 769             xsignal0 (Qarith_error);
 770
 771           IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
 772           FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
 773           return arg;
 774         }
 775
 776       i1 = XINT (arg);
 777       i2 = XINT (divisor);
 778
 779       if (i2 == 0)
 780         xsignal0 (Qarith_error);
 781
 782       XSETINT (arg, (*int_round2) (i1, i2));
 783       return arg;
 784     }
 785
 786   if (FLOATP (arg))
 787     {
 788       double d;
 789
 790       IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
 791       FLOAT_TO_INT (d, arg, name, arg);
 792     }
 793
 794   return arg;
 795 }
 796
 797 /* With C's /, the result is implementation-defined if either operand
 798    is negative, so take care with negative operands in the following
 799    integer functions.  */
 800
 801 static EMACS_INT
 802 ceiling2 (EMACS_INT i1, EMACS_INT i2)
 803 {
 804   return (i2 < 0
 805           ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
 806           : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
 807 }
 808
 809 static EMACS_INT
 810 floor2 (EMACS_INT i1, EMACS_INT i2)
 811 {
 812   return (i2 < 0
 813           ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 814           : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 815 }
 816
 817 static EMACS_INT
 818 truncate2 (EMACS_INT i1, EMACS_INT i2)
 819 {
 820   return (i2 < 0
 821           ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
 822           : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
 823 }
 824
 825 static EMACS_INT
 826 round2 (EMACS_INT i1, EMACS_INT i2)
 827 {
 828   /* The C language's division operator gives us one remainder R, but
 829      we want the remainder R1 on the other side of 0 if R1 is closer
 830      to 0 than R is; because we want to round to even, we also want R1
 831      if R and R1 are the same distance from 0 and if C's quotient is
 832      odd.  */
 833   EMACS_INT q = i1 / i2;
 834   EMACS_INT r = i1 % i2;
 835   EMACS_INT abs_r = r < 0 ? -r : r;
 836   EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
 837   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 838 }
 839
 840 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 841    if `rint' exists but does not work right.  */
 842 #ifdef HAVE_RINT
 843 #define emacs_rint rint
 844 #else
 845 static double
 846 emacs_rint (double d)
 847 {
 848   return floor (d + 0.5);
 849 }
 850 #endif
 851
 852 static double
 853 double_identity (double d)
 854 {
 855   return d;
 856 }
 857
 858 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 859        doc: /* Return the smallest integer no less than ARG.
 860 This rounds the value towards +inf.
 861 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 862   (Lisp_Object arg, Lisp_Object divisor)
 863 {
 864   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 865 }
 866
 867 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 868        doc: /* Return the largest integer no greater than ARG.
 869 This rounds the value towards -inf.
 870 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 871   (Lisp_Object arg, Lisp_Object divisor)
 872 {
 873   return rounding_driver (arg, divisor, floor, floor2, "floor");
 874 }
 875
 876 DEFUN ("round", Fround, Sround, 1, 2, 0,
 877        doc: /* Return the nearest integer to ARG.
 878 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 879
 880 Rounding a value equidistant between two integers may choose the
 881 integer closer to zero, or it may prefer an even integer, depending on
 882 your machine.  For example, \(round 2.5\) can return 3 on some
 883 systems, but 2 on others.  */)
 884   (Lisp_Object arg, Lisp_Object divisor)
 885 {
 886   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 887 }
 888
 889 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 890        doc: /* Truncate a floating point number to an int.
 891 Rounds ARG toward zero.
 892 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 893   (Lisp_Object arg, Lisp_Object divisor)
 894 {
 895   return rounding_driver (arg, divisor, double_identity, truncate2,
 896                           "truncate");
 897 }
 898
 899
 900 Lisp_Object
 901 fmod_float (Lisp_Object x, Lisp_Object y)
 902 {
 903   double f1, f2;
 904
 905   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 906   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 907
 908   if (! IEEE_FLOATING_POINT && f2 == 0)
 909     xsignal0 (Qarith_error);
 910
 911   /* If the "remainder" comes out with the wrong sign, fix it.  */
 912   IN_FLOAT2 ((f1 = fmod (f1, f2),
 913               f1 = (f2 < 0 ? f1 > 0 : f1 < 0) ? f1 + f2 : f1),
 914              "mod", x, y);
 915   return make_float (f1);
 916 }
 917 \f
 918 /* It's not clear these are worth adding.  */
 919
 920 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 921        doc: /* Return the smallest integer no less than ARG, as a float.
 922 \(Round toward +inf.\)  */)
 923   (register Lisp_Object arg)
 924 {
 925   double d = extract_float (arg);
 926   IN_FLOAT (d = ceil (d), "fceiling", arg);
 927   return make_float (d);
 928 }
 929
 930 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 931        doc: /* Return the largest integer no greater than ARG, as a float.
 932 \(Round towards -inf.\)  */)
 933   (register Lisp_Object arg)
 934 {
 935   double d = extract_float (arg);
 936   IN_FLOAT (d = floor (d), "ffloor", arg);
 937   return make_float (d);
 938 }
 939
 940 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 941        doc: /* Return the nearest integer to ARG, as a float.  */)
 942   (register Lisp_Object arg)
 943 {
 944   double d = extract_float (arg);
 945   IN_FLOAT (d = emacs_rint (d), "fround", arg);
 946   return make_float (d);
 947 }
 948
 949 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 950        doc: /* Truncate a floating point number to an integral float value.
 951 Rounds the value toward zero.  */)
 952   (register Lisp_Object arg)
 953 {
 954   double d = extract_float (arg);
 955   if (d >= 0.0)
 956     IN_FLOAT (d = floor (d), "ftruncate", arg);
 957   else
 958     IN_FLOAT (d = ceil (d), "ftruncate", arg);
 959   return make_float (d);
 960 }
 961 \f
 962 #ifdef FLOAT_CATCH_SIGILL
 963 static void
 964 float_error (signo)
 965      int signo;
 966 {
 967   if (! in_float)
 968     fatal_error_signal (signo);
 969
 970 #ifdef BSD_SYSTEM
 971   sigsetmask (SIGEMPTYMASK);
 972 #else
 973   /* Must reestablish handler each time it is called.  */
 974   signal (SIGILL, float_error);
 975 #endif /* BSD_SYSTEM */
 976
 977   SIGNAL_THREAD_CHECK (signo);
 978   in_float = 0;
 979
 980   xsignal1 (Qarith_error, float_error_arg);
 981 }
 982
 983 /* Another idea was to replace the library function `infnan'
 984    where SIGILL is signaled.  */
 985
 986 #endif /* FLOAT_CATCH_SIGILL */
 987
 988 #ifdef HAVE_MATHERR
 989 int
 990 matherr (struct exception *x)
 991 {
 992   Lisp_Object args;
 993   const char *name = x->name;
 994
 995   if (! in_float)
 996     /* Not called from emacs-lisp float routines; do the default thing. */
 997     return 0;
 998   if (!strcmp (x->name, "pow"))
 999     name = "expt";
1000
1001   args
1002     = Fcons (build_string (name),
1003              Fcons (make_float (x->arg1),
1004                     ((!strcmp (name, "log") || !strcmp (name, "pow"))
1005                      ? Fcons (make_float (x->arg2), Qnil)
1006                      : Qnil)));
1007   switch (x->type)
1008     {
1009     case DOMAIN:        xsignal (Qdomain_error, args);          break;
1010     case SING:          xsignal (Qsingularity_error, args);     break;
1011     case OVERFLOW:      xsignal (Qoverflow_error, args);        break;
1012     case UNDERFLOW:     xsignal (Qunderflow_error, args);       break;
1013     default:            xsignal (Qarith_error, args);           break;
1014     }
1015   return (1);   /* don't set errno or print a message */
1016 }
1017 #endif /* HAVE_MATHERR */
1018
1019 void
1020 init_floatfns (void)
1021 {
1022 #ifdef FLOAT_CATCH_SIGILL
1023   signal (SIGILL, float_error);
1024 #endif
1025   in_float = 0;
1026 }
1027
1028 void
1029 syms_of_floatfns (void)
1030 {
1031   defsubr (&Sacos);
1032   defsubr (&Sasin);
1033   defsubr (&Satan);
1034   defsubr (&Scos);
1035   defsubr (&Ssin);
1036   defsubr (&Stan);
1037 #if defined HAVE_ISNAN && defined HAVE_COPYSIGN
1038   defsubr (&Sisnan);
1039   defsubr (&Scopysign);
1040   defsubr (&Sfrexp);
1041   defsubr (&Sldexp);
1042 #endif
1043 #if 0
1044   defsubr (&Sacosh);
1045   defsubr (&Sasinh);
1046   defsubr (&Satanh);
1047   defsubr (&Scosh);
1048   defsubr (&Ssinh);
1049   defsubr (&Stanh);
1050   defsubr (&Sbessel_y0);
1051   defsubr (&Sbessel_y1);
1052   defsubr (&Sbessel_yn);
1053   defsubr (&Sbessel_j0);
1054   defsubr (&Sbessel_j1);
1055   defsubr (&Sbessel_jn);
1056   defsubr (&Serf);
1057   defsubr (&Serfc);
1058   defsubr (&Slog_gamma);
1059   defsubr (&Scube_root);
1060 #endif
1061   defsubr (&Sfceiling);
1062   defsubr (&Sffloor);
1063   defsubr (&Sfround);
1064   defsubr (&Sftruncate);
1065   defsubr (&Sexp);
1066   defsubr (&Sexpt);
1067   defsubr (&Slog);
1068   defsubr (&Slog10);
1069   defsubr (&Ssqrt);
1070
1071   defsubr (&Sabs);
1072   defsubr (&Sfloat);
1073   defsubr (&Slogb);
1074   defsubr (&Sceiling);
1075   defsubr (&Sfloor);
1076   defsubr (&Sround);
1077   defsubr (&Struncate);
1078 }