git.chise.org Git - chise/xemacs-chise.git.1/blob - src/floatfns.c

   1 /* Primitive operations on floating point for XEmacs Lisp interpreter.
   2    Copyright (C) 1988, 1993, 1994 Free Software Foundation, Inc.
   3
   4 This file is part of XEmacs.
   5
   6 XEmacs is free software; you can redistribute it and/or modify it
   7 under the terms of the GNU General Public License as published by the
   8 Free Software Foundation; either version 2, or (at your option) any
   9 later version.
  10
  11 XEmacs is distributed in the hope that it will be useful, but WITHOUT
  12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with XEmacs; 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 /* Synched up with: FSF 19.30. */
  22
  23 /* ANSI C requires only these float functions:
  24    acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
  25    frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
  26
  27    Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
  28    Define HAVE_CBRT if you have cbrt().
  29    Define HAVE_RINT if you have rint().
  30    If you don't define these, then the appropriate routines will be simulated.
  31
  32    Define HAVE_MATHERR if on a system supporting the SysV matherr() callback.
  33    (This should happen automatically.)
  34
  35    Define FLOAT_CHECK_ERRNO if the float library routines set errno.
  36    This has no effect if HAVE_MATHERR is defined.
  37
  38    Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
  39    (What systems actually do this?  Let me know. -jwz)
  40
  41    Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
  42    either setting errno, or signalling SIGFPE/SIGILL.  Otherwise, domain and
  43    range checking will happen before calling the float routines.  This has
  44    no effect if HAVE_MATHERR is defined (since matherr will be called when
  45    a domain error occurs).
  46  */
  47
  48 #include <config.h>
  49 #include "lisp.h"
  50 #include "syssignal.h"
  51
  52 #ifdef LISP_FLOAT_TYPE
  53
  54 /* Need to define a differentiating symbol -- see sysfloat.h */
  55 #define THIS_FILENAME floatfns
  56 #include "sysfloat.h"
  57
  58 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
  59    if `rint' exists but does not work right.  */
  60 #ifdef HAVE_RINT
  61 #define emacs_rint rint
  62 #else
  63 static double
  64 emacs_rint (double x)
  65 {
  66   double r = floor (x + 0.5);
  67   double diff = fabs (r - x);
  68   /* Round to even and correct for any roundoff errors.  */
  69   if (diff >= 0.5 && (diff > 0.5 || r != 2.0 * floor (r / 2.0)))
  70     r += r < x ? 1.0 : -1.0;
  71   return r;
  72 }
  73 #endif
  74
  75 /* Nonzero while executing in floating point.
  76    This tells float_error what to do.  */
  77 static int in_float;
  78
  79 /* If an argument is out of range for a mathematical function,
  80    here is the actual argument value to use in the error message.  */
  81 static Lisp_Object float_error_arg, float_error_arg2;
  82 static const char *float_error_fn_name;
  83
  84 /* Evaluate the floating point expression D, recording NUM
  85    as the original argument for error messages.
  86    D is normally an assignment expression.
  87    Handle errors which may result in signals or may set errno.
  88
  89    Note that float_error may be declared to return void, so you can't
  90    just cast the zero after the colon to (SIGTYPE) to make the types
  91    check properly.  */
  92 #ifdef FLOAT_CHECK_ERRNO
  93 #define IN_FLOAT(d, name, num)                          \
  94   do {                                                  \
  95     float_error_arg = num;                              \
  96     float_error_fn_name = name;                         \
  97     in_float = 1; errno = 0; (d); in_float = 0;         \
  98     if (errno != 0) in_float_error ();                  \
  99   } while (0)
 100 #define IN_FLOAT2(d, name, num, num2)                   \
 101   do {                                                  \
 102     float_error_arg = num;                              \
 103     float_error_arg2 = num2;                            \
 104     float_error_fn_name = name;                         \
 105     in_float = 2; errno = 0; (d); in_float = 0;         \
 106     if (errno != 0) in_float_error ();                  \
 107   } while (0)
 108 #else
 109 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
 110 #define IN_FLOAT2(d, name, num, num2) (in_float = 2, (d), in_float = 0)
 111 #endif
 112
 113
 114 #define arith_error(op,arg) \
 115   Fsignal (Qarith_error, list2 (build_string (op), arg))
 116 #define range_error(op,arg) \
 117   Fsignal (Qrange_error, list2 (build_string (op), arg))
 118 #define range_error2(op,a1,a2) \
 119   Fsignal (Qrange_error, list3 (build_string (op), a1, a2))
 120 #define domain_error(op,arg) \
 121   Fsignal (Qdomain_error, list2 (build_string (op), arg))
 122 #define domain_error2(op,a1,a2) \
 123   Fsignal (Qdomain_error, list3 (build_string (op), a1, a2))
 124
 125
 126 /* Convert float to Lisp Integer if it fits, else signal a range
 127    error using the given arguments.  */
 128 static Lisp_Object
 129 float_to_int (double x, const char *name, Lisp_Object num, Lisp_Object num2)
 130 {
 131   if (x >= ((EMACS_INT) 1 << (VALBITS-1))
 132       || x <= - ((EMACS_INT) 1 << (VALBITS-1)) - (EMACS_INT) 1)
 133   {
 134     if (!UNBOUNDP (num2))
 135       range_error2 (name, num, num2);
 136     else
 137       range_error (name, num);
 138   }
 139   return (make_int ((EMACS_INT) x));
 140 }
 141
 142
 143 static void
 144 in_float_error (void)
 145 {
 146   switch (errno)
 147   {
 148   case 0:
 149     break;
 150   case EDOM:
 151     if (in_float == 2)
 152       domain_error2 (float_error_fn_name, float_error_arg, float_error_arg2);
 153     else
 154       domain_error (float_error_fn_name, float_error_arg);
 155     break;
 156   case ERANGE:
 157     range_error (float_error_fn_name, float_error_arg);
 158     break;
 159   default:
 160     arith_error (float_error_fn_name, float_error_arg);
 161     break;
 162   }
 163 }
 164
 165 \f
 166 static Lisp_Object
 167 mark_float (Lisp_Object obj)
 168 {
 169   return Qnil;
 170 }
 171
 172 static int
 173 float_equal (Lisp_Object obj1, Lisp_Object obj2, int depth)
 174 {
 175   return (extract_float (obj1) == extract_float (obj2));
 176 }
 177
 178 static unsigned long
 179 float_hash (Lisp_Object obj, int depth)
 180 {
 181   /* mod the value down to 32-bit range */
 182   /* #### change for 64-bit machines */
 183   return (unsigned long) fmod (extract_float (obj), 4e9);
 184 }
 185
 186 static const struct lrecord_description float_description[] = {
 187   { XD_END }
 188 };
 189
 190 DEFINE_BASIC_LRECORD_IMPLEMENTATION ("float", float,
 191                                      mark_float, print_float, 0, float_equal,
 192                                      float_hash, float_description,
 193                                      Lisp_Float);
 194 \f
 195 /* Extract a Lisp number as a `double', or signal an error.  */
 196
 197 double
 198 extract_float (Lisp_Object num)
 199 {
 200   if (FLOATP (num))
 201     return XFLOAT_DATA (num);
 202
 203   if (INTP (num))
 204     return (double) XINT (num);
 205
 206   return extract_float (wrong_type_argument (Qnumberp, num));
 207 }
 208 #endif /* LISP_FLOAT_TYPE */
 209
 210 \f
 211 /* Trig functions.  */
 212 #ifdef LISP_FLOAT_TYPE
 213
 214 DEFUN ("acos", Facos, 1, 1, 0, /*
 215 Return the inverse cosine of NUMBER.
 216 */
 217        (number))
 218 {
 219   double d = extract_float (number);
 220 #ifdef FLOAT_CHECK_DOMAIN
 221   if (d > 1.0 || d < -1.0)
 222     domain_error ("acos", number);
 223 #endif
 224   IN_FLOAT (d = acos (d), "acos", number);
 225   return make_float (d);
 226 }
 227
 228 DEFUN ("asin", Fasin, 1, 1, 0, /*
 229 Return the inverse sine of NUMBER.
 230 */
 231        (number))
 232 {
 233   double d = extract_float (number);
 234 #ifdef FLOAT_CHECK_DOMAIN
 235   if (d > 1.0 || d < -1.0)
 236     domain_error ("asin", number);
 237 #endif
 238   IN_FLOAT (d = asin (d), "asin", number);
 239   return make_float (d);
 240 }
 241
 242 DEFUN ("atan", Fatan, 1, 2, 0, /*
 243 Return the inverse tangent of NUMBER.
 244 If optional second argument NUMBER2 is provided,
 245 return atan2 (NUMBER, NUMBER2).
 246 */
 247        (number, number2))
 248 {
 249   double d = extract_float (number);
 250
 251   if (NILP (number2))
 252     IN_FLOAT (d = atan (d), "atan", number);
 253   else
 254     {
 255       double d2 = extract_float (number2);
 256 #ifdef FLOAT_CHECK_DOMAIN
 257       if (d == 0.0 && d2 == 0.0)
 258         domain_error2 ("atan", number, number2);
 259 #endif
 260       IN_FLOAT2 (d = atan2 (d, d2), "atan", number, number2);
 261     }
 262   return make_float (d);
 263 }
 264
 265 DEFUN ("cos", Fcos, 1, 1, 0, /*
 266 Return the cosine of NUMBER.
 267 */
 268        (number))
 269 {
 270   double d = extract_float (number);
 271   IN_FLOAT (d = cos (d), "cos", number);
 272   return make_float (d);
 273 }
 274
 275 DEFUN ("sin", Fsin, 1, 1, 0, /*
 276 Return the sine of NUMBER.
 277 */
 278        (number))
 279 {
 280   double d = extract_float (number);
 281   IN_FLOAT (d = sin (d), "sin", number);
 282   return make_float (d);
 283 }
 284
 285 DEFUN ("tan", Ftan, 1, 1, 0, /*
 286 Return the tangent of NUMBER.
 287 */
 288        (number))
 289 {
 290   double d = extract_float (number);
 291   double c = cos (d);
 292 #ifdef FLOAT_CHECK_DOMAIN
 293   if (c == 0.0)
 294     domain_error ("tan", number);
 295 #endif
 296   IN_FLOAT (d = (sin (d) / c), "tan", number);
 297   return make_float (d);
 298 }
 299 #endif /* LISP_FLOAT_TYPE (trig functions) */
 300
 301 \f
 302 /* Bessel functions */
 303 #if 0 /* Leave these out unless we find there's a reason for them.  */
 304 /* #ifdef LISP_FLOAT_TYPE */
 305
 306 DEFUN ("bessel-j0", Fbessel_j0, 1, 1, 0, /*
 307 Return the bessel function j0 of NUMBER.
 308 */
 309        (number))
 310 {
 311   double d = extract_float (number);
 312   IN_FLOAT (d = j0 (d), "bessel-j0", number);
 313   return make_float (d);
 314 }
 315
 316 DEFUN ("bessel-j1", Fbessel_j1, 1, 1, 0, /*
 317 Return the bessel function j1 of NUMBER.
 318 */
 319        (number))
 320 {
 321   double d = extract_float (number);
 322   IN_FLOAT (d = j1 (d), "bessel-j1", number);
 323   return make_float (d);
 324 }
 325
 326 DEFUN ("bessel-jn", Fbessel_jn, 2, 2, 0, /*
 327 Return the order N bessel function output jn of NUMBER.
 328 The first number (the order) is truncated to an integer.
 329 */
 330        (number1, number2))
 331 {
 332   int i1 = extract_float (number1);
 333   double f2 = extract_float (number2);
 334
 335   IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", number1);
 336   return make_float (f2);
 337 }
 338
 339 DEFUN ("bessel-y0", Fbessel_y0, 1, 1, 0, /*
 340 Return the bessel function y0 of NUMBER.
 341 */
 342        (number))
 343 {
 344   double d = extract_float (number);
 345   IN_FLOAT (d = y0 (d), "bessel-y0", number);
 346   return make_float (d);
 347 }
 348
 349 DEFUN ("bessel-y1", Fbessel_y1, 1, 1, 0, /*
 350 Return the bessel function y1 of NUMBER.
 351 */
 352        (number))
 353 {
 354   double d = extract_float (number);
 355   IN_FLOAT (d = y1 (d), "bessel-y0", number);
 356   return make_float (d);
 357 }
 358
 359 DEFUN ("bessel-yn", Fbessel_yn, 2, 2, 0, /*
 360 Return the order N bessel function output yn of NUMBER.
 361 The first number (the order) is truncated to an integer.
 362 */
 363        (number1, number2))
 364 {
 365   int i1 = extract_float (number1);
 366   double f2 = extract_float (number2);
 367
 368   IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", number1);
 369   return make_float (f2);
 370 }
 371
 372 #endif /* 0 (bessel functions) */
 373 \f
 374 /* Error functions. */
 375 #if 0 /* Leave these out unless we see they are worth having.  */
 376 /* #ifdef LISP_FLOAT_TYPE */
 377
 378 DEFUN ("erf", Ferf, 1, 1, 0, /*
 379 Return the mathematical error function of NUMBER.
 380 */
 381        (number))
 382 {
 383   double d = extract_float (number);
 384   IN_FLOAT (d = erf (d), "erf", number);
 385   return make_float (d);
 386 }
 387
 388 DEFUN ("erfc", Ferfc, 1, 1, 0, /*
 389 Return the complementary error function of NUMBER.
 390 */
 391        (number))
 392 {
 393   double d = extract_float (number);
 394   IN_FLOAT (d = erfc (d), "erfc", number);
 395   return make_float (d);
 396 }
 397
 398 DEFUN ("log-gamma", Flog_gamma, 1, 1, 0, /*
 399 Return the log gamma of NUMBER.
 400 */
 401        (number))
 402 {
 403   double d = extract_float (number);
 404   IN_FLOAT (d = lgamma (d), "log-gamma", number);
 405   return make_float (d);
 406 }
 407
 408 #endif /* 0 (error functions) */
 409
 410 \f
 411 /* Root and Log functions. */
 412
 413 #ifdef LISP_FLOAT_TYPE
 414 DEFUN ("exp", Fexp, 1, 1, 0, /*
 415 Return the exponential base e of NUMBER.
 416 */
 417        (number))
 418 {
 419   double d = extract_float (number);
 420 #ifdef FLOAT_CHECK_DOMAIN
 421   if (d > 709.7827)   /* Assume IEEE doubles here */
 422     range_error ("exp", number);
 423   else if (d < -709.0)
 424     return make_float (0.0);
 425   else
 426 #endif
 427     IN_FLOAT (d = exp (d), "exp", number);
 428   return make_float (d);
 429 }
 430 #endif /* LISP_FLOAT_TYPE */
 431
 432
 433 DEFUN ("expt", Fexpt, 2, 2, 0, /*
 434 Return the exponential NUMBER1 ** NUMBER2.
 435 */
 436        (number1, number2))
 437 {
 438   if (INTP (number1) && /* common lisp spec */
 439       INTP (number2)) /* don't promote, if both are ints */
 440     {
 441       EMACS_INT retval;
 442       EMACS_INT x = XINT (number1);
 443       EMACS_INT y = XINT (number2);
 444
 445       if (y < 0)
 446         {
 447           if (x == 1)
 448             retval = 1;
 449           else if (x == -1)
 450             retval = (y & 1) ? -1 : 1;
 451           else
 452             retval = 0;
 453         }
 454       else
 455         {
 456           retval = 1;
 457           while (y > 0)
 458             {
 459               if (y & 1)
 460                 retval *= x;
 461               x *= x;
 462               y = (EMACS_UINT) y >> 1;
 463             }
 464         }
 465       return make_int (retval);
 466     }
 467
 468 #ifdef LISP_FLOAT_TYPE
 469   {
 470     double f1 = extract_float (number1);
 471     double f2 = extract_float (number2);
 472     /* Really should check for overflow, too */
 473     if (f1 == 0.0 && f2 == 0.0)
 474       f1 = 1.0;
 475 # ifdef FLOAT_CHECK_DOMAIN
 476     else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor(f2)))
 477       domain_error2 ("expt", number1, number2);
 478 # endif /* FLOAT_CHECK_DOMAIN */
 479     IN_FLOAT2 (f1 = pow (f1, f2), "expt", number1, number2);
 480     return make_float (f1);
 481   }
 482 #else
 483   CHECK_INT_OR_FLOAT (number1);
 484   CHECK_INT_OR_FLOAT (number2);
 485   return Fexpt (number1, number2);
 486 #endif /* LISP_FLOAT_TYPE */
 487 }
 488
 489 #ifdef LISP_FLOAT_TYPE
 490 DEFUN ("log", Flog, 1, 2, 0, /*
 491 Return the natural logarithm of NUMBER.
 492 If second optional argument BASE is given, return the logarithm of
 493 NUMBER using that base.
 494 */
 495        (number, base))
 496 {
 497   double d = extract_float (number);
 498 #ifdef FLOAT_CHECK_DOMAIN
 499   if (d <= 0.0)
 500     domain_error2 ("log", number, base);
 501 #endif
 502   if (NILP (base))
 503     IN_FLOAT (d = log (d), "log", number);
 504   else
 505     {
 506       double b = extract_float (base);
 507 #ifdef FLOAT_CHECK_DOMAIN
 508       if (b <= 0.0 || b == 1.0)
 509         domain_error2 ("log", number, base);
 510 #endif
 511       if (b == 10.0)
 512         IN_FLOAT2 (d = log10 (d), "log", number, base);
 513       else
 514         IN_FLOAT2 (d = (log (d) / log (b)), "log", number, base);
 515     }
 516   return make_float (d);
 517 }
 518
 519
 520 DEFUN ("log10", Flog10, 1, 1, 0, /*
 521 Return the logarithm base 10 of NUMBER.
 522 */
 523        (number))
 524 {
 525   double d = extract_float (number);
 526 #ifdef FLOAT_CHECK_DOMAIN
 527   if (d <= 0.0)
 528     domain_error ("log10", number);
 529 #endif
 530   IN_FLOAT (d = log10 (d), "log10", number);
 531   return make_float (d);
 532 }
 533
 534
 535 DEFUN ("sqrt", Fsqrt, 1, 1, 0, /*
 536 Return the square root of NUMBER.
 537 */
 538        (number))
 539 {
 540   double d = extract_float (number);
 541 #ifdef FLOAT_CHECK_DOMAIN
 542   if (d < 0.0)
 543     domain_error ("sqrt", number);
 544 #endif
 545   IN_FLOAT (d = sqrt (d), "sqrt", number);
 546   return make_float (d);
 547 }
 548
 549
 550 DEFUN ("cube-root", Fcube_root, 1, 1, 0, /*
 551 Return the cube root of NUMBER.
 552 */
 553        (number))
 554 {
 555   double d = extract_float (number);
 556 #ifdef HAVE_CBRT
 557   IN_FLOAT (d = cbrt (d), "cube-root", number);
 558 #else
 559   if (d >= 0.0)
 560     IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", number);
 561   else
 562     IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", number);
 563 #endif
 564   return make_float (d);
 565 }
 566 #endif /* LISP_FLOAT_TYPE */
 567
 568 \f
 569 /* Inverse trig functions. */
 570 #ifdef LISP_FLOAT_TYPE
 571 /* #if 0  Not clearly worth adding...  */
 572
 573 DEFUN ("acosh", Facosh, 1, 1, 0, /*
 574 Return the inverse hyperbolic cosine of NUMBER.
 575 */
 576        (number))
 577 {
 578   double d = extract_float (number);
 579 #ifdef FLOAT_CHECK_DOMAIN
 580   if (d < 1.0)
 581     domain_error ("acosh", number);
 582 #endif
 583 #ifdef HAVE_INVERSE_HYPERBOLIC
 584   IN_FLOAT (d = acosh (d), "acosh", number);
 585 #else
 586   IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", number);
 587 #endif
 588   return make_float (d);
 589 }
 590
 591 DEFUN ("asinh", Fasinh, 1, 1, 0, /*
 592 Return the inverse hyperbolic sine of NUMBER.
 593 */
 594        (number))
 595 {
 596   double d = extract_float (number);
 597 #ifdef HAVE_INVERSE_HYPERBOLIC
 598   IN_FLOAT (d = asinh (d), "asinh", number);
 599 #else
 600   IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", number);
 601 #endif
 602   return make_float (d);
 603 }
 604
 605 DEFUN ("atanh", Fatanh, 1, 1, 0, /*
 606 Return the inverse hyperbolic tangent of NUMBER.
 607 */
 608        (number))
 609 {
 610   double d = extract_float (number);
 611 #ifdef FLOAT_CHECK_DOMAIN
 612   if (d >= 1.0 || d <= -1.0)
 613     domain_error ("atanh", number);
 614 #endif
 615 #ifdef HAVE_INVERSE_HYPERBOLIC
 616   IN_FLOAT (d = atanh (d), "atanh", number);
 617 #else
 618   IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", number);
 619 #endif
 620   return make_float (d);
 621 }
 622
 623 DEFUN ("cosh", Fcosh, 1, 1, 0, /*
 624 Return the hyperbolic cosine of NUMBER.
 625 */
 626        (number))
 627 {
 628   double d = extract_float (number);
 629 #ifdef FLOAT_CHECK_DOMAIN
 630   if (d > 710.0 || d < -710.0)
 631     range_error ("cosh", number);
 632 #endif
 633   IN_FLOAT (d = cosh (d), "cosh", number);
 634   return make_float (d);
 635 }
 636
 637 DEFUN ("sinh", Fsinh, 1, 1, 0, /*
 638 Return the hyperbolic sine of NUMBER.
 639 */
 640        (number))
 641 {
 642   double d = extract_float (number);
 643 #ifdef FLOAT_CHECK_DOMAIN
 644   if (d > 710.0 || d < -710.0)
 645     range_error ("sinh", number);
 646 #endif
 647   IN_FLOAT (d = sinh (d), "sinh", number);
 648   return make_float (d);
 649 }
 650
 651 DEFUN ("tanh", Ftanh, 1, 1, 0, /*
 652 Return the hyperbolic tangent of NUMBER.
 653 */
 654        (number))
 655 {
 656   double d = extract_float (number);
 657   IN_FLOAT (d = tanh (d), "tanh", number);
 658   return make_float (d);
 659 }
 660 #endif /* LISP_FLOAT_TYPE (inverse trig functions) */
 661 \f
 662 /* Rounding functions */
 663
 664 DEFUN ("abs", Fabs, 1, 1, 0, /*
 665 Return the absolute value of NUMBER.
 666 */
 667        (number))
 668 {
 669 #ifdef LISP_FLOAT_TYPE
 670   if (FLOATP (number))
 671     {
 672       IN_FLOAT (number = make_float (fabs (XFLOAT_DATA (number))),
 673                 "abs", number);
 674       return number;
 675     }
 676 #endif /* LISP_FLOAT_TYPE */
 677
 678   if (INTP (number))
 679     return (XINT (number) >= 0) ? number : make_int (- XINT (number));
 680
 681   return Fabs (wrong_type_argument (Qnumberp, number));
 682 }
 683
 684 #ifdef LISP_FLOAT_TYPE
 685 DEFUN ("float", Ffloat, 1, 1, 0, /*
 686 Return the floating point number numerically equal to NUMBER.
 687 */
 688        (number))
 689 {
 690   if (INTP (number))
 691     return make_float ((double) XINT (number));
 692
 693   if (FLOATP (number))          /* give 'em the same float back */
 694     return number;
 695
 696   return Ffloat (wrong_type_argument (Qnumberp, number));
 697 }
 698 #endif /* LISP_FLOAT_TYPE */
 699
 700
 701 #ifdef LISP_FLOAT_TYPE
 702 DEFUN ("logb", Flogb, 1, 1, 0, /*
 703 Return largest integer <= the base 2 log of the magnitude of NUMBER.
 704 This is the same as the exponent of a float.
 705 */
 706        (number))
 707 {
 708   double f = extract_float (number);
 709
 710   if (f == 0.0)
 711     return make_int (- (EMACS_INT)(((EMACS_UINT) 1) << (VALBITS - 1))); /* most-negative-fixnum */
 712 #ifdef HAVE_LOGB
 713   {
 714     Lisp_Object val;
 715     IN_FLOAT (val = make_int ((EMACS_INT) logb (f)), "logb", number);
 716     return val;
 717   }
 718 #else
 719 #ifdef HAVE_FREXP
 720   {
 721     int exqp;
 722     IN_FLOAT (frexp (f, &exqp), "logb", number);
 723     return make_int (exqp - 1);
 724   }
 725 #else
 726   {
 727     int i;
 728     double d;
 729     EMACS_INT val;
 730     if (f < 0.0)
 731       f = -f;
 732     val = -1;
 733     while (f < 0.5)
 734       {
 735         for (i = 1, d = 0.5; d * d >= f; i += i)
 736           d *= d;
 737         f /= d;
 738         val -= i;
 739       }
 740     while (f >= 1.0)
 741       {
 742         for (i = 1, d = 2.0; d * d <= f; i += i)
 743           d *= d;
 744         f /= d;
 745         val += i;
 746       }
 747     return make_int (val);
 748   }
 749 #endif /* ! HAVE_FREXP */
 750 #endif /* ! HAVE_LOGB */
 751 }
 752 #endif /* LISP_FLOAT_TYPE */
 753
 754
 755 DEFUN ("ceiling", Fceiling, 1, 1, 0, /*
 756 Return the smallest integer no less than NUMBER.  (Round toward +inf.)
 757 */
 758        (number))
 759 {
 760 #ifdef LISP_FLOAT_TYPE
 761   if (FLOATP (number))
 762     {
 763       double d;
 764       IN_FLOAT ((d = ceil (XFLOAT_DATA (number))), "ceiling", number);
 765       return (float_to_int (d, "ceiling", number, Qunbound));
 766     }
 767 #endif /* LISP_FLOAT_TYPE */
 768
 769   if (INTP (number))
 770     return number;
 771
 772   return Fceiling (wrong_type_argument (Qnumberp, number));
 773 }
 774
 775
 776 DEFUN ("floor", Ffloor, 1, 2, 0, /*
 777 Return the largest integer no greater than NUMBER.  (Round towards -inf.)
 778 With optional second argument DIVISOR, return the largest integer no
 779 greater than NUMBER/DIVISOR.
 780 */
 781        (number, divisor))
 782 {
 783   CHECK_INT_OR_FLOAT (number);
 784
 785   if (! NILP (divisor))
 786     {
 787       EMACS_INT i1, i2;
 788
 789       CHECK_INT_OR_FLOAT (divisor);
 790
 791 #ifdef LISP_FLOAT_TYPE
 792       if (FLOATP (number) || FLOATP (divisor))
 793         {
 794           double f1 = extract_float (number);
 795           double f2 = extract_float (divisor);
 796
 797           if (f2 == 0)
 798             Fsignal (Qarith_error, Qnil);
 799
 800           IN_FLOAT2 (f1 = floor (f1 / f2), "floor", number, divisor);
 801           return float_to_int (f1, "floor", number, divisor);
 802         }
 803 #endif /* LISP_FLOAT_TYPE */
 804
 805       i1 = XINT (number);
 806       i2 = XINT (divisor);
 807
 808       if (i2 == 0)
 809         Fsignal (Qarith_error, Qnil);
 810
 811       /* With C's /, the result is implementation-defined if either operand
 812          is negative, so use only nonnegative operands.  */
 813       i1 = (i2 < 0
 814             ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
 815             : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
 816
 817       return (make_int (i1));
 818     }
 819
 820 #ifdef LISP_FLOAT_TYPE
 821   if (FLOATP (number))
 822     {
 823       double d;
 824       IN_FLOAT ((d = floor (XFLOAT_DATA (number))), "floor", number);
 825       return (float_to_int (d, "floor", number, Qunbound));
 826     }
 827 #endif /* LISP_FLOAT_TYPE */
 828
 829   return number;
 830 }
 831
 832 DEFUN ("round", Fround, 1, 1, 0, /*
 833 Return the nearest integer to NUMBER.
 834 */
 835        (number))
 836 {
 837 #ifdef LISP_FLOAT_TYPE
 838   if (FLOATP (number))
 839     {
 840       double d;
 841       /* Screw the prevailing rounding mode.  */
 842       IN_FLOAT ((d = emacs_rint (XFLOAT_DATA (number))), "round", number);
 843       return (float_to_int (d, "round", number, Qunbound));
 844     }
 845 #endif /* LISP_FLOAT_TYPE */
 846
 847   if (INTP (number))
 848     return number;
 849
 850   return Fround (wrong_type_argument (Qnumberp, number));
 851 }
 852
 853 DEFUN ("truncate", Ftruncate, 1, 1, 0, /*
 854 Truncate a floating point number to an integer.
 855 Rounds the value toward zero.
 856 */
 857        (number))
 858 {
 859 #ifdef LISP_FLOAT_TYPE
 860   if (FLOATP (number))
 861     return float_to_int (XFLOAT_DATA (number), "truncate", number, Qunbound);
 862 #endif /* LISP_FLOAT_TYPE */
 863
 864   if (INTP (number))
 865     return number;
 866
 867   return Ftruncate (wrong_type_argument (Qnumberp, number));
 868 }
 869 \f
 870 /* Float-rounding functions. */
 871 #ifdef LISP_FLOAT_TYPE
 872 /* #if 1  It's not clear these are worth adding... */
 873
 874 DEFUN ("fceiling", Ffceiling, 1, 1, 0, /*
 875 Return the smallest integer no less than NUMBER, as a float.
 876 \(Round toward +inf.\)
 877 */
 878        (number))
 879 {
 880   double d = extract_float (number);
 881   IN_FLOAT (d = ceil (d), "fceiling", number);
 882   return make_float (d);
 883 }
 884
 885 DEFUN ("ffloor", Fffloor, 1, 1, 0, /*
 886 Return the largest integer no greater than NUMBER, as a float.
 887 \(Round towards -inf.\)
 888 */
 889        (number))
 890 {
 891   double d = extract_float (number);
 892   IN_FLOAT (d = floor (d), "ffloor", number);
 893   return make_float (d);
 894 }
 895
 896 DEFUN ("fround", Ffround, 1, 1, 0, /*
 897 Return the nearest integer to NUMBER, as a float.
 898 */
 899        (number))
 900 {
 901   double d = extract_float (number);
 902   IN_FLOAT (d = emacs_rint (d), "fround", number);
 903   return make_float (d);
 904 }
 905
 906 DEFUN ("ftruncate", Fftruncate, 1, 1, 0, /*
 907 Truncate a floating point number to an integral float value.
 908 Rounds the value toward zero.
 909 */
 910        (number))
 911 {
 912   double d = extract_float (number);
 913   if (d >= 0.0)
 914     IN_FLOAT (d = floor (d), "ftruncate", number);
 915   else
 916     IN_FLOAT (d = ceil (d), "ftruncate", number);
 917   return make_float (d);
 918 }
 919
 920 #endif /* LISP_FLOAT_TYPE (float-rounding functions) */
 921
 922 \f
 923 #ifdef LISP_FLOAT_TYPE
 924 #ifdef FLOAT_CATCH_SIGILL
 925 static SIGTYPE
 926 float_error (int signo)
 927 {
 928   if (! in_float)
 929     fatal_error_signal (signo);
 930
 931   EMACS_REESTABLISH_SIGNAL (signo, arith_error);
 932   EMACS_UNBLOCK_SIGNAL (signo);
 933
 934   in_float = 0;
 935
 936   /* Was Fsignal(), but it just doesn't make sense for an error
 937      occurring inside a signal handler to be restartable, considering
 938      that anything could happen when the error is signaled and trapped
 939      and considering the asynchronous nature of signal handlers. */
 940   signal_error (Qarith_error, list1 (float_error_arg));
 941 }
 942
 943 /* Another idea was to replace the library function `infnan'
 944    where SIGILL is signaled.  */
 945
 946 #endif /* FLOAT_CATCH_SIGILL */
 947
 948 /* In C++, it is impossible to determine what type matherr expects
 949    without some more configure magic.
 950    We shouldn't be using matherr anyways - it's a non-standard SYSVism. */
 951 #if defined (HAVE_MATHERR) && !defined(__cplusplus)
 952 int
 953 matherr (struct exception *x)
 954 {
 955   Lisp_Object args;
 956   if (! in_float)
 957     /* Not called from emacs-lisp float routines; do the default thing. */
 958     return 0;
 959
 960   /* if (!strcmp (x->name, "pow")) x->name = "expt"; */
 961
 962   args = Fcons (build_string (x->name),
 963                 Fcons (make_float (x->arg1),
 964                        ((in_float == 2)
 965                         ? Fcons (make_float (x->arg2), Qnil)
 966                         : Qnil)));
 967   switch (x->type)
 968     {
 969     case DOMAIN:    Fsignal (Qdomain_error,      args); break;
 970     case SING:      Fsignal (Qsingularity_error, args); break;
 971     case OVERFLOW:  Fsignal (Qoverflow_error,    args); break;
 972     case UNDERFLOW: Fsignal (Qunderflow_error,   args); break;
 973     default:        Fsignal (Qarith_error,       args); break;
 974     }
 975   return 1;     /* don't set errno or print a message */
 976 }
 977 #endif /* HAVE_MATHERR */
 978 #endif /* LISP_FLOAT_TYPE */
 979
 980 \f
 981 void
 982 init_floatfns_very_early (void)
 983 {
 984 #ifdef LISP_FLOAT_TYPE
 985 # ifdef FLOAT_CATCH_SIGILL
 986   signal (SIGILL, float_error);
 987 # endif
 988   in_float = 0;
 989 #endif /* LISP_FLOAT_TYPE */
 990 }
 991
 992 void
 993 syms_of_floatfns (void)
 994 {
 995   INIT_LRECORD_IMPLEMENTATION (float);
 996
 997   /* Trig functions.  */
 998
 999 #ifdef LISP_FLOAT_TYPE
1000   DEFSUBR (Facos);
1001   DEFSUBR (Fasin);
1002   DEFSUBR (Fatan);
1003   DEFSUBR (Fcos);
1004   DEFSUBR (Fsin);
1005   DEFSUBR (Ftan);
1006 #endif /* LISP_FLOAT_TYPE */
1007
1008   /* Bessel functions */
1009
1010 #if 0
1011   DEFSUBR (Fbessel_y0);
1012   DEFSUBR (Fbessel_y1);
1013   DEFSUBR (Fbessel_yn);
1014   DEFSUBR (Fbessel_j0);
1015   DEFSUBR (Fbessel_j1);
1016   DEFSUBR (Fbessel_jn);
1017 #endif /* 0 */
1018
1019   /* Error functions. */
1020
1021 #if 0
1022   DEFSUBR (Ferf);
1023   DEFSUBR (Ferfc);
1024   DEFSUBR (Flog_gamma);
1025 #endif /* 0 */
1026
1027   /* Root and Log functions. */
1028
1029 #ifdef LISP_FLOAT_TYPE
1030   DEFSUBR (Fexp);
1031 #endif /* LISP_FLOAT_TYPE */
1032   DEFSUBR (Fexpt);
1033 #ifdef LISP_FLOAT_TYPE
1034   DEFSUBR (Flog);
1035   DEFSUBR (Flog10);
1036   DEFSUBR (Fsqrt);
1037   DEFSUBR (Fcube_root);
1038 #endif /* LISP_FLOAT_TYPE */
1039
1040   /* Inverse trig functions. */
1041
1042 #ifdef LISP_FLOAT_TYPE
1043   DEFSUBR (Facosh);
1044   DEFSUBR (Fasinh);
1045   DEFSUBR (Fatanh);
1046   DEFSUBR (Fcosh);
1047   DEFSUBR (Fsinh);
1048   DEFSUBR (Ftanh);
1049 #endif /* LISP_FLOAT_TYPE */
1050
1051   /* Rounding functions */
1052
1053   DEFSUBR (Fabs);
1054 #ifdef LISP_FLOAT_TYPE
1055   DEFSUBR (Ffloat);
1056   DEFSUBR (Flogb);
1057 #endif /* LISP_FLOAT_TYPE */
1058   DEFSUBR (Fceiling);
1059   DEFSUBR (Ffloor);
1060   DEFSUBR (Fround);
1061   DEFSUBR (Ftruncate);
1062
1063   /* Float-rounding functions. */
1064
1065 #ifdef LISP_FLOAT_TYPE
1066   DEFSUBR (Ffceiling);
1067   DEFSUBR (Fffloor);
1068   DEFSUBR (Ffround);
1069   DEFSUBR (Fftruncate);
1070 #endif /* LISP_FLOAT_TYPE */
1071 }
1072
1073 void
1074 vars_of_floatfns (void)
1075 {
1076 #ifdef LISP_FLOAT_TYPE
1077   Fprovide (intern ("lisp-float-type"));
1078 #endif
1079 }