src/floatfns.c

   1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
   2    Copyright (C) 1988, 1993, 1994, 1999 Free Software Foundation, Inc.
   3
   4 This file is part of GNU Emacs.
   5
   6 GNU Emacs is free software; you can redistribute it and/or modify
   7 it under the terms of the GNU General Public License as published by
   8 the Free Software Foundation; either version 2, or (at your option)
   9 any later version.
  10
  11 GNU Emacs is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 GNU General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with GNU Emacs; see the file COPYING.  If not, write to
  18 the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19 Boston, MA 02111-1307, USA.  */
  20
  21
  22 /* ANSI C requires only these float functions:
  23    acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
  24    frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
  25
  26    Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
  27    Define HAVE_CBRT if you have cbrt.
  28    Define HAVE_RINT if you have a working rint.
  29    If you don't define these, then the appropriate routines will be simulated.
  30
  31    Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
  32    (This should happen automatically.)
  33
  34    Define FLOAT_CHECK_ERRNO if the float library routines set errno.
  35    This has no effect if HAVE_MATHERR is defined.
  36
  37    Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
  38    (What systems actually do this?  Please let us know.)
  39
  40    Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
  41    either setting errno, or signaling SIGFPE/SIGILL.  Otherwise, domain and
  42    range checking will happen before calling the float routines.  This has
  43    no effect if HAVE_MATHERR is defined (since matherr will be called when
  44    a domain error occurs.)
  45  */
  46
  47 #include <config.h>
  48 #include <signal.h>
  49 #include "lisp.h"
  50 #include "syssignal.h"
  51
  52 #if STDC_HEADERS
  53 #include <float.h>
  54 #endif
  55
  56 /* If IEEE_FLOATING_POINT isn't defined, default it from FLT_*. */
  57 #ifndef IEEE_FLOATING_POINT
  58 #if (FLT_RADIX == 2 && FLT_MANT_DIG == 24 \
  59      && FLT_MIN_EXP == -125 && FLT_MAX_EXP == 128)
  60 #define IEEE_FLOATING_POINT 1
  61 #else
  62 #define IEEE_FLOATING_POINT 0
  63 #endif
  64 #endif
  65
  66 /* Work around a problem that happens because math.h on hpux 7
  67    defines two static variables--which, in Emacs, are not really static,
  68    because `static' is defined as nothing.  The problem is that they are
  69    defined both here and in lread.c.
  70    These macros prevent the name conflict.  */
  71 #if defined (HPUX) && !defined (HPUX8)
  72 #define _MAXLDBL floatfns_maxldbl
  73 #define _NMAXLDBL floatfns_nmaxldbl
  74 #endif
  75
  76 #include <math.h>
  77
  78 /* This declaration is omitted on some systems, like Ultrix.  */
  79 #if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
  80 extern double logb ();
  81 #endif /* not HPUX and HAVE_LOGB and no logb macro */
  82
  83 #if defined(DOMAIN) && defined(SING) && defined(OVERFLOW)
  84     /* If those are defined, then this is probably a `matherr' machine. */
  85 # ifndef HAVE_MATHERR
  86 #  define HAVE_MATHERR
  87 # endif
  88 #endif
  89
  90 #ifdef NO_MATHERR
  91 #undef HAVE_MATHERR
  92 #endif
  93
  94 #ifdef HAVE_MATHERR
  95 # ifdef FLOAT_CHECK_ERRNO
  96 #  undef FLOAT_CHECK_ERRNO
  97 # endif
  98 # ifdef FLOAT_CHECK_DOMAIN
  99 #  undef FLOAT_CHECK_DOMAIN
 100 # endif
 101 #endif
 102
 103 #ifndef NO_FLOAT_CHECK_ERRNO
 104 #define FLOAT_CHECK_ERRNO
 105 #endif
 106
 107 #ifdef FLOAT_CHECK_ERRNO
 108 # include <errno.h>
 109
 110 #ifndef USE_CRT_DLL
 111 extern int errno;
 112 #endif
 113 #endif
 114
 115 /* Avoid traps on VMS from sinh and cosh.
 116    All the other functions set errno instead.  */
 117
 118 #ifdef VMS
 119 #undef cosh
 120 #undef sinh
 121 #define cosh(x) ((exp(x)+exp(-x))*0.5)
 122 #define sinh(x) ((exp(x)-exp(-x))*0.5)
 123 #endif /* VMS */
 124
 125 #ifdef FLOAT_CATCH_SIGILL
 126 static SIGTYPE float_error ();
 127 #endif
 128
 129 /* Nonzero while executing in floating point.
 130    This tells float_error what to do.  */
 131
 132 static int in_float;
 133
 134 /* If an argument is out of range for a mathematical function,
 135    here is the actual argument value to use in the error message.
 136    These variables are used only across the floating point library call
 137    so there is no need to staticpro them.  */
 138
 139 static Lisp_Object float_error_arg, float_error_arg2;
 140
 141 static char *float_error_fn_name;
 142
 143 /* Evaluate the floating point expression D, recording NUM
 144    as the original argument for error messages.
 145    D is normally an assignment expression.
 146    Handle errors which may result in signals or may set errno.
 147
 148    Note that float_error may be declared to return void, so you can't
 149    just cast the zero after the colon to (SIGTYPE) to make the types
 150    check properly.  */
 151
 152 #ifdef FLOAT_CHECK_ERRNO
 153 #define IN_FLOAT(d, name, num)                          \
 154   do {                                                  \
 155     float_error_arg = num;                              \
 156     float_error_fn_name = name;                         \
 157     in_float = 1; errno = 0; (d); in_float = 0;         \
 158     switch (errno) {                                    \
 159     case 0: break;                                      \
 160     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 161     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 162     default:     arith_error (float_error_fn_name, float_error_arg);    \
 163     }                                                   \
 164   } while (0)
 165 #define IN_FLOAT2(d, name, num, num2)                   \
 166   do {                                                  \
 167     float_error_arg = num;                              \
 168     float_error_arg2 = num2;                            \
 169     float_error_fn_name = name;                         \
 170     in_float = 1; errno = 0; (d); in_float = 0;         \
 171     switch (errno) {                                    \
 172     case 0: break;                                      \
 173     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 174     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 175     default:     arith_error (float_error_fn_name, float_error_arg);    \
 176     }                                                   \
 177   } while (0)
 178 #else
 179 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
 180 #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
 181 #endif
 182
 183 /* Convert float to Lisp_Int if it fits, else signal a range error
 184    using the given arguments.  */
 185 #define FLOAT_TO_INT(x, i, name, num)                                   \
 186   do                                                                    \
 187     {                                                                   \
 188       if (FIXNUM_OVERFLOW_P (x))                                        \
 189         range_error (name, num);                                        \
 190       XSETINT (i,  (EMACS_INT)(x));                                     \
 191     }                                                                   \
 192   while (0)
 193 #define FLOAT_TO_INT2(x, i, name, num1, num2)                           \
 194   do                                                                    \
 195     {                                                                   \
 196       if (FIXNUM_OVERFLOW_P (x))                                        \
 197         range_error2 (name, num1, num2);                                \
 198       XSETINT (i,  (EMACS_INT)(x));                                     \
 199     }                                                                   \
 200   while (0)
 201
 202 #define arith_error(op,arg) \
 203   Fsignal (Qarith_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
 204 #define range_error(op,arg) \
 205   Fsignal (Qrange_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
 206 #define range_error2(op,a1,a2) \
 207   Fsignal (Qrange_error, Fcons (build_string ((op)), \
 208                                 Fcons ((a1), Fcons ((a2), Qnil))))
 209 #define domain_error(op,arg) \
 210   Fsignal (Qdomain_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
 211 #define domain_error2(op,a1,a2) \
 212   Fsignal (Qdomain_error, Fcons (build_string ((op)), \
 213                                  Fcons ((a1), Fcons ((a2), Qnil))))
 214
 215 /* Extract a Lisp number as a `double', or signal an error.  */
 216
 217 double
 218 extract_float (num)
 219      Lisp_Object num;
 220 {
 221   CHECK_NUMBER_OR_FLOAT (num);
 222
 223   if (FLOATP (num))
 224     return XFLOAT_DATA (num);
 225   return (double) XINT (num);
 226 }
 227 \f
 228 /* Trig functions.  */
 229
 230 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
 231        doc: /* Return the inverse cosine of ARG.  */)
 232      (arg)
 233      register Lisp_Object arg;
 234 {
 235   double d = extract_float (arg);
 236 #ifdef FLOAT_CHECK_DOMAIN
 237   if (d > 1.0 || d < -1.0)
 238     domain_error ("acos", arg);
 239 #endif
 240   IN_FLOAT (d = acos (d), "acos", arg);
 241   return make_float (d);
 242 }
 243
 244 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
 245        doc: /* Return the inverse sine of ARG.  */)
 246      (arg)
 247      register Lisp_Object arg;
 248 {
 249   double d = extract_float (arg);
 250 #ifdef FLOAT_CHECK_DOMAIN
 251   if (d > 1.0 || d < -1.0)
 252     domain_error ("asin", arg);
 253 #endif
 254   IN_FLOAT (d = asin (d), "asin", arg);
 255   return make_float (d);
 256 }
 257
 258 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
 259        doc: /* Return the inverse tangent of the arguments.
 260 If only one argument Y is given, return the inverse tangent of Y.
 261 If two arguments Y and X are given, return the inverse tangent of Y
 262 divided by X, i.e. the angle in radians between the vector (X, Y)
 263 and the x-axis.  */)
 264      (y, x)
 265      register Lisp_Object y, x;
 266 {
 267   double d = extract_float (y);
 268
 269   if (NILP (x))
 270     IN_FLOAT (d = atan (d), "atan", y);
 271   else
 272     {
 273       double d2 = extract_float (x);
 274
 275       IN_FLOAT2 (d = atan2 (d, d2), "atan", y, x);
 276     }
 277   return make_float (d);
 278 }
 279
 280 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 281        doc: /* Return the cosine of ARG.  */)
 282      (arg)
 283      register Lisp_Object arg;
 284 {
 285   double d = extract_float (arg);
 286   IN_FLOAT (d = cos (d), "cos", arg);
 287   return make_float (d);
 288 }
 289
 290 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 291        doc: /* Return the sine of ARG.  */)
 292      (arg)
 293      register Lisp_Object arg;
 294 {
 295   double d = extract_float (arg);
 296   IN_FLOAT (d = sin (d), "sin", arg);
 297   return make_float (d);
 298 }
 299
 300 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 301        doc: /* Return the tangent of ARG.  */)
 302      (arg)
 303      register Lisp_Object arg;
 304 {
 305   double d = extract_float (arg);
 306   double c = cos (d);
 307 #ifdef FLOAT_CHECK_DOMAIN
 308   if (c == 0.0)
 309     domain_error ("tan", arg);
 310 #endif
 311   IN_FLOAT (d = sin (d) / c, "tan", arg);
 312   return make_float (d);
 313 }
 314 \f
 315 #if 0 /* Leave these out unless we find there's a reason for them.  */
 316
 317 DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
 318        doc: /* Return the bessel function j0 of ARG.  */)
 319      (arg)
 320      register Lisp_Object arg;
 321 {
 322   double d = extract_float (arg);
 323   IN_FLOAT (d = j0 (d), "bessel-j0", arg);
 324   return make_float (d);
 325 }
 326
 327 DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
 328        doc: /* Return the bessel function j1 of ARG.  */)
 329      (arg)
 330      register Lisp_Object arg;
 331 {
 332   double d = extract_float (arg);
 333   IN_FLOAT (d = j1 (d), "bessel-j1", arg);
 334   return make_float (d);
 335 }
 336
 337 DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
 338        doc: /* Return the order N bessel function output jn of ARG.
 339 The first arg (the order) is truncated to an integer.  */)
 340      (n, arg)
 341      register Lisp_Object n, arg;
 342 {
 343   int i1 = extract_float (n);
 344   double f2 = extract_float (arg);
 345
 346   IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
 347   return make_float (f2);
 348 }
 349
 350 DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
 351        doc: /* Return the bessel function y0 of ARG.  */)
 352      (arg)
 353      register Lisp_Object arg;
 354 {
 355   double d = extract_float (arg);
 356   IN_FLOAT (d = y0 (d), "bessel-y0", arg);
 357   return make_float (d);
 358 }
 359
 360 DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
 361        doc: /* Return the bessel function y1 of ARG.  */)
 362      (arg)
 363      register Lisp_Object arg;
 364 {
 365   double d = extract_float (arg);
 366   IN_FLOAT (d = y1 (d), "bessel-y0", arg);
 367   return make_float (d);
 368 }
 369
 370 DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
 371        doc: /* Return the order N bessel function output yn of ARG.
 372 The first arg (the order) is truncated to an integer.  */)
 373      (n, arg)
 374      register Lisp_Object n, arg;
 375 {
 376   int i1 = extract_float (n);
 377   double f2 = extract_float (arg);
 378
 379   IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
 380   return make_float (f2);
 381 }
 382
 383 #endif
 384 \f
 385 #if 0 /* Leave these out unless we see they are worth having.  */
 386
 387 DEFUN ("erf", Ferf, Serf, 1, 1, 0,
 388        doc: /* Return the mathematical error function of ARG.  */)
 389      (arg)
 390      register Lisp_Object arg;
 391 {
 392   double d = extract_float (arg);
 393   IN_FLOAT (d = erf (d), "erf", arg);
 394   return make_float (d);
 395 }
 396
 397 DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
 398        doc: /* Return the complementary error function of ARG.  */)
 399      (arg)
 400      register Lisp_Object arg;
 401 {
 402   double d = extract_float (arg);
 403   IN_FLOAT (d = erfc (d), "erfc", arg);
 404   return make_float (d);
 405 }
 406
 407 DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
 408        doc: /* Return the log gamma of ARG.  */)
 409      (arg)
 410      register Lisp_Object arg;
 411 {
 412   double d = extract_float (arg);
 413   IN_FLOAT (d = lgamma (d), "log-gamma", arg);
 414   return make_float (d);
 415 }
 416
 417 DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
 418        doc: /* Return the cube root of ARG.  */)
 419      (arg)
 420      register Lisp_Object arg;
 421 {
 422   double d = extract_float (arg);
 423 #ifdef HAVE_CBRT
 424   IN_FLOAT (d = cbrt (d), "cube-root", arg);
 425 #else
 426   if (d >= 0.0)
 427     IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
 428   else
 429     IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
 430 #endif
 431   return make_float (d);
 432 }
 433
 434 #endif
 435 \f
 436 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 437        doc: /* Return the exponential base e of ARG.  */)
 438      (arg)
 439      register Lisp_Object arg;
 440 {
 441   double d = extract_float (arg);
 442 #ifdef FLOAT_CHECK_DOMAIN
 443   if (d > 709.7827)   /* Assume IEEE doubles here */
 444     range_error ("exp", arg);
 445   else if (d < -709.0)
 446     return make_float (0.0);
 447   else
 448 #endif
 449     IN_FLOAT (d = exp (d), "exp", arg);
 450   return make_float (d);
 451 }
 452
 453 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 454        doc: /* Return the exponential ARG1 ** ARG2.  */)
 455      (arg1, arg2)
 456      register Lisp_Object arg1, arg2;
 457 {
 458   double f1, f2;
 459
 460   CHECK_NUMBER_OR_FLOAT (arg1);
 461   CHECK_NUMBER_OR_FLOAT (arg2);
 462   if (INTEGERP (arg1)     /* common lisp spec */
 463       && INTEGERP (arg2)) /* don't promote, if both are ints */
 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 second 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 = -(VALMASK >> 1);
 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             Fsignal (Qarith_error, Qnil);
 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         Fsignal (Qarith_error, Qnil);
 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     Fsignal (Qarith_error, Qnil);
 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   in_float = 0;
 985
 986   Fsignal (Qarith_error, Fcons (float_error_arg, Qnil));
 987 }
 988
 989 /* Another idea was to replace the library function `infnan'
 990    where SIGILL is signaled.  */
 991
 992 #endif /* FLOAT_CATCH_SIGILL */
 993
 994 #ifdef HAVE_MATHERR
 995 int
 996 matherr (x)
 997      struct exception *x;
 998 {
 999   Lisp_Object args;
1000   if (! in_float)
1001     /* Not called from emacs-lisp float routines; do the default thing. */
1002     return 0;
1003   if (!strcmp (x->name, "pow"))
1004     x->name = "expt";
1005
1006   args
1007     = Fcons (build_string (x->name),
1008              Fcons (make_float (x->arg1),
1009                     ((!strcmp (x->name, "log") || !strcmp (x->name, "pow"))
1010                      ? Fcons (make_float (x->arg2), Qnil)
1011                      : Qnil)));
1012   switch (x->type)
1013     {
1014     case DOMAIN:        Fsignal (Qdomain_error, args);          break;
1015     case SING:          Fsignal (Qsingularity_error, args);     break;
1016     case OVERFLOW:      Fsignal (Qoverflow_error, args);        break;
1017     case UNDERFLOW:     Fsignal (Qunderflow_error, args);       break;
1018     default:            Fsignal (Qarith_error, args);           break;
1019     }
1020   return (1);   /* don't set errno or print a message */
1021 }
1022 #endif /* HAVE_MATHERR */
1023
1024 void
1025 init_floatfns ()
1026 {
1027 #ifdef FLOAT_CATCH_SIGILL
1028   signal (SIGILL, float_error);
1029 #endif
1030   in_float = 0;
1031 }
1032
1033 void
1034 syms_of_floatfns ()
1035 {
1036   defsubr (&Sacos);
1037   defsubr (&Sasin);
1038   defsubr (&Satan);
1039   defsubr (&Scos);
1040   defsubr (&Ssin);
1041   defsubr (&Stan);
1042 #if 0
1043   defsubr (&Sacosh);
1044   defsubr (&Sasinh);
1045   defsubr (&Satanh);
1046   defsubr (&Scosh);
1047   defsubr (&Ssinh);
1048   defsubr (&Stanh);
1049   defsubr (&Sbessel_y0);
1050   defsubr (&Sbessel_y1);
1051   defsubr (&Sbessel_yn);
1052   defsubr (&Sbessel_j0);
1053   defsubr (&Sbessel_j1);
1054   defsubr (&Sbessel_jn);
1055   defsubr (&Serf);
1056   defsubr (&Serfc);
1057   defsubr (&Slog_gamma);
1058   defsubr (&Scube_root);
1059 #endif
1060   defsubr (&Sfceiling);
1061   defsubr (&Sffloor);
1062   defsubr (&Sfround);
1063   defsubr (&Sftruncate);
1064   defsubr (&Sexp);
1065   defsubr (&Sexpt);
1066   defsubr (&Slog);
1067   defsubr (&Slog10);
1068   defsubr (&Ssqrt);
1069
1070   defsubr (&Sabs);
1071   defsubr (&Sfloat);
1072   defsubr (&Slogb);
1073   defsubr (&Sceiling);
1074   defsubr (&Sfloor);
1075   defsubr (&Sround);
1076   defsubr (&Struncate);
1077 }