src/floatfns.c

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