src/floatfns.c

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