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