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