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