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 static SIGTYPE float_error ();
 126
 127 /* Nonzero while executing in floating point.
 128    This tells float_error what to do.  */
 129
 130 static int in_float;
 131
 132 /* If an argument is out of range for a mathematical function,
 133    here is the actual argument value to use in the error message.
 134    These variables are used only across the floating point library call
 135    so there is no need to staticpro them.  */
 136
 137 static Lisp_Object float_error_arg, float_error_arg2;
 138
 139 static char *float_error_fn_name;
 140
 141 /* Evaluate the floating point expression D, recording NUM
 142    as the original argument for error messages.
 143    D is normally an assignment expression.
 144    Handle errors which may result in signals or may set errno.
 145
 146    Note that float_error may be declared to return void, so you can't
 147    just cast the zero after the colon to (SIGTYPE) to make the types
 148    check properly.  */
 149
 150 #ifdef FLOAT_CHECK_ERRNO
 151 #define IN_FLOAT(d, name, num)                          \
 152   do {                                                  \
 153     float_error_arg = num;                              \
 154     float_error_fn_name = name;                         \
 155     in_float = 1; errno = 0; (d); in_float = 0;         \
 156     switch (errno) {                                    \
 157     case 0: break;                                      \
 158     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 159     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 160     default:     arith_error (float_error_fn_name, float_error_arg);    \
 161     }                                                   \
 162   } while (0)
 163 #define IN_FLOAT2(d, name, num, num2)                   \
 164   do {                                                  \
 165     float_error_arg = num;                              \
 166     float_error_arg2 = num2;                            \
 167     float_error_fn_name = name;                         \
 168     in_float = 1; errno = 0; (d); in_float = 0;         \
 169     switch (errno) {                                    \
 170     case 0: break;                                      \
 171     case EDOM:   domain_error (float_error_fn_name, float_error_arg);   \
 172     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
 173     default:     arith_error (float_error_fn_name, float_error_arg);    \
 174     }                                                   \
 175   } while (0)
 176 #else
 177 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
 178 #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
 179 #endif
 180
 181 /* Convert float to Lisp_Int if it fits, else signal a range error
 182    using the given arguments.  */
 183 #define FLOAT_TO_INT(x, i, name, num)                                   \
 184   do                                                                    \
 185     {                                                                   \
 186       if ((x) >= (((EMACS_INT) 1) << (VALBITS-1)) ||                    \
 187           (x) <= - (((EMACS_INT) 1) << (VALBITS-1)) - 1)                \
 188         range_error (name, num);                                        \
 189       XSETINT (i,  (EMACS_INT)(x));                                     \
 190     }                                                                   \
 191   while (0)
 192 #define FLOAT_TO_INT2(x, i, name, num1, num2)                           \
 193   do                                                                    \
 194     {                                                                   \
 195       if ((x) >= (((EMACS_INT) 1) << (VALBITS-1)) ||                    \
 196           (x) <= - (((EMACS_INT) 1) << (VALBITS-1)) - 1)                \
 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, 0);
 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   "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   "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, 1, 0,
 259   "Return the inverse tangent of ARG.")
 260   (arg)
 261      register Lisp_Object arg;
 262 {
 263   double d = extract_float (arg);
 264   IN_FLOAT (d = atan (d), "atan", arg);
 265   return make_float (d);
 266 }
 267
 268 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
 269   "Return the cosine of ARG.")
 270   (arg)
 271      register Lisp_Object arg;
 272 {
 273   double d = extract_float (arg);
 274   IN_FLOAT (d = cos (d), "cos", arg);
 275   return make_float (d);
 276 }
 277
 278 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
 279   "Return the sine of ARG.")
 280   (arg)
 281      register Lisp_Object arg;
 282 {
 283   double d = extract_float (arg);
 284   IN_FLOAT (d = sin (d), "sin", arg);
 285   return make_float (d);
 286 }
 287
 288 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
 289   "Return the tangent of ARG.")
 290   (arg)
 291      register Lisp_Object arg;
 292 {
 293   double d = extract_float (arg);
 294   double c = cos (d);
 295 #ifdef FLOAT_CHECK_DOMAIN
 296   if (c == 0.0)
 297     domain_error ("tan", arg);
 298 #endif
 299   IN_FLOAT (d = sin (d) / c, "tan", arg);
 300   return make_float (d);
 301 }
 302 \f
 303 #if 0 /* Leave these out unless we find there's a reason for them.  */
 304
 305 DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
 306   "Return the bessel function j0 of ARG.")
 307   (arg)
 308      register Lisp_Object arg;
 309 {
 310   double d = extract_float (arg);
 311   IN_FLOAT (d = j0 (d), "bessel-j0", arg);
 312   return make_float (d);
 313 }
 314
 315 DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
 316   "Return the bessel function j1 of ARG.")
 317   (arg)
 318      register Lisp_Object arg;
 319 {
 320   double d = extract_float (arg);
 321   IN_FLOAT (d = j1 (d), "bessel-j1", arg);
 322   return make_float (d);
 323 }
 324
 325 DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
 326   "Return the order N bessel function output jn of ARG.\n\
 327 The first arg (the order) is truncated to an integer.")
 328   (n, arg)
 329      register Lisp_Object n, arg;
 330 {
 331   int i1 = extract_float (n);
 332   double f2 = extract_float (arg);
 333
 334   IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
 335   return make_float (f2);
 336 }
 337
 338 DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
 339   "Return the bessel function y0 of ARG.")
 340   (arg)
 341      register Lisp_Object arg;
 342 {
 343   double d = extract_float (arg);
 344   IN_FLOAT (d = y0 (d), "bessel-y0", arg);
 345   return make_float (d);
 346 }
 347
 348 DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
 349   "Return the bessel function y1 of ARG.")
 350   (arg)
 351      register Lisp_Object arg;
 352 {
 353   double d = extract_float (arg);
 354   IN_FLOAT (d = y1 (d), "bessel-y0", arg);
 355   return make_float (d);
 356 }
 357
 358 DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
 359   "Return the order N bessel function output yn of ARG.\n\
 360 The first arg (the order) is truncated to an integer.")
 361   (n, arg)
 362      register Lisp_Object n, arg;
 363 {
 364   int i1 = extract_float (n);
 365   double f2 = extract_float (arg);
 366
 367   IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
 368   return make_float (f2);
 369 }
 370
 371 #endif
 372 \f
 373 #if 0 /* Leave these out unless we see they are worth having.  */
 374
 375 DEFUN ("erf", Ferf, Serf, 1, 1, 0,
 376   "Return the mathematical error function of ARG.")
 377   (arg)
 378      register Lisp_Object arg;
 379 {
 380   double d = extract_float (arg);
 381   IN_FLOAT (d = erf (d), "erf", arg);
 382   return make_float (d);
 383 }
 384
 385 DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
 386   "Return the complementary error function of ARG.")
 387   (arg)
 388      register Lisp_Object arg;
 389 {
 390   double d = extract_float (arg);
 391   IN_FLOAT (d = erfc (d), "erfc", arg);
 392   return make_float (d);
 393 }
 394
 395 DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
 396   "Return the log gamma of ARG.")
 397   (arg)
 398      register Lisp_Object arg;
 399 {
 400   double d = extract_float (arg);
 401   IN_FLOAT (d = lgamma (d), "log-gamma", arg);
 402   return make_float (d);
 403 }
 404
 405 DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
 406   "Return the cube root of ARG.")
 407   (arg)
 408      register Lisp_Object arg;
 409 {
 410   double d = extract_float (arg);
 411 #ifdef HAVE_CBRT
 412   IN_FLOAT (d = cbrt (d), "cube-root", arg);
 413 #else
 414   if (d >= 0.0)
 415     IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
 416   else
 417     IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
 418 #endif
 419   return make_float (d);
 420 }
 421
 422 #endif
 423 \f
 424 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
 425   "Return the exponential base e of ARG.")
 426   (arg)
 427      register Lisp_Object arg;
 428 {
 429   double d = extract_float (arg);
 430 #ifdef FLOAT_CHECK_DOMAIN
 431   if (d > 709.7827)   /* Assume IEEE doubles here */
 432     range_error ("exp", arg);
 433   else if (d < -709.0)
 434     return make_float (0.0);
 435   else
 436 #endif
 437     IN_FLOAT (d = exp (d), "exp", arg);
 438   return make_float (d);
 439 }
 440
 441 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
 442   "Return the exponential ARG1 ** ARG2.")
 443   (arg1, arg2)
 444      register Lisp_Object arg1, arg2;
 445 {
 446   double f1, f2;
 447
 448   CHECK_NUMBER_OR_FLOAT (arg1, 0);
 449   CHECK_NUMBER_OR_FLOAT (arg2, 0);
 450   if (INTEGERP (arg1)     /* common lisp spec */
 451       && INTEGERP (arg2)) /* don't promote, if both are ints */
 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 (f1 = pow (f1, f2), "expt", arg1, arg2);
 492   return make_float (f1);
 493 }
 494
 495 DEFUN ("log", Flog, Slog, 1, 2, 0,
 496   "Return the natural logarithm of ARG.\n\
 497 If second optional argument BASE is given, return log ARG using that base.")
 498   (arg, base)
 499      register Lisp_Object arg, base;
 500 {
 501   double d = extract_float (arg);
 502
 503 #ifdef FLOAT_CHECK_DOMAIN
 504   if (d <= 0.0)
 505     domain_error2 ("log", arg, base);
 506 #endif
 507   if (NILP (base))
 508     IN_FLOAT (d = log (d), "log", arg);
 509   else
 510     {
 511       double b = extract_float (base);
 512
 513 #ifdef FLOAT_CHECK_DOMAIN
 514       if (b <= 0.0 || b == 1.0)
 515         domain_error2 ("log", arg, base);
 516 #endif
 517       if (b == 10.0)
 518         IN_FLOAT2 (d = log10 (d), "log", arg, base);
 519       else
 520         IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
 521     }
 522   return make_float (d);
 523 }
 524
 525 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
 526   "Return the logarithm base 10 of ARG.")
 527   (arg)
 528      register Lisp_Object arg;
 529 {
 530   double d = extract_float (arg);
 531 #ifdef FLOAT_CHECK_DOMAIN
 532   if (d <= 0.0)
 533     domain_error ("log10", arg);
 534 #endif
 535   IN_FLOAT (d = log10 (d), "log10", arg);
 536   return make_float (d);
 537 }
 538
 539 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
 540   "Return the square root of ARG.")
 541   (arg)
 542      register Lisp_Object arg;
 543 {
 544   double d = extract_float (arg);
 545 #ifdef FLOAT_CHECK_DOMAIN
 546   if (d < 0.0)
 547     domain_error ("sqrt", arg);
 548 #endif
 549   IN_FLOAT (d = sqrt (d), "sqrt", arg);
 550   return make_float (d);
 551 }
 552 \f
 553 #if 0 /* Not clearly worth adding.  */
 554
 555 DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
 556   "Return the inverse hyperbolic cosine of ARG.")
 557   (arg)
 558      register Lisp_Object arg;
 559 {
 560   double d = extract_float (arg);
 561 #ifdef FLOAT_CHECK_DOMAIN
 562   if (d < 1.0)
 563     domain_error ("acosh", arg);
 564 #endif
 565 #ifdef HAVE_INVERSE_HYPERBOLIC
 566   IN_FLOAT (d = acosh (d), "acosh", arg);
 567 #else
 568   IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
 569 #endif
 570   return make_float (d);
 571 }
 572
 573 DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
 574   "Return the inverse hyperbolic sine of ARG.")
 575   (arg)
 576      register Lisp_Object arg;
 577 {
 578   double d = extract_float (arg);
 579 #ifdef HAVE_INVERSE_HYPERBOLIC
 580   IN_FLOAT (d = asinh (d), "asinh", arg);
 581 #else
 582   IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
 583 #endif
 584   return make_float (d);
 585 }
 586
 587 DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
 588   "Return the inverse hyperbolic tangent of ARG.")
 589   (arg)
 590      register Lisp_Object arg;
 591 {
 592   double d = extract_float (arg);
 593 #ifdef FLOAT_CHECK_DOMAIN
 594   if (d >= 1.0 || d <= -1.0)
 595     domain_error ("atanh", arg);
 596 #endif
 597 #ifdef HAVE_INVERSE_HYPERBOLIC
 598   IN_FLOAT (d = atanh (d), "atanh", arg);
 599 #else
 600   IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
 601 #endif
 602   return make_float (d);
 603 }
 604
 605 DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
 606   "Return the hyperbolic cosine of ARG.")
 607   (arg)
 608      register Lisp_Object arg;
 609 {
 610   double d = extract_float (arg);
 611 #ifdef FLOAT_CHECK_DOMAIN
 612   if (d > 710.0 || d < -710.0)
 613     range_error ("cosh", arg);
 614 #endif
 615   IN_FLOAT (d = cosh (d), "cosh", arg);
 616   return make_float (d);
 617 }
 618
 619 DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
 620   "Return the hyperbolic sine of ARG.")
 621   (arg)
 622      register Lisp_Object arg;
 623 {
 624   double d = extract_float (arg);
 625 #ifdef FLOAT_CHECK_DOMAIN
 626   if (d > 710.0 || d < -710.0)
 627     range_error ("sinh", arg);
 628 #endif
 629   IN_FLOAT (d = sinh (d), "sinh", arg);
 630   return make_float (d);
 631 }
 632
 633 DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
 634   "Return the hyperbolic tangent of ARG.")
 635   (arg)
 636      register Lisp_Object arg;
 637 {
 638   double d = extract_float (arg);
 639   IN_FLOAT (d = tanh (d), "tanh", arg);
 640   return make_float (d);
 641 }
 642 #endif
 643 \f
 644 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
 645   "Return the absolute value of ARG.")
 646   (arg)
 647      register Lisp_Object arg;
 648 {
 649   CHECK_NUMBER_OR_FLOAT (arg, 0);
 650
 651   if (FLOATP (arg))
 652     IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
 653   else if (XINT (arg) < 0)
 654     XSETINT (arg, - XINT (arg));
 655
 656   return arg;
 657 }
 658
 659 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
 660   "Return the floating point number equal to ARG.")
 661   (arg)
 662      register Lisp_Object arg;
 663 {
 664   CHECK_NUMBER_OR_FLOAT (arg, 0);
 665
 666   if (INTEGERP (arg))
 667     return make_float ((double) XINT (arg));
 668   else                          /* give 'em the same float back */
 669     return arg;
 670 }
 671
 672 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
 673   "Returns largest integer <= the base 2 log of the magnitude of ARG.\n\
 674 This is the same as the exponent of a float.")
 675      (arg)
 676      Lisp_Object arg;
 677 {
 678   Lisp_Object val;
 679   EMACS_INT value;
 680   double f = extract_float (arg);
 681
 682   if (f == 0.0)
 683     value = -(VALMASK >> 1);
 684   else
 685     {
 686 #ifdef HAVE_LOGB
 687       IN_FLOAT (value = logb (f), "logb", arg);
 688 #else
 689 #ifdef HAVE_FREXP
 690       int ivalue;
 691       IN_FLOAT (frexp (f, &ivalue), "logb", arg);
 692       value = ivalue - 1;
 693 #else
 694       int i;
 695       double d;
 696       if (f < 0.0)
 697         f = -f;
 698       value = -1;
 699       while (f < 0.5)
 700         {
 701           for (i = 1, d = 0.5; d * d >= f; i += i)
 702             d *= d;
 703           f /= d;
 704           value -= i;
 705         }
 706       while (f >= 1.0)
 707         {
 708           for (i = 1, d = 2.0; d * d <= f; i += i)
 709             d *= d;
 710           f /= d;
 711           value += i;
 712         }
 713 #endif
 714 #endif
 715     }
 716   XSETINT (val, value);
 717   return val;
 718 }
 719
 720
 721 /* the rounding functions  */
 722
 723 static Lisp_Object
 724 rounding_driver (arg, divisor, double_round, int_round2, name)
 725      register Lisp_Object arg, divisor;
 726      double (*double_round) ();
 727      EMACS_INT (*int_round2) ();
 728      char *name;
 729 {
 730   CHECK_NUMBER_OR_FLOAT (arg, 0);
 731
 732   if (! NILP (divisor))
 733     {
 734       EMACS_INT i1, i2;
 735
 736       CHECK_NUMBER_OR_FLOAT (divisor, 1);
 737
 738       if (FLOATP (arg) || FLOATP (divisor))
 739         {
 740           double f1, f2;
 741
 742           f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
 743           f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
 744           if (! IEEE_FLOATING_POINT && f2 == 0)
 745             Fsignal (Qarith_error, Qnil);
 746
 747           IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
 748           FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
 749           return arg;
 750         }
 751
 752       i1 = XINT (arg);
 753       i2 = XINT (divisor);
 754
 755       if (i2 == 0)
 756         Fsignal (Qarith_error, Qnil);
 757
 758       XSETINT (arg, (*int_round2) (i1, i2));
 759       return arg;
 760     }
 761
 762   if (FLOATP (arg))
 763     {
 764       double d;
 765
 766       IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
 767       FLOAT_TO_INT (d, arg, name, arg);
 768     }
 769
 770   return arg;
 771 }
 772
 773 /* With C's /, the result is implementation-defined if either operand
 774    is negative, so take care with negative operands in the following
 775    integer functions.  */
 776
 777 static EMACS_INT
 778 ceiling2 (i1, i2)
 779      EMACS_INT i1, i2;
 780 {
 781   return (i2 < 0
 782           ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
 783           : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
 784 }
 785
 786 static EMACS_INT
 787 floor2 (i1, i2)
 788      EMACS_INT i1, i2;
 789 {
 790   return (i2 < 0
 791           ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 792           : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 793 }
 794
 795 static EMACS_INT
 796 truncate2 (i1, i2)
 797      EMACS_INT i1, i2;
 798 {
 799   return (i2 < 0
 800           ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
 801           : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
 802 }
 803
 804 static EMACS_INT
 805 round2 (i1, i2)
 806      EMACS_INT i1, i2;
 807 {
 808   /* The C language's division operator gives us one remainder R, but
 809      we want the remainder R1 on the other side of 0 if R1 is closer
 810      to 0 than R is; because we want to round to even, we also want R1
 811      if R and R1 are the same distance from 0 and if C's quotient is
 812      odd.  */
 813   EMACS_INT q = i1 / i2;
 814   EMACS_INT r = i1 % i2;
 815   EMACS_INT abs_r = r < 0 ? -r : r;
 816   EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
 817   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
 818 }
 819
 820 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
 821    if `rint' exists but does not work right.  */
 822 #ifdef HAVE_RINT
 823 #define emacs_rint rint
 824 #else
 825 static double
 826 emacs_rint (d)
 827      double d;
 828 {
 829   return floor (d + 0.5);
 830 }
 831 #endif
 832
 833 static double
 834 double_identity (d)
 835      double d;
 836 {
 837   return d;
 838 }
 839
 840 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
 841   "Return the smallest integer no less than ARG.  (Round toward +inf.)\n\
 842 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.")
 843   (arg, divisor)
 844      Lisp_Object arg, divisor;
 845 {
 846   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
 847 }
 848
 849 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
 850   "Return the largest integer no greater than ARG.  (Round towards -inf.)\n\
 851 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.")
 852   (arg, divisor)
 853      Lisp_Object arg, divisor;
 854 {
 855   return rounding_driver (arg, divisor, floor, floor2, "floor");
 856 }
 857
 858 DEFUN ("round", Fround, Sround, 1, 2, 0,
 859   "Return the nearest integer to ARG.\n\
 860 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.")
 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        "Truncate a floating point number to an int.\n\
 869 Rounds ARG toward zero.\n\
 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     Fsignal (Qarith_error, Qnil);
 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   "Return the smallest integer no less than ARG, as a float.\n\
 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   "Return the largest integer no greater than ARG, as a float.\n\
 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   "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   "Truncate a floating point number to an integral float value.\n\
 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 #ifdef BSD4_1
 956   sigrelse (SIGILL);
 957 #else /* not BSD4_1 */
 958   sigsetmask (SIGEMPTYMASK);
 959 #endif /* not BSD4_1 */
 960 #else
 961   /* Must reestablish handler each time it is called.  */
 962   signal (SIGILL, float_error);
 963 #endif /* BSD_SYSTEM */
 964
 965   in_float = 0;
 966
 967   Fsignal (Qarith_error, Fcons (float_error_arg, Qnil));
 968 }
 969
 970 /* Another idea was to replace the library function `infnan'
 971    where SIGILL is signaled.  */
 972
 973 #endif /* FLOAT_CATCH_SIGILL */
 974
 975 #ifdef HAVE_MATHERR
 976 int
 977 matherr (x)
 978      struct exception *x;
 979 {
 980   Lisp_Object args;
 981   if (! in_float)
 982     /* Not called from emacs-lisp float routines; do the default thing. */
 983     return 0;
 984   if (!strcmp (x->name, "pow"))
 985     x->name = "expt";
 986
 987   args
 988     = Fcons (build_string (x->name),
 989              Fcons (make_float (x->arg1),
 990                     ((!strcmp (x->name, "log") || !strcmp (x->name, "pow"))
 991                      ? Fcons (make_float (x->arg2), Qnil)
 992                      : Qnil)));
 993   switch (x->type)
 994     {
 995     case DOMAIN:        Fsignal (Qdomain_error, args);          break;
 996     case SING:          Fsignal (Qsingularity_error, args);     break;
 997     case OVERFLOW:      Fsignal (Qoverflow_error, args);        break;
 998     case UNDERFLOW:     Fsignal (Qunderflow_error, args);       break;
 999     default:            Fsignal (Qarith_error, args);           break;
1000     }
1001   return (1);   /* don't set errno or print a message */
1002 }
1003 #endif /* HAVE_MATHERR */
1004
1005 void
1006 init_floatfns ()
1007 {
1008 #ifdef FLOAT_CATCH_SIGILL
1009   signal (SIGILL, float_error);
1010 #endif
1011   in_float = 0;
1012 }
1013
1014 void
1015 syms_of_floatfns ()
1016 {
1017   defsubr (&Sacos);
1018   defsubr (&Sasin);
1019   defsubr (&Satan);
1020   defsubr (&Scos);
1021   defsubr (&Ssin);
1022   defsubr (&Stan);
1023 #if 0
1024   defsubr (&Sacosh);
1025   defsubr (&Sasinh);
1026   defsubr (&Satanh);
1027   defsubr (&Scosh);
1028   defsubr (&Ssinh);
1029   defsubr (&Stanh);
1030   defsubr (&Sbessel_y0);
1031   defsubr (&Sbessel_y1);
1032   defsubr (&Sbessel_yn);
1033   defsubr (&Sbessel_j0);
1034   defsubr (&Sbessel_j1);
1035   defsubr (&Sbessel_jn);
1036   defsubr (&Serf);
1037   defsubr (&Serfc);
1038   defsubr (&Slog_gamma);
1039   defsubr (&Scube_root);
1040 #endif
1041   defsubr (&Sfceiling);
1042   defsubr (&Sffloor);
1043   defsubr (&Sfround);
1044   defsubr (&Sftruncate);
1045   defsubr (&Sexp);
1046   defsubr (&Sexpt);
1047   defsubr (&Slog);
1048   defsubr (&Slog10);
1049   defsubr (&Ssqrt);
1050
1051   defsubr (&Sabs);
1052   defsubr (&Sfloat);
1053   defsubr (&Slogb);
1054   defsubr (&Sceiling);
1055   defsubr (&Sfloor);
1056   defsubr (&Sround);
1057   defsubr (&Struncate);
1058 }