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