src/floatfns.c

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