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