src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2    Copyright (C) 1988, 1993 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, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  19
  20
  21 /* ANSI C requires only these float functions:
  22    acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
  23    frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
  24
  25    Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
  26    Define HAVE_CBRT if you have cbrt.
  27    Define HAVE_RINT if you have rint.
  28    If you don't define these, then the appropriate routines will be simulated.
  29
  30    Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
  31    (This should happen automatically.)
  32
  33    Define FLOAT_CHECK_ERRNO if the float library routines set errno.
  34    This has no effect if HAVE_MATHERR is defined.
  35
  36    Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
  37    (What systems actually do this?  Please let us know.)
  38
  39    Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
  40    either setting errno, or signalling SIGFPE/SIGILL.  Otherwise, domain and
  41    range checking will happen before calling the float routines.  This has
  42    no effect if HAVE_MATHERR is defined (since matherr will be called when
  43    a domain error occurs.)
  44  */
  45
  46 #include <signal.h>
  47
  48 #include <config.h>
  49 #include "lisp.h"
  50 #include "syssignal.h"
  51
  52 Lisp_Object Qarith_error;
  53
  54 #ifdef LISP_FLOAT_TYPE
  55
  56 #include <math.h>
  57
  58 #ifndef hpux
  59 /* These declarations are omitted on some systems, like Ultrix.  */
  60 extern double logb ();
  61 #endif
  62
  63 #if defined(DOMAIN) && defined(SING) && defined(OVERFLOW)
  64     /* If those are defined, then this is probably a `matherr' machine. */
  65 # ifndef HAVE_MATHERR
  66 #  define HAVE_MATHERR
  67 # endif
  68 #endif
  69
  70 #ifdef NO_MATHERR
  71 #undef HAVE_MATHERR
  72 #endif
  73
  74 #ifdef HAVE_MATHERR
  75 # ifdef FLOAT_CHECK_ERRNO
  76 #  undef FLOAT_CHECK_ERRNO
  77 # endif
  78 # ifdef FLOAT_CHECK_DOMAIN
  79 #  undef FLOAT_CHECK_DOMAIN
  80 # endif
  81 #endif
  82
  83 #ifndef NO_FLOAT_CHECK_ERRNO
  84 #define FLOAT_CHECK_ERRNO
  85 #endif
  86
  87 #ifdef FLOAT_CHECK_ERRNO
  88 # include <errno.h>
  89
  90 extern int errno;
  91 #endif
  92
  93 /* Avoid traps on VMS from sinh and cosh.
  94    All the other functions set errno instead.  */
  95
  96 #ifdef VMS
  97 #undef cosh
  98 #undef sinh
  99 #define cosh(x) ((exp(x)+exp(-x))*0.5)
 100 #define sinh(x) ((exp(x)-exp(-x))*0.5)
 101 #endif /* VMS */
 102
 103 #ifndef HAVE_RINT
 104 #define rint(x) (floor((x)+0.5))
 105 #endif
 106
 107 static SIGTYPE float_error ();
 108
 109 /* Nonzero while executing in floating point.
 110    This tells float_error what to do.  */
 111
 112 static int in_float;
 113
 114 /* If an argument is out of range for a mathematical function,
 115    here is the actual argument value to use in the error message.  */
 116
 117 static Lisp_Object float_error_arg, float_error_arg2;
 118
 119 static char *float_error_fn_name;
 120
 121 /* Evaluate the floating point expression D, recording NUM
 122    as the original argument for error messages.
 123    D is normally an assignment expression.
 124    Handle errors which may result in signals or may set errno.
 125
 126    Note that float_error may be declared to return void, so you can't
 127    just cast the zero after the colon to (SIGTYPE) to make the types
 128    check properly.  */
 129
 130 #ifdef FLOAT_CHECK_ERRNO
 131 #define IN_FLOAT(d, name, num)                          \
 132   do {                                                  \
 133     float_error_arg = num;                              \
 134     float_error_fn_name = name;                         \
 135     in_float = 1; errno = 0; (d); in_float = 0;         \
 136     switch (errno) {                                    \
 137     case 0: break;                                      \
 138     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 139     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 140     default:     arith_error (float_error_fn_name, float_error_arg);    \
 141     }                                                   \
 142   } while (0)
 143 #define IN_FLOAT2(d, name, num, num2)                   \
 144   do {                                                  \
 145     float_error_arg = num;                              \
 146     float_error_arg2 = num2;                            \
 147     float_error_fn_name = name;                         \
 148     in_float = 1; errno = 0; (d); in_float = 0;         \
 149     switch (errno) {                                    \
 150     case 0: break;                                      \
 151     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 152     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 153     default:     arith_error (float_error_fn_name, float_error_arg);    \
 154     }                                                   \
 155   } while (0)
 156 #else
 157 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
 158 #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
 159 #endif
 160
 161 #define arith_error(op,arg) \
 162   Fsignal (Qarith_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
 163 #define range_error(op,arg) \
 164   Fsignal (Qrange_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
 165 #define domain_error(op,arg) \
 166   Fsignal (Qdomain_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
 167 #define domain_error2(op,a1,a2) \
 168   Fsignal (Qdomain_error, Fcons (build_string ((op)), Fcons ((a1), Fcons ((a2), Qnil))))
 169
 170 /* Extract a Lisp number as a `double', or signal an error.  */
 171
 172 double
 173 extract_float (num)
 174      Lisp_Object num;
 175 {
 176   CHECK_NUMBER_OR_FLOAT (num, 0);
 177
 178   if (XTYPE (num) == Lisp_Float)
 179     return XFLOAT (num)->data;
 180   return (double) XINT (num);
 181 }
 182 \f
 183 /* Trig functions.  */
 184
 185 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
 186   "Return the inverse cosine of ARG.")
 187   (arg)
 188      register Lisp_Object arg;
 189 {
 190   double d = extract_float (arg);
 191 #ifdef FLOAT_CHECK_DOMAIN
 192   if (d > 1.0 || d < -1.0)
 193     domain_error ("acos", arg);
 194 #endif
 195   IN_FLOAT (d = acos (d), "acos", arg);
 196   return make_float (d);
 197 }
 198
 199 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
 200   "Return the inverse sine of ARG.")
 201   (arg)
 202      register Lisp_Object arg;
 203 {
 204   double d = extract_float (arg);
 205 #ifdef FLOAT_CHECK_DOMAIN
 206   if (d > 1.0 || d < -1.0)
 207     domain_error ("asin", arg);
 208 #endif
 209   IN_FLOAT (d = asin (d), "asin", arg);
 210   return make_float (d);
 211 }
 212
 213 DEFUN ("atan", Fatan, Satan, 1, 1, 0,
 214   "Return the inverse tangent of ARG.")
 215   (arg)
 216      register Lisp_Object arg;
 217 {
 218   double d = extract_float (arg);
 219   IN_FLOAT (d = atan (d), "atan", arg);
 220   return make_float (d);
 221 }
 222
 223 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 224   "Return the cosine of ARG.")
 225   (arg)
 226      register Lisp_Object arg;
 227 {
 228   double d = extract_float (arg);
 229   IN_FLOAT (d = cos (d), "cos", arg);
 230   return make_float (d);
 231 }
 232
 233 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 234   "Return the sine of ARG.")
 235   (arg)
 236      register Lisp_Object arg;
 237 {
 238   double d = extract_float (arg);
 239   IN_FLOAT (d = sin (d), "sin", arg);
 240   return make_float (d);
 241 }
 242
 243 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 244   "Return the tangent of ARG.")
 245   (arg)
 246      register Lisp_Object arg;
 247 {
 248   double d = extract_float (arg);
 249   double c = cos (d);
 250 #ifdef FLOAT_CHECK_DOMAIN
 251   if (c == 0.0)
 252     domain_error ("tan", arg);
 253 #endif
 254   IN_FLOAT (d = sin (d) / c, "tan", arg);
 255   return make_float (d);
 256 }
 257 \f
 258 #if 0 /* Leave these out unless we find there's a reason for them.  */
 259
 260 DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
 261   "Return the bessel function j0 of ARG.")
 262   (arg)
 263      register Lisp_Object arg;
 264 {
 265   double d = extract_float (arg);
 266   IN_FLOAT (d = j0 (d), "bessel-j0", arg);
 267   return make_float (d);
 268 }
 269
 270 DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
 271   "Return the bessel function j1 of ARG.")
 272   (arg)
 273      register Lisp_Object arg;
 274 {
 275   double d = extract_float (arg);
 276   IN_FLOAT (d = j1 (d), "bessel-j1", arg);
 277   return make_float (d);
 278 }
 279
 280 DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
 281   "Return the order N bessel function output jn of ARG.\n\
 282 The first arg (the order) is truncated to an integer.")
 283   (arg1, arg2)
 284      register Lisp_Object arg1, arg2;
 285 {
 286   int i1 = extract_float (arg1);
 287   double f2 = extract_float (arg2);
 288
 289   IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", arg1);
 290   return make_float (f2);
 291 }
 292
 293 DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
 294   "Return the bessel function y0 of ARG.")
 295   (arg)
 296      register Lisp_Object arg;
 297 {
 298   double d = extract_float (arg);
 299   IN_FLOAT (d = y0 (d), "bessel-y0", arg);
 300   return make_float (d);
 301 }
 302
 303 DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
 304   "Return the bessel function y1 of ARG.")
 305   (arg)
 306      register Lisp_Object arg;
 307 {
 308   double d = extract_float (arg);
 309   IN_FLOAT (d = y1 (d), "bessel-y0", arg);
 310   return make_float (d);
 311 }
 312
 313 DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
 314   "Return the order N bessel function output yn of ARG.\n\
 315 The first arg (the order) is truncated to an integer.")
 316   (arg1, arg2)
 317      register Lisp_Object arg1, arg2;
 318 {
 319   int i1 = extract_float (arg1);
 320   double f2 = extract_float (arg2);
 321
 322   IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", arg1);
 323   return make_float (f2);
 324 }
 325
 326 #endif
 327 \f
 328 #if 0 /* Leave these out unless we see they are worth having.  */
 329
 330 DEFUN ("erf", Ferf, Serf, 1, 1, 0,
 331   "Return the mathematical error function of ARG.")
 332   (arg)
 333      register Lisp_Object arg;
 334 {
 335   double d = extract_float (arg);
 336   IN_FLOAT (d = erf (d), "erf", arg);
 337   return make_float (d);
 338 }
 339
 340 DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
 341   "Return the complementary error function of ARG.")
 342   (arg)
 343      register Lisp_Object arg;
 344 {
 345   double d = extract_float (arg);
 346   IN_FLOAT (d = erfc (d), "erfc", arg);
 347   return make_float (d);
 348 }
 349
 350 DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
 351   "Return the log gamma of ARG.")
 352   (arg)
 353      register Lisp_Object arg;
 354 {
 355   double d = extract_float (arg);
 356   IN_FLOAT (d = lgamma (d), "log-gamma", arg);
 357   return make_float (d);
 358 }
 359
 360 DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
 361   "Return the cube root of ARG.")
 362   (arg)
 363      register Lisp_Object arg;
 364 {
 365   double d = extract_float (arg);
 366 #ifdef HAVE_CBRT
 367   IN_FLOAT (d = cbrt (d), "cube-root", arg);
 368 #else
 369   if (d >= 0.0)
 370     IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
 371   else
 372     IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
 373 #endif
 374   return make_float (d);
 375 }
 376
 377 #endif
 378 \f
 379 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 380   "Return the exponential base e of ARG.")
 381   (arg)
 382      register Lisp_Object arg;
 383 {
 384   double d = extract_float (arg);
 385 #ifdef FLOAT_CHECK_DOMAIN
 386   if (d > 709.7827)   /* Assume IEEE doubles here */
 387     range_error ("exp", arg);
 388   else if (d < -709.0)
 389     return make_float (0.0);
 390   else
 391 #endif
 392     IN_FLOAT (d = exp (d), "exp", arg);
 393   return make_float (d);
 394 }
 395
 396 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 397   "Return the exponential X ** Y.")
 398   (arg1, arg2)
 399      register Lisp_Object arg1, arg2;
 400 {
 401   double f1, f2;
 402
 403   CHECK_NUMBER_OR_FLOAT (arg1, 0);
 404   CHECK_NUMBER_OR_FLOAT (arg2, 0);
 405   if (XTYPE (arg1) == Lisp_Int     /* common lisp spec */
 406       && XTYPE (arg2) == Lisp_Int) /* don't promote, if both are ints */
 407     {                           /* this can be improved by pre-calculating */
 408       int acc, x, y;            /* some binary powers of x then accumulating */
 409       Lisp_Object val;
 410
 411       x = XINT (arg1);
 412       y = XINT (arg2);
 413       acc = 1;
 414
 415       if (y < 0)
 416         {
 417           if (x == 1)
 418             acc = 1;
 419           else if (x == -1)
 420             acc = (y & 1) ? -1 : 1;
 421           else
 422             acc = 0;
 423         }
 424       else
 425         {
 426           for (; y > 0; y--)
 427           while (y > 0)
 428             {
 429               if (y & 1)
 430                 acc *= x;
 431               x *= x;
 432               y = (unsigned)y >> 1;
 433             }
 434         }
 435       XSET (val, Lisp_Int, acc);
 436       return val;
 437     }
 438   f1 = (XTYPE (arg1) == Lisp_Float) ? XFLOAT (arg1)->data : XINT (arg1);
 439   f2 = (XTYPE (arg2) == Lisp_Float) ? XFLOAT (arg2)->data : XINT (arg2);
 440   /* Really should check for overflow, too */
 441   if (f1 == 0.0 && f2 == 0.0)
 442     f1 = 1.0;
 443 #ifdef FLOAT_CHECK_DOMAIN
 444   else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor(f2)))
 445     domain_error2 ("expt", arg1, arg2);
 446 #endif
 447   IN_FLOAT2 (f1 = pow (f1, f2), "expt", arg1, arg2);
 448   return make_float (f1);
 449 }
 450
 451 DEFUN ("log", Flog, Slog, 1, 2, 0,
 452   "Return the natural logarithm of ARG.\n\
 453 If second optional argument BASE is given, return log ARG using that base.")
 454   (arg, base)
 455      register Lisp_Object arg, base;
 456 {
 457   double d = extract_float (arg);
 458
 459 #ifdef FLOAT_CHECK_DOMAIN
 460   if (d <= 0.0)
 461     domain_error2 ("log", arg, base);
 462 #endif
 463   if (NILP (base))
 464     IN_FLOAT (d = log (d), "log", arg);
 465   else
 466     {
 467       double b = extract_float (base);
 468
 469 #ifdef FLOAT_CHECK_DOMAIN
 470       if (b <= 0.0 || b == 1.0)
 471         domain_error2 ("log", arg, base);
 472 #endif
 473       if (b == 10.0)
 474         IN_FLOAT2 (d = log10 (d), "log", arg, base);
 475       else
 476         IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
 477     }
 478   return make_float (d);
 479 }
 480
 481 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
 482   "Return the logarithm base 10 of ARG.")
 483   (arg)
 484      register Lisp_Object arg;
 485 {
 486   double d = extract_float (arg);
 487 #ifdef FLOAT_CHECK_DOMAIN
 488   if (d <= 0.0)
 489     domain_error ("log10", arg);
 490 #endif
 491   IN_FLOAT (d = log10 (d), "log10", arg);
 492   return make_float (d);
 493 }
 494
 495 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 496   "Return the square root of ARG.")
 497   (arg)
 498      register Lisp_Object arg;
 499 {
 500   double d = extract_float (arg);
 501 #ifdef FLOAT_CHECK_DOMAIN
 502   if (d < 0.0)
 503     domain_error ("sqrt", arg);
 504 #endif
 505   IN_FLOAT (d = sqrt (d), "sqrt", arg);
 506   return make_float (d);
 507 }
 508 \f
 509 #if 0 /* Not clearly worth adding.  */
 510
 511 DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
 512   "Return the inverse hyperbolic cosine of ARG.")
 513   (arg)
 514      register Lisp_Object arg;
 515 {
 516   double d = extract_float (arg);
 517 #ifdef FLOAT_CHECK_DOMAIN
 518   if (d < 1.0)
 519     domain_error ("acosh", arg);
 520 #endif
 521 #ifdef HAVE_INVERSE_HYPERBOLIC
 522   IN_FLOAT (d = acosh (d), "acosh", arg);
 523 #else
 524   IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
 525 #endif
 526   return make_float (d);
 527 }
 528
 529 DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
 530   "Return the inverse hyperbolic sine of ARG.")
 531   (arg)
 532      register Lisp_Object arg;
 533 {
 534   double d = extract_float (arg);
 535 #ifdef HAVE_INVERSE_HYPERBOLIC
 536   IN_FLOAT (d = asinh (d), "asinh", arg);
 537 #else
 538   IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
 539 #endif
 540   return make_float (d);
 541 }
 542
 543 DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
 544   "Return the inverse hyperbolic tangent of ARG.")
 545   (arg)
 546      register Lisp_Object arg;
 547 {
 548   double d = extract_float (arg);
 549 #ifdef FLOAT_CHECK_DOMAIN
 550   if (d >= 1.0 || d <= -1.0)
 551     domain_error ("atanh", arg);
 552 #endif
 553 #ifdef HAVE_INVERSE_HYPERBOLIC
 554   IN_FLOAT (d = atanh (d), "atanh", arg);
 555 #else
 556   IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
 557 #endif
 558   return make_float (d);
 559 }
 560
 561 DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
 562   "Return the hyperbolic cosine of ARG.")
 563   (arg)
 564      register Lisp_Object arg;
 565 {
 566   double d = extract_float (arg);
 567 #ifdef FLOAT_CHECK_DOMAIN
 568   if (d > 710.0 || d < -710.0)
 569     range_error ("cosh", arg);
 570 #endif
 571   IN_FLOAT (d = cosh (d), "cosh", arg);
 572   return make_float (d);
 573 }
 574
 575 DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
 576   "Return the hyperbolic sine of ARG.")
 577   (arg)
 578      register Lisp_Object arg;
 579 {
 580   double d = extract_float (arg);
 581 #ifdef FLOAT_CHECK_DOMAIN
 582   if (d > 710.0 || d < -710.0)
 583     range_error ("sinh", arg);
 584 #endif
 585   IN_FLOAT (d = sinh (d), "sinh", arg);
 586   return make_float (d);
 587 }
 588
 589 DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
 590   "Return the hyperbolic tangent of ARG.")
 591   (arg)
 592      register Lisp_Object arg;
 593 {
 594   double d = extract_float (arg);
 595   IN_FLOAT (d = tanh (d), "tanh", arg);
 596   return make_float (d);
 597 }
 598 #endif
 599 \f
 600 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 601   "Return the absolute value of ARG.")
 602   (arg)
 603      register Lisp_Object arg;
 604 {
 605   CHECK_NUMBER_OR_FLOAT (arg, 0);
 606
 607   if (XTYPE (arg) == Lisp_Float)
 608     IN_FLOAT (arg = make_float (fabs (XFLOAT (arg)->data)), "abs", arg);
 609   else if (XINT (arg) < 0)
 610     XSETINT (arg, - XFASTINT (arg));
 611
 612   return arg;
 613 }
 614
 615 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 616   "Return the floating point number equal to ARG.")
 617   (arg)
 618      register Lisp_Object arg;
 619 {
 620   CHECK_NUMBER_OR_FLOAT (arg, 0);
 621
 622   if (XTYPE (arg) == Lisp_Int)
 623     return make_float ((double) XINT (arg));
 624   else                          /* give 'em the same float back */
 625     return arg;
 626 }
 627
 628 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 629   "Returns largest integer <= the base 2 log of the magnitude of ARG.\n\
 630 This is the same as the exponent of a float.")
 631      (arg)
 632      Lisp_Object arg;
 633 {
 634   Lisp_Object val;
 635   int value;
 636   double f = extract_float (arg);
 637
 638 #ifdef HAVE_LOGB
 639   IN_FLOAT (value = logb (f), "logb", arg);
 640   XSET (val, Lisp_Int, value);
 641 #else
 642 #ifdef HAVE_FREXP
 643   {
 644     int exp;
 645
 646     IN_FLOAT (frexp (f, &exp), "logb", arg);
 647     XSET (val, Lisp_Int, exp-1);
 648   }
 649 #else
 650   /* Would someone like to write code to emulate logb?  */
 651   error ("`logb' not implemented on this operating system");
 652 #endif
 653 #endif
 654
 655   return val;
 656 }
 657
 658 /* the rounding functions  */
 659
 660 DEFUN ("ceiling", Fceiling, Sceiling, 1, 1, 0,
 661   "Return the smallest integer no less than ARG.  (Round toward +inf.)")
 662   (arg)
 663      register Lisp_Object arg;
 664 {
 665   CHECK_NUMBER_OR_FLOAT (arg, 0);
 666
 667   if (XTYPE (arg) == Lisp_Float)
 668     IN_FLOAT (XSET (arg, Lisp_Int, ceil (XFLOAT (arg)->data)), "ceiling", arg);
 669
 670   return arg;
 671 }
 672
 673 #endif /* LISP_FLOAT_TYPE */
 674
 675
 676 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 677   "Return the largest integer no greater than ARG.  (Round towards -inf.)\n\
 678 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.")
 679   (arg, divisor)
 680      register Lisp_Object arg, divisor;
 681 {
 682   CHECK_NUMBER_OR_FLOAT (arg, 0);
 683
 684   if (! NILP (divisor))
 685     {
 686       int i1, i2;
 687
 688       CHECK_NUMBER_OR_FLOAT (divisor, 1);
 689
 690 #ifdef LISP_FLOAT_TYPE
 691       if (XTYPE (arg) == Lisp_Float || XTYPE (divisor) == Lisp_Float)
 692         {
 693           double f1, f2;
 694
 695           f1 = XTYPE (arg) == Lisp_Float ? XFLOAT (arg)->data : XINT (arg);
 696           f2 = (XTYPE (divisor) == Lisp_Float
 697                 ? XFLOAT (divisor)->data : XINT (divisor));
 698           if (f2 == 0)
 699             Fsignal (Qarith_error, Qnil);
 700
 701           IN_FLOAT2 (XSET (arg, Lisp_Int, floor (f1 / f2)),
 702                      "floor", arg, divisor);
 703           return arg;
 704         }
 705 #endif
 706
 707       i1 = XINT (arg);
 708       i2 = XINT (divisor);
 709
 710       if (i2 == 0)
 711         Fsignal (Qarith_error, Qnil);
 712
 713       /* With C's /, the result is implementation-defined if either operand
 714          is negative, so use only nonnegative operands.  */
 715       i1 = (i2 < 0
 716             ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 717             : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 718
 719       XSET (arg, Lisp_Int, i1);
 720       return arg;
 721     }
 722
 723 #ifdef LISP_FLOAT_TYPE
 724   if (XTYPE (arg) == Lisp_Float)
 725     IN_FLOAT (XSET (arg, Lisp_Int, floor (XFLOAT (arg)->data)), "floor", arg);
 726 #endif
 727
 728   return arg;
 729 }
 730
 731 #ifdef LISP_FLOAT_TYPE
 732
 733 DEFUN ("round", Fround, Sround, 1, 1, 0,
 734   "Return the nearest integer to ARG.")
 735   (arg)
 736      register Lisp_Object arg;
 737 {
 738   CHECK_NUMBER_OR_FLOAT (arg, 0);
 739
 740   if (XTYPE (arg) == Lisp_Float)
 741     /* Screw the prevailing rounding mode.  */
 742     IN_FLOAT (XSET (arg, Lisp_Int, rint (XFLOAT (arg)->data)), "round", arg);
 743
 744   return arg;
 745 }
 746
 747 DEFUN ("truncate", Ftruncate, Struncate, 1, 1, 0,
 748        "Truncate a floating point number to an int.\n\
 749 Rounds the value toward zero.")
 750   (arg)
 751      register Lisp_Object arg;
 752 {
 753   CHECK_NUMBER_OR_FLOAT (arg, 0);
 754
 755   if (XTYPE (arg) == Lisp_Float)
 756     XSET (arg, Lisp_Int, (int) XFLOAT (arg)->data);
 757
 758   return arg;
 759 }
 760 \f
 761 #if 0
 762 /* It's not clear these are worth adding.  */
 763
 764 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
 765   "Return the smallest integer no less than ARG, as a float.\n\
 766 \(Round toward +inf.\)")
 767   (arg)
 768      register Lisp_Object arg;
 769 {
 770   double d = extract_float (arg);
 771   IN_FLOAT (d = ceil (d), "fceiling", arg);
 772   return make_float (d);
 773 }
 774
 775 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
 776   "Return the largest integer no greater than ARG, as a float.\n\
 777 \(Round towards -inf.\)")
 778   (arg)
 779      register Lisp_Object arg;
 780 {
 781   double d = extract_float (arg);
 782   IN_FLOAT (d = floor (d), "ffloor", arg);
 783   return make_float (d);
 784 }
 785
 786 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
 787   "Return the nearest integer to ARG, as a float.")
 788   (arg)
 789      register Lisp_Object arg;
 790 {
 791   double d = extract_float (arg);
 792   IN_FLOAT (d = rint (XFLOAT (arg)->data), "fround", arg);
 793   return make_float (d);
 794 }
 795
 796 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
 797        "Truncate a floating point number to an integral float value.\n\
 798 Rounds the value toward zero.")
 799   (arg)
 800      register Lisp_Object arg;
 801 {
 802   double d = extract_float (arg);
 803   if (d >= 0.0)
 804     IN_FLOAT (d = floor (d), "ftruncate", arg);
 805   else
 806     IN_FLOAT (d = ceil (d), arg);
 807   return make_float (d);
 808 }
 809 #endif
 810 \f
 811 #ifdef FLOAT_CATCH_SIGILL
 812 static SIGTYPE
 813 float_error (signo)
 814      int signo;
 815 {
 816   if (! in_float)
 817     fatal_error_signal (signo);
 818
 819 #ifdef BSD
 820 #ifdef BSD4_1
 821   sigrelse (SIGILL);
 822 #else /* not BSD4_1 */
 823   sigsetmask (SIGEMPTYMASK);
 824 #endif /* not BSD4_1 */
 825 #else
 826   /* Must reestablish handler each time it is called.  */
 827   signal (SIGILL, float_error);
 828 #endif /* BSD */
 829
 830   in_float = 0;
 831
 832   Fsignal (Qarith_error, Fcons (float_error_arg, Qnil));
 833 }
 834
 835 /* Another idea was to replace the library function `infnan'
 836    where SIGILL is signaled.  */
 837
 838 #endif /* FLOAT_CATCH_SIGILL */
 839
 840 #ifdef HAVE_MATHERR
 841 int
 842 matherr (x)
 843      struct exception *x;
 844 {
 845   Lisp_Object args;
 846   if (! in_float)
 847     /* Not called from emacs-lisp float routines; do the default thing. */
 848     return 0;
 849   if (!strcmp (x->name, "pow"))
 850     x->name = "expt";
 851
 852   args
 853     = Fcons (build_string (x->name),
 854              Fcons (make_float (x->arg1),
 855                     ((!strcmp (x->name, "log") || !strcmp (x->name, "pow"))
 856                      ? Fcons (make_float (x->arg2), Qnil)
 857                      : Qnil)));
 858   switch (x->type)
 859     {
 860     case DOMAIN:        Fsignal (Qdomain_error, args);          break;
 861     case SING:          Fsignal (Qsingularity_error, args);     break;
 862     case OVERFLOW:      Fsignal (Qoverflow_error, args);        break;
 863     case UNDERFLOW:     Fsignal (Qunderflow_error, args);       break;
 864     default:            Fsignal (Qarith_error, args);           break;
 865     }
 866   return (1);   /* don't set errno or print a message */
 867 }
 868 #endif /* HAVE_MATHERR */
 869
 870 init_floatfns ()
 871 {
 872 #ifdef FLOAT_CATCH_SIGILL
 873   signal (SIGILL, float_error);
 874 #endif
 875   in_float = 0;
 876 }
 877
 878 #else /* not LISP_FLOAT_TYPE */
 879
 880 init_floatfns ()
 881 {}
 882
 883 #endif /* not LISP_FLOAT_TYPE */
 884
 885 syms_of_floatfns ()
 886 {
 887 #ifdef LISP_FLOAT_TYPE
 888   defsubr (&Sacos);
 889   defsubr (&Sasin);
 890   defsubr (&Satan);
 891   defsubr (&Scos);
 892   defsubr (&Ssin);
 893   defsubr (&Stan);
 894 #if 0
 895   defsubr (&Sacosh);
 896   defsubr (&Sasinh);
 897   defsubr (&Satanh);
 898   defsubr (&Scosh);
 899   defsubr (&Ssinh);
 900   defsubr (&Stanh);
 901   defsubr (&Sbessel_y0);
 902   defsubr (&Sbessel_y1);
 903   defsubr (&Sbessel_yn);
 904   defsubr (&Sbessel_j0);
 905   defsubr (&Sbessel_j1);
 906   defsubr (&Sbessel_jn);
 907   defsubr (&Serf);
 908   defsubr (&Serfc);
 909   defsubr (&Slog_gamma);
 910   defsubr (&Scube_root);
 911   defsubr (&Sfceiling);
 912   defsubr (&Sffloor);
 913   defsubr (&Sfround);
 914   defsubr (&Sftruncate);
 915 #endif
 916   defsubr (&Sexp);
 917   defsubr (&Sexpt);
 918   defsubr (&Slog);
 919   defsubr (&Slog10);
 920   defsubr (&Ssqrt);
 921
 922   defsubr (&Sabs);
 923   defsubr (&Sfloat);
 924   defsubr (&Slogb);
 925   defsubr (&Sceiling);
 926   defsubr (&Sround);
 927   defsubr (&Struncate);
 928 #endif /* LISP_FLOAT_TYPE */
 929   defsubr (&Sfloor);
 930 }