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