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