src/floatfns.c

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