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