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