src/floatfns.c

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