src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2
   3 Copyright (C) 1988, 1993-1994, 1999, 2001-2012
   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 #include <float.h>
  57 /* If IEEE_FLOATING_POINT isn't defined, default it from FLT_*. */
  58 #ifndef IEEE_FLOATING_POINT
  59 #if (FLT_RADIX == 2 && FLT_MANT_DIG == 24 \
  60      && FLT_MIN_EXP == -125 && FLT_MAX_EXP == 128)
  61 #define IEEE_FLOATING_POINT 1
  62 #else
  63 #define IEEE_FLOATING_POINT 0
  64 #endif
  65 #endif
  66
  67 #include <math.h>
  68
  69 /* This declaration is omitted on some systems, like Ultrix.  */
  70 #if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
  71 extern double logb (double);
  72 #endif /* not HPUX and HAVE_LOGB and no logb macro */
  73
  74 #if defined (DOMAIN) && defined (SING) && defined (OVERFLOW)
  75     /* If those are defined, then this is probably a `matherr' machine. */
  76 # ifndef HAVE_MATHERR
  77 #  define HAVE_MATHERR
  78 # endif
  79 #endif
  80
  81 #ifdef NO_MATHERR
  82 #undef HAVE_MATHERR
  83 #endif
  84
  85 #ifdef HAVE_MATHERR
  86 # ifdef FLOAT_CHECK_ERRNO
  87 #  undef FLOAT_CHECK_ERRNO
  88 # endif
  89 # ifdef FLOAT_CHECK_DOMAIN
  90 #  undef FLOAT_CHECK_DOMAIN
  91 # endif
  92 #endif
  93
  94 #ifndef NO_FLOAT_CHECK_ERRNO
  95 #define FLOAT_CHECK_ERRNO
  96 #endif
  97
  98 #ifdef FLOAT_CHECK_ERRNO
  99 # include <errno.h>
 100 #endif
 101
 102 #ifdef FLOAT_CATCH_SIGILL
 103 static void float_error ();
 104 #endif
 105
 106 /* Nonzero while executing in floating point.
 107    This tells float_error what to do.  */
 108
 109 static int in_float;
 110
 111 /* If an argument is out of range for a mathematical function,
 112    here is the actual argument value to use in the error message.
 113    These variables are used only across the floating point library call
 114    so there is no need to staticpro them.  */
 115
 116 static Lisp_Object float_error_arg, float_error_arg2;
 117
 118 static const char *float_error_fn_name;
 119
 120 /* Evaluate the floating point expression D, recording NUM
 121    as the original argument for error messages.
 122    D is normally an assignment expression.
 123    Handle errors which may result in signals or may set errno.
 124
 125    Note that float_error may be declared to return void, so you can't
 126    just cast the zero after the colon to (void) to make the types
 127    check properly.  */
 128
 129 #ifdef FLOAT_CHECK_ERRNO
 130 #define IN_FLOAT(d, name, num)                          \
 131   do {                                                  \
 132     float_error_arg = num;                              \
 133     float_error_fn_name = name;                         \
 134     in_float = 1; errno = 0; (d); in_float = 0;         \
 135     switch (errno) {                                    \
 136     case 0: break;                                      \
 137     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 138     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 139     default:     arith_error (float_error_fn_name, float_error_arg);    \
 140     }                                                   \
 141   } while (0)
 142 #define IN_FLOAT2(d, name, num, num2)                   \
 143   do {                                                  \
 144     float_error_arg = num;                              \
 145     float_error_arg2 = num2;                            \
 146     float_error_fn_name = name;                         \
 147     in_float = 1; errno = 0; (d); in_float = 0;         \
 148     switch (errno) {                                    \
 149     case 0: break;                                      \
 150     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 151     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 152     default:     arith_error (float_error_fn_name, float_error_arg);    \
 153     }                                                   \
 154   } while (0)
 155 #else
 156 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
 157 #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
 158 #endif
 159
 160 /* Convert float to Lisp_Int if it fits, else signal a range error
 161    using the given arguments.  */
 162 #define FLOAT_TO_INT(x, i, name, num)                                   \
 163   do                                                                    \
 164     {                                                                   \
 165       if (FIXNUM_OVERFLOW_P (x))                                        \
 166         range_error (name, num);                                        \
 167       XSETINT (i,  (EMACS_INT)(x));                                     \
 168     }                                                                   \
 169   while (0)
 170 #define FLOAT_TO_INT2(x, i, name, num1, num2)                           \
 171   do                                                                    \
 172     {                                                                   \
 173       if (FIXNUM_OVERFLOW_P (x))                                        \
 174         range_error2 (name, num1, num2);                                \
 175       XSETINT (i,  (EMACS_INT)(x));                                     \
 176     }                                                                   \
 177   while (0)
 178
 179 #define arith_error(op,arg) \
 180   xsignal2 (Qarith_error, build_string ((op)), (arg))
 181 #define range_error(op,arg) \
 182   xsignal2 (Qrange_error, build_string ((op)), (arg))
 183 #define range_error2(op,a1,a2) \
 184   xsignal3 (Qrange_error, build_string ((op)), (a1), (a2))
 185 #define domain_error(op,arg) \
 186   xsignal2 (Qdomain_error, build_string ((op)), (arg))
 187 #ifdef FLOAT_CHECK_DOMAIN
 188 #define domain_error2(op,a1,a2) \
 189   xsignal3 (Qdomain_error, build_string ((op)), (a1), (a2))
 190 #endif
 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 #undef isnan
 286 #define isnan(x) ((x) != (x))
 287
 288 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
 289        doc: /* Return non nil iff argument X is a NaN.  */)
 290   (Lisp_Object x)
 291 {
 292   CHECK_FLOAT (x);
 293   return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
 294 }
 295
 296 #ifdef HAVE_COPYSIGN
 297 DEFUN ("copysign", Fcopysign, Scopysign, 2, 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 y;              /* some binary powers of x then accumulating */
 488       EMACS_UINT acc, x;  /* Unsigned so that overflow is well defined.  */
 489       Lisp_Object val;
 490
 491       x = XINT (arg1);
 492       y = XINT (arg2);
 493       acc = (y & 1 ? x : 1);
 494
 495       while ((y >>= 1) != 0)
 496         {
 497           x *= x;
 498           if (y & 1)
 499             acc *= x;
 500         }
 501       XSETINT (val, acc);
 502       return val;
 503     }
 504   f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
 505   f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
 506   /* Really should check for overflow, too */
 507   if (f1 == 0.0 && f2 == 0.0)
 508     f1 = 1.0;
 509 #ifdef FLOAT_CHECK_DOMAIN
 510   else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor (f2)))
 511     domain_error2 ("expt", arg1, arg2);
 512 #endif
 513   IN_FLOAT2 (f3 = pow (f1, f2), "expt", arg1, arg2);
 514   /* Check for overflow in the result.  */
 515   if (f1 != 0.0 && f3 == 0.0)
 516     range_error ("expt", arg1);
 517   return make_float (f3);
 518 }
 519
 520 DEFUN ("log", Flog, Slog, 1, 2, 0,
 521        doc: /* Return the natural logarithm of ARG.
 522 If the optional argument BASE is given, return log ARG using that base.  */)
 523   (register Lisp_Object arg, Lisp_Object base)
 524 {
 525   double d = extract_float (arg);
 526
 527 #ifdef FLOAT_CHECK_DOMAIN
 528   if (d <= 0.0)
 529     domain_error2 ("log", arg, base);
 530 #endif
 531   if (NILP (base))
 532     IN_FLOAT (d = log (d), "log", arg);
 533   else
 534     {
 535       double b = extract_float (base);
 536
 537 #ifdef FLOAT_CHECK_DOMAIN
 538       if (b <= 0.0 || b == 1.0)
 539         domain_error2 ("log", arg, base);
 540 #endif
 541       if (b == 10.0)
 542         IN_FLOAT2 (d = log10 (d), "log", arg, base);
 543       else
 544         IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
 545     }
 546   return make_float (d);
 547 }
 548
 549 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
 550        doc: /* Return the logarithm base 10 of ARG.  */)
 551   (register Lisp_Object arg)
 552 {
 553   double d = extract_float (arg);
 554 #ifdef FLOAT_CHECK_DOMAIN
 555   if (d <= 0.0)
 556     domain_error ("log10", arg);
 557 #endif
 558   IN_FLOAT (d = log10 (d), "log10", arg);
 559   return make_float (d);
 560 }
 561
 562 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 563        doc: /* Return the square root of ARG.  */)
 564   (register Lisp_Object arg)
 565 {
 566   double d = extract_float (arg);
 567 #ifdef FLOAT_CHECK_DOMAIN
 568   if (d < 0.0)
 569     domain_error ("sqrt", arg);
 570 #endif
 571   IN_FLOAT (d = sqrt (d), "sqrt", arg);
 572   return make_float (d);
 573 }
 574 \f
 575 #if 0 /* Not clearly worth adding.  */
 576
 577 DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
 578        doc: /* Return the inverse hyperbolic cosine of ARG.  */)
 579   (register Lisp_Object arg)
 580 {
 581   double d = extract_float (arg);
 582 #ifdef FLOAT_CHECK_DOMAIN
 583   if (d < 1.0)
 584     domain_error ("acosh", arg);
 585 #endif
 586 #ifdef HAVE_INVERSE_HYPERBOLIC
 587   IN_FLOAT (d = acosh (d), "acosh", arg);
 588 #else
 589   IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
 590 #endif
 591   return make_float (d);
 592 }
 593
 594 DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
 595        doc: /* Return the inverse hyperbolic sine of ARG.  */)
 596   (register Lisp_Object arg)
 597 {
 598   double d = extract_float (arg);
 599 #ifdef HAVE_INVERSE_HYPERBOLIC
 600   IN_FLOAT (d = asinh (d), "asinh", arg);
 601 #else
 602   IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
 603 #endif
 604   return make_float (d);
 605 }
 606
 607 DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
 608        doc: /* Return the inverse hyperbolic tangent of ARG.  */)
 609   (register Lisp_Object arg)
 610 {
 611   double d = extract_float (arg);
 612 #ifdef FLOAT_CHECK_DOMAIN
 613   if (d >= 1.0 || d <= -1.0)
 614     domain_error ("atanh", arg);
 615 #endif
 616 #ifdef HAVE_INVERSE_HYPERBOLIC
 617   IN_FLOAT (d = atanh (d), "atanh", arg);
 618 #else
 619   IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
 620 #endif
 621   return make_float (d);
 622 }
 623
 624 DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
 625        doc: /* Return the hyperbolic cosine of ARG.  */)
 626   (register Lisp_Object arg)
 627 {
 628   double d = extract_float (arg);
 629 #ifdef FLOAT_CHECK_DOMAIN
 630   if (d > 710.0 || d < -710.0)
 631     range_error ("cosh", arg);
 632 #endif
 633   IN_FLOAT (d = cosh (d), "cosh", arg);
 634   return make_float (d);
 635 }
 636
 637 DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
 638        doc: /* Return the hyperbolic sine of ARG.  */)
 639   (register Lisp_Object arg)
 640 {
 641   double d = extract_float (arg);
 642 #ifdef FLOAT_CHECK_DOMAIN
 643   if (d > 710.0 || d < -710.0)
 644     range_error ("sinh", arg);
 645 #endif
 646   IN_FLOAT (d = sinh (d), "sinh", arg);
 647   return make_float (d);
 648 }
 649
 650 DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
 651        doc: /* Return the hyperbolic tangent of ARG.  */)
 652   (register Lisp_Object arg)
 653 {
 654   double d = extract_float (arg);
 655   IN_FLOAT (d = tanh (d), "tanh", arg);
 656   return make_float (d);
 657 }
 658 #endif
 659 \f
 660 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 661        doc: /* Return the absolute value of ARG.  */)
 662   (register Lisp_Object arg)
 663 {
 664   CHECK_NUMBER_OR_FLOAT (arg);
 665
 666   if (FLOATP (arg))
 667     arg = make_float (fabs (XFLOAT_DATA (arg)));
 668   else if (XINT (arg) < 0)
 669     XSETINT (arg, - XINT (arg));
 670
 671   return arg;
 672 }
 673
 674 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 675        doc: /* Return the floating point number equal to ARG.  */)
 676   (register Lisp_Object arg)
 677 {
 678   CHECK_NUMBER_OR_FLOAT (arg);
 679
 680   if (INTEGERP (arg))
 681     return make_float ((double) XINT (arg));
 682   else                          /* give 'em the same float back */
 683     return arg;
 684 }
 685
 686 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 687        doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
 688 This is the same as the exponent of a float.  */)
 689   (Lisp_Object arg)
 690 {
 691   Lisp_Object val;
 692   EMACS_INT value;
 693   double f = extract_float (arg);
 694
 695   if (f == 0.0)
 696     value = MOST_NEGATIVE_FIXNUM;
 697   else
 698     {
 699 #ifdef HAVE_LOGB
 700       IN_FLOAT (value = logb (f), "logb", arg);
 701 #else
 702 #ifdef HAVE_FREXP
 703       int ivalue;
 704       IN_FLOAT (frexp (f, &ivalue), "logb", arg);
 705       value = ivalue - 1;
 706 #else
 707       int i;
 708       double d;
 709       if (f < 0.0)
 710         f = -f;
 711       value = -1;
 712       while (f < 0.5)
 713         {
 714           for (i = 1, d = 0.5; d * d >= f; i += i)
 715             d *= d;
 716           f /= d;
 717           value -= i;
 718         }
 719       while (f >= 1.0)
 720         {
 721           for (i = 1, d = 2.0; d * d <= f; i += i)
 722             d *= d;
 723           f /= d;
 724           value += i;
 725         }
 726 #endif
 727 #endif
 728     }
 729   XSETINT (val, value);
 730   return val;
 731 }
 732
 733
 734 /* the rounding functions  */
 735
 736 static Lisp_Object
 737 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
 738                  double (*double_round) (double),
 739                  EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
 740                  const char *name)
 741 {
 742   CHECK_NUMBER_OR_FLOAT (arg);
 743
 744   if (! NILP (divisor))
 745     {
 746       EMACS_INT i1, i2;
 747
 748       CHECK_NUMBER_OR_FLOAT (divisor);
 749
 750       if (FLOATP (arg) || FLOATP (divisor))
 751         {
 752           double f1, f2;
 753
 754           f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 755           f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
 756           if (! IEEE_FLOATING_POINT && f2 == 0)
 757             xsignal0 (Qarith_error);
 758
 759           IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
 760           FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
 761           return arg;
 762         }
 763
 764       i1 = XINT (arg);
 765       i2 = XINT (divisor);
 766
 767       if (i2 == 0)
 768         xsignal0 (Qarith_error);
 769
 770       XSETINT (arg, (*int_round2) (i1, i2));
 771       return arg;
 772     }
 773
 774   if (FLOATP (arg))
 775     {
 776       double d;
 777
 778       IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
 779       FLOAT_TO_INT (d, arg, name, arg);
 780     }
 781
 782   return arg;
 783 }
 784
 785 /* With C's /, the result is implementation-defined if either operand
 786    is negative, so take care with negative operands in the following
 787    integer functions.  */
 788
 789 static EMACS_INT
 790 ceiling2 (EMACS_INT i1, EMACS_INT i2)
 791 {
 792   return (i2 < 0
 793           ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
 794           : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
 795 }
 796
 797 static EMACS_INT
 798 floor2 (EMACS_INT i1, EMACS_INT i2)
 799 {
 800   return (i2 < 0
 801           ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 802           : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 803 }
 804
 805 static EMACS_INT
 806 truncate2 (EMACS_INT i1, EMACS_INT i2)
 807 {
 808   return (i2 < 0
 809           ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
 810           : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
 811 }
 812
 813 static EMACS_INT
 814 round2 (EMACS_INT i1, EMACS_INT i2)
 815 {
 816   /* The C language's division operator gives us one remainder R, but
 817      we want the remainder R1 on the other side of 0 if R1 is closer
 818      to 0 than R is; because we want to round to even, we also want R1
 819      if R and R1 are the same distance from 0 and if C's quotient is
 820      odd.  */
 821   EMACS_INT q = i1 / i2;
 822   EMACS_INT r = i1 % i2;
 823   EMACS_INT abs_r = r < 0 ? -r : r;
 824   EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
 825   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 826 }
 827
 828 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 829    if `rint' exists but does not work right.  */
 830 #ifdef HAVE_RINT
 831 #define emacs_rint rint
 832 #else
 833 static double
 834 emacs_rint (double d)
 835 {
 836   return floor (d + 0.5);
 837 }
 838 #endif
 839
 840 static double
 841 double_identity (double d)
 842 {
 843   return d;
 844 }
 845
 846 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 847        doc: /* Return the smallest integer no less than ARG.
 848 This rounds the value towards +inf.
 849 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
 850   (Lisp_Object arg, Lisp_Object divisor)
 851 {
 852   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 853 }
 854
 855 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 856        doc: /* Return the largest integer no greater than ARG.
 857 This rounds the value towards -inf.
 858 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
 859   (Lisp_Object arg, Lisp_Object divisor)
 860 {
 861   return rounding_driver (arg, divisor, floor, floor2, "floor");
 862 }
 863
 864 DEFUN ("round", Fround, Sround, 1, 2, 0,
 865        doc: /* Return the nearest integer to ARG.
 866 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
 867
 868 Rounding a value equidistant between two integers may choose the
 869 integer closer to zero, or it may prefer an even integer, depending on
 870 your machine.  For example, \(round 2.5\) can return 3 on some
 871 systems, but 2 on others.  */)
 872   (Lisp_Object arg, Lisp_Object divisor)
 873 {
 874   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
 875 }
 876
 877 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
 878        doc: /* Truncate a floating point number to an int.
 879 Rounds ARG toward zero.
 880 With optional DIVISOR, truncate ARG/DIVISOR.  */)
 881   (Lisp_Object arg, Lisp_Object divisor)
 882 {
 883   return rounding_driver (arg, divisor, double_identity, truncate2,
 884                           "truncate");
 885 }
 886
 887
 888 Lisp_Object
 889 fmod_float (Lisp_Object x, Lisp_Object y)
 890 {
 891   double f1, f2;
 892
 893   f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
 894   f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
 895
 896   if (! IEEE_FLOATING_POINT && f2 == 0)
 897     xsignal0 (Qarith_error);
 898
 899   /* If the "remainder" comes out with the wrong sign, fix it.  */
 900   IN_FLOAT2 ((f1 = fmod (f1, f2),
 901               f1 = (f2 < 0 ? f1 > 0 : f1 < 0) ? f1 + f2 : f1),
 902              "mod", x, y);
 903   return make_float (f1);
 904 }
 905 \f
 906 /* It's not clear these are worth adding.  */
 907
 908 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 909        doc: /* Return the smallest integer no less than ARG, as a float.
 910 \(Round toward +inf.\)  */)
 911   (register Lisp_Object arg)
 912 {
 913   double d = extract_float (arg);
 914   IN_FLOAT (d = ceil (d), "fceiling", arg);
 915   return make_float (d);
 916 }
 917
 918 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 919        doc: /* Return the largest integer no greater than ARG, as a float.
 920 \(Round towards -inf.\)  */)
 921   (register Lisp_Object arg)
 922 {
 923   double d = extract_float (arg);
 924   IN_FLOAT (d = floor (d), "ffloor", arg);
 925   return make_float (d);
 926 }
 927
 928 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 929        doc: /* Return the nearest integer to ARG, as a float.  */)
 930   (register Lisp_Object arg)
 931 {
 932   double d = extract_float (arg);
 933   IN_FLOAT (d = emacs_rint (d), "fround", arg);
 934   return make_float (d);
 935 }
 936
 937 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 938        doc: /* Truncate a floating point number to an integral float value.
 939 Rounds the value toward zero.  */)
 940   (register Lisp_Object arg)
 941 {
 942   double d = extract_float (arg);
 943   if (d >= 0.0)
 944     IN_FLOAT (d = floor (d), "ftruncate", arg);
 945   else
 946     IN_FLOAT (d = ceil (d), "ftruncate", arg);
 947   return make_float (d);
 948 }
 949 \f
 950 #ifdef FLOAT_CATCH_SIGILL
 951 static void
 952 float_error (int signo)
 953 {
 954   if (! in_float)
 955     fatal_error_signal (signo);
 956
 957 #ifdef BSD_SYSTEM
 958   sigsetmask (SIGEMPTYMASK);
 959 #else
 960   /* Must reestablish handler each time it is called.  */
 961   signal (SIGILL, float_error);
 962 #endif /* BSD_SYSTEM */
 963
 964   SIGNAL_THREAD_CHECK (signo);
 965   in_float = 0;
 966
 967   xsignal1 (Qarith_error, float_error_arg);
 968 }
 969
 970 /* Another idea was to replace the library function `infnan'
 971    where SIGILL is signaled.  */
 972
 973 #endif /* FLOAT_CATCH_SIGILL */
 974
 975 #ifdef HAVE_MATHERR
 976 int
 977 matherr (struct exception *x)
 978 {
 979   Lisp_Object args;
 980   const char *name = x->name;
 981
 982   if (! in_float)
 983     /* Not called from emacs-lisp float routines; do the default thing. */
 984     return 0;
 985   if (!strcmp (x->name, "pow"))
 986     name = "expt";
 987
 988   args
 989     = Fcons (build_string (name),
 990              Fcons (make_float (x->arg1),
 991                     ((!strcmp (name, "log") || !strcmp (name, "pow"))
 992                      ? Fcons (make_float (x->arg2), Qnil)
 993                      : Qnil)));
 994   switch (x->type)
 995     {
 996     case DOMAIN:        xsignal (Qdomain_error, args);          break;
 997     case SING:          xsignal (Qsingularity_error, args);     break;
 998     case OVERFLOW:      xsignal (Qoverflow_error, args);        break;
 999     case UNDERFLOW:     xsignal (Qunderflow_error, args);       break;
1000     default:            xsignal (Qarith_error, args);           break;
1001     }
1002   return (1);   /* don't set errno or print a message */
1003 }
1004 #endif /* HAVE_MATHERR */
1005
1006 void
1007 init_floatfns (void)
1008 {
1009 #ifdef FLOAT_CATCH_SIGILL
1010   signal (SIGILL, float_error);
1011 #endif
1012   in_float = 0;
1013 }
1014
1015 void
1016 syms_of_floatfns (void)
1017 {
1018   defsubr (&Sacos);
1019   defsubr (&Sasin);
1020   defsubr (&Satan);
1021   defsubr (&Scos);
1022   defsubr (&Ssin);
1023   defsubr (&Stan);
1024   defsubr (&Sisnan);
1025 #ifdef HAVE_COPYSIGN
1026   defsubr (&Scopysign);
1027   defsubr (&Sfrexp);
1028   defsubr (&Sldexp);
1029 #endif
1030 #if 0
1031   defsubr (&Sacosh);
1032   defsubr (&Sasinh);
1033   defsubr (&Satanh);
1034   defsubr (&Scosh);
1035   defsubr (&Ssinh);
1036   defsubr (&Stanh);
1037   defsubr (&Sbessel_y0);
1038   defsubr (&Sbessel_y1);
1039   defsubr (&Sbessel_yn);
1040   defsubr (&Sbessel_j0);
1041   defsubr (&Sbessel_j1);
1042   defsubr (&Sbessel_jn);
1043   defsubr (&Serf);
1044   defsubr (&Serfc);
1045   defsubr (&Slog_gamma);
1046   defsubr (&Scube_root);
1047 #endif
1048   defsubr (&Sfceiling);
1049   defsubr (&Sffloor);
1050   defsubr (&Sfround);
1051   defsubr (&Sftruncate);
1052   defsubr (&Sexp);
1053   defsubr (&Sexpt);
1054   defsubr (&Slog);
1055   defsubr (&Slog10);
1056   defsubr (&Ssqrt);
1057
1058   defsubr (&Sabs);
1059   defsubr (&Sfloat);
1060   defsubr (&Slogb);
1061   defsubr (&Sceiling);
1062   defsubr (&Sfloor);
1063   defsubr (&Sround);
1064   defsubr (&Struncate);
1065 }