src/floatfns.c

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