os.texi "Startup Summary" small changes
[emacs.git] / src / floatfns.c
blob305c78cae63b4bd7892d20604873698cc25a2cb8
1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
3 Copyright (C) 1988, 1993-1994, 1999, 2001-2012
4 Free Software Foundation, Inc.
6 Author: Wolfgang Rupprecht
7 (according to ack.texi)
9 This file is part of GNU Emacs.
11 GNU Emacs is free software: you can redistribute it and/or modify
12 it under the terms of the GNU General Public License as published by
13 the Free Software Foundation, either version 3 of the License, or
14 (at your option) any later version.
16 GNU Emacs is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 GNU General Public License for more details.
21 You should have received a copy of the GNU General Public License
22 along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>. */
25 /* ANSI C requires only these float functions:
26 acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
27 frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
29 Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
30 Define HAVE_CBRT if you have cbrt.
31 Define HAVE_RINT if you have a working rint.
32 If you don't define these, then the appropriate routines will be simulated.
34 Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
35 (This should happen automatically.)
37 Define FLOAT_CHECK_ERRNO if the float library routines set errno.
38 This has no effect if HAVE_MATHERR is defined.
40 Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
41 (What systems actually do this? Please let us know.)
43 Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
44 either setting errno, or signaling SIGFPE/SIGILL. Otherwise, domain and
45 range checking will happen before calling the float routines. This has
46 no effect if HAVE_MATHERR is defined (since matherr will be called when
47 a domain error occurs.)
50 #include <config.h>
51 #include <signal.h>
52 #include <setjmp.h>
53 #include "lisp.h"
54 #include "syssignal.h"
56 #include <float.h>
57 /* If IEEE_FLOATING_POINT isn't defined, default it from FLT_*. */
58 #ifndef IEEE_FLOATING_POINT
59 #if (FLT_RADIX == 2 && FLT_MANT_DIG == 24 \
60 && FLT_MIN_EXP == -125 && FLT_MAX_EXP == 128)
61 #define IEEE_FLOATING_POINT 1
62 #else
63 #define IEEE_FLOATING_POINT 0
64 #endif
65 #endif
67 #include <math.h>
69 /* This declaration is omitted on some systems, like Ultrix. */
70 #if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
71 extern double logb (double);
72 #endif /* not HPUX and HAVE_LOGB and no logb macro */
74 #if defined (DOMAIN) && defined (SING) && defined (OVERFLOW)
75 /* If those are defined, then this is probably a `matherr' machine. */
76 # ifndef HAVE_MATHERR
77 # define HAVE_MATHERR
78 # endif
79 #endif
81 #ifdef NO_MATHERR
82 #undef HAVE_MATHERR
83 #endif
85 #ifdef HAVE_MATHERR
86 # ifdef FLOAT_CHECK_ERRNO
87 # undef FLOAT_CHECK_ERRNO
88 # endif
89 # ifdef FLOAT_CHECK_DOMAIN
90 # undef FLOAT_CHECK_DOMAIN
91 # endif
92 #endif
94 #ifndef NO_FLOAT_CHECK_ERRNO
95 #define FLOAT_CHECK_ERRNO
96 #endif
98 #ifdef FLOAT_CHECK_ERRNO
99 # include <errno.h>
100 #endif
102 #ifdef FLOAT_CATCH_SIGILL
103 static void float_error ();
104 #endif
106 /* Nonzero while executing in floating point.
107 This tells float_error what to do. */
109 static int in_float;
111 /* If an argument is out of range for a mathematical function,
112 here is the actual argument value to use in the error message.
113 These variables are used only across the floating point library call
114 so there is no need to staticpro them. */
116 static Lisp_Object float_error_arg, float_error_arg2;
118 static const char *float_error_fn_name;
120 /* Evaluate the floating point expression D, recording NUM
121 as the original argument for error messages.
122 D is normally an assignment expression.
123 Handle errors which may result in signals or may set errno.
125 Note that float_error may be declared to return void, so you can't
126 just cast the zero after the colon to (void) to make the types
127 check properly. */
129 #ifdef FLOAT_CHECK_ERRNO
130 #define IN_FLOAT(d, name, num) \
131 do { \
132 float_error_arg = num; \
133 float_error_fn_name = name; \
134 in_float = 1; errno = 0; (d); in_float = 0; \
135 switch (errno) { \
136 case 0: break; \
137 case EDOM: domain_error (float_error_fn_name, float_error_arg); \
138 case ERANGE: range_error (float_error_fn_name, float_error_arg); \
139 default: arith_error (float_error_fn_name, float_error_arg); \
141 } while (0)
142 #define IN_FLOAT2(d, name, num, num2) \
143 do { \
144 float_error_arg = num; \
145 float_error_arg2 = num2; \
146 float_error_fn_name = name; \
147 in_float = 1; errno = 0; (d); in_float = 0; \
148 switch (errno) { \
149 case 0: break; \
150 case EDOM: domain_error (float_error_fn_name, float_error_arg); \
151 case ERANGE: range_error (float_error_fn_name, float_error_arg); \
152 default: arith_error (float_error_fn_name, float_error_arg); \
154 } while (0)
155 #else
156 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
157 #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
158 #endif
160 /* Convert float to Lisp_Int if it fits, else signal a range error
161 using the given arguments. */
162 #define FLOAT_TO_INT(x, i, name, num) \
163 do \
165 if (FIXNUM_OVERFLOW_P (x)) \
166 range_error (name, num); \
167 XSETINT (i, (EMACS_INT)(x)); \
169 while (0)
170 #define FLOAT_TO_INT2(x, i, name, num1, num2) \
171 do \
173 if (FIXNUM_OVERFLOW_P (x)) \
174 range_error2 (name, num1, num2); \
175 XSETINT (i, (EMACS_INT)(x)); \
177 while (0)
179 #define arith_error(op,arg) \
180 xsignal2 (Qarith_error, build_string ((op)), (arg))
181 #define range_error(op,arg) \
182 xsignal2 (Qrange_error, build_string ((op)), (arg))
183 #define range_error2(op,a1,a2) \
184 xsignal3 (Qrange_error, build_string ((op)), (a1), (a2))
185 #define domain_error(op,arg) \
186 xsignal2 (Qdomain_error, build_string ((op)), (arg))
187 #ifdef FLOAT_CHECK_DOMAIN
188 #define domain_error2(op,a1,a2) \
189 xsignal3 (Qdomain_error, build_string ((op)), (a1), (a2))
190 #endif
192 /* Extract a Lisp number as a `double', or signal an error. */
194 double
195 extract_float (Lisp_Object num)
197 CHECK_NUMBER_OR_FLOAT (num);
199 if (FLOATP (num))
200 return XFLOAT_DATA (num);
201 return (double) XINT (num);
204 /* Trig functions. */
206 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
207 doc: /* Return the inverse cosine of ARG. */)
208 (register Lisp_Object arg)
210 double d = extract_float (arg);
211 #ifdef FLOAT_CHECK_DOMAIN
212 if (d > 1.0 || d < -1.0)
213 domain_error ("acos", arg);
214 #endif
215 IN_FLOAT (d = acos (d), "acos", arg);
216 return make_float (d);
219 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
220 doc: /* Return the inverse sine of ARG. */)
221 (register Lisp_Object arg)
223 double d = extract_float (arg);
224 #ifdef FLOAT_CHECK_DOMAIN
225 if (d > 1.0 || d < -1.0)
226 domain_error ("asin", arg);
227 #endif
228 IN_FLOAT (d = asin (d), "asin", arg);
229 return make_float (d);
232 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
233 doc: /* Return the inverse tangent of the arguments.
234 If only one argument Y is given, return the inverse tangent of Y.
235 If two arguments Y and X are given, return the inverse tangent of Y
236 divided by X, i.e. the angle in radians between the vector (X, Y)
237 and the x-axis. */)
238 (register Lisp_Object y, Lisp_Object x)
240 double d = extract_float (y);
242 if (NILP (x))
243 IN_FLOAT (d = atan (d), "atan", y);
244 else
246 double d2 = extract_float (x);
248 IN_FLOAT2 (d = atan2 (d, d2), "atan", y, x);
250 return make_float (d);
253 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
254 doc: /* Return the cosine of ARG. */)
255 (register Lisp_Object arg)
257 double d = extract_float (arg);
258 IN_FLOAT (d = cos (d), "cos", arg);
259 return make_float (d);
262 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
263 doc: /* Return the sine of ARG. */)
264 (register Lisp_Object arg)
266 double d = extract_float (arg);
267 IN_FLOAT (d = sin (d), "sin", arg);
268 return make_float (d);
271 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
272 doc: /* Return the tangent of ARG. */)
273 (register Lisp_Object arg)
275 double d = extract_float (arg);
276 double c = cos (d);
277 #ifdef FLOAT_CHECK_DOMAIN
278 if (c == 0.0)
279 domain_error ("tan", arg);
280 #endif
281 IN_FLOAT (d = sin (d) / c, "tan", arg);
282 return make_float (d);
285 #undef isnan
286 #define isnan(x) ((x) != (x))
288 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
289 doc: /* Return non nil iff argument X is a NaN. */)
290 (Lisp_Object x)
292 CHECK_FLOAT (x);
293 return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
296 #ifdef HAVE_COPYSIGN
297 DEFUN ("copysign", Fcopysign, Scopysign, 2, 2, 0,
298 doc: /* Copy sign of X2 to value of X1, and return the result.
299 Cause an error if X1 or X2 is not a float. */)
300 (Lisp_Object x1, Lisp_Object x2)
302 double f1, f2;
304 CHECK_FLOAT (x1);
305 CHECK_FLOAT (x2);
307 f1 = XFLOAT_DATA (x1);
308 f2 = XFLOAT_DATA (x2);
310 return make_float (copysign (f1, f2));
313 DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
314 doc: /* Get significand and exponent of a floating point number.
315 Breaks the floating point number X into its binary significand SGNFCAND
316 \(a floating point value between 0.5 (included) and 1.0 (excluded))
317 and an integral exponent EXP for 2, such that:
319 X = SGNFCAND * 2^EXP
321 The function returns the cons cell (SGNFCAND . EXP).
322 If X is zero, both parts (SGNFCAND and EXP) are zero. */)
323 (Lisp_Object x)
325 double f = XFLOATINT (x);
327 if (f == 0.0)
328 return Fcons (make_float (0.0), make_number (0));
329 else
331 int exponent;
332 double sgnfcand = frexp (f, &exponent);
333 return Fcons (make_float (sgnfcand), make_number (exponent));
337 DEFUN ("ldexp", Fldexp, Sldexp, 1, 2, 0,
338 doc: /* Construct number X from significand SGNFCAND and exponent EXP.
339 Returns the floating point value resulting from multiplying SGNFCAND
340 (the significand) by 2 raised to the power of EXP (the exponent). */)
341 (Lisp_Object sgnfcand, Lisp_Object exponent)
343 CHECK_NUMBER (exponent);
344 return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exponent)));
346 #endif
348 #if 0 /* Leave these out unless we find there's a reason for them. */
350 DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
351 doc: /* Return the bessel function j0 of ARG. */)
352 (register Lisp_Object arg)
354 double d = extract_float (arg);
355 IN_FLOAT (d = j0 (d), "bessel-j0", arg);
356 return make_float (d);
359 DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
360 doc: /* Return the bessel function j1 of ARG. */)
361 (register Lisp_Object arg)
363 double d = extract_float (arg);
364 IN_FLOAT (d = j1 (d), "bessel-j1", arg);
365 return make_float (d);
368 DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
369 doc: /* Return the order N bessel function output jn of ARG.
370 The first arg (the order) is truncated to an integer. */)
371 (register Lisp_Object n, Lisp_Object arg)
373 int i1 = extract_float (n);
374 double f2 = extract_float (arg);
376 IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
377 return make_float (f2);
380 DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
381 doc: /* Return the bessel function y0 of ARG. */)
382 (register Lisp_Object arg)
384 double d = extract_float (arg);
385 IN_FLOAT (d = y0 (d), "bessel-y0", arg);
386 return make_float (d);
389 DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
390 doc: /* Return the bessel function y1 of ARG. */)
391 (register Lisp_Object arg)
393 double d = extract_float (arg);
394 IN_FLOAT (d = y1 (d), "bessel-y0", arg);
395 return make_float (d);
398 DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
399 doc: /* Return the order N bessel function output yn of ARG.
400 The first arg (the order) is truncated to an integer. */)
401 (register Lisp_Object n, Lisp_Object arg)
403 int i1 = extract_float (n);
404 double f2 = extract_float (arg);
406 IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
407 return make_float (f2);
410 #endif
412 #if 0 /* Leave these out unless we see they are worth having. */
414 DEFUN ("erf", Ferf, Serf, 1, 1, 0,
415 doc: /* Return the mathematical error function of ARG. */)
416 (register Lisp_Object arg)
418 double d = extract_float (arg);
419 IN_FLOAT (d = erf (d), "erf", arg);
420 return make_float (d);
423 DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
424 doc: /* Return the complementary error function of ARG. */)
425 (register Lisp_Object arg)
427 double d = extract_float (arg);
428 IN_FLOAT (d = erfc (d), "erfc", arg);
429 return make_float (d);
432 DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
433 doc: /* Return the log gamma of ARG. */)
434 (register Lisp_Object arg)
436 double d = extract_float (arg);
437 IN_FLOAT (d = lgamma (d), "log-gamma", arg);
438 return make_float (d);
441 DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
442 doc: /* Return the cube root of ARG. */)
443 (register Lisp_Object arg)
445 double d = extract_float (arg);
446 #ifdef HAVE_CBRT
447 IN_FLOAT (d = cbrt (d), "cube-root", arg);
448 #else
449 if (d >= 0.0)
450 IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
451 else
452 IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
453 #endif
454 return make_float (d);
457 #endif
459 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
460 doc: /* Return the exponential base e of ARG. */)
461 (register Lisp_Object arg)
463 double d = extract_float (arg);
464 #ifdef FLOAT_CHECK_DOMAIN
465 if (d > 709.7827) /* Assume IEEE doubles here */
466 range_error ("exp", arg);
467 else if (d < -709.0)
468 return make_float (0.0);
469 else
470 #endif
471 IN_FLOAT (d = exp (d), "exp", arg);
472 return make_float (d);
475 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
476 doc: /* Return the exponential ARG1 ** ARG2. */)
477 (register Lisp_Object arg1, Lisp_Object arg2)
479 double f1, f2, f3;
481 CHECK_NUMBER_OR_FLOAT (arg1);
482 CHECK_NUMBER_OR_FLOAT (arg2);
483 if (INTEGERP (arg1) /* common lisp spec */
484 && INTEGERP (arg2) /* don't promote, if both are ints, and */
485 && 0 <= XINT (arg2)) /* we are sure the result is not fractional */
486 { /* this can be improved by pre-calculating */
487 EMACS_INT acc, x, y; /* some binary powers of x then accumulating */
488 Lisp_Object val;
490 x = XINT (arg1);
491 y = XINT (arg2);
492 acc = 1;
494 if (y < 0)
496 if (x == 1)
497 acc = 1;
498 else if (x == -1)
499 acc = (y & 1) ? -1 : 1;
500 else
501 acc = 0;
503 else
505 while (y > 0)
507 if (y & 1)
508 acc *= x;
509 x *= x;
510 y >>= 1;
513 XSETINT (val, acc);
514 return val;
516 f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
517 f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
518 /* Really should check for overflow, too */
519 if (f1 == 0.0 && f2 == 0.0)
520 f1 = 1.0;
521 #ifdef FLOAT_CHECK_DOMAIN
522 else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor (f2)))
523 domain_error2 ("expt", arg1, arg2);
524 #endif
525 IN_FLOAT2 (f3 = pow (f1, f2), "expt", arg1, arg2);
526 /* Check for overflow in the result. */
527 if (f1 != 0.0 && f3 == 0.0)
528 range_error ("expt", arg1);
529 return make_float (f3);
532 DEFUN ("log", Flog, Slog, 1, 2, 0,
533 doc: /* Return the natural logarithm of ARG.
534 If the optional argument BASE is given, return log ARG using that base. */)
535 (register Lisp_Object arg, Lisp_Object base)
537 double d = extract_float (arg);
539 #ifdef FLOAT_CHECK_DOMAIN
540 if (d <= 0.0)
541 domain_error2 ("log", arg, base);
542 #endif
543 if (NILP (base))
544 IN_FLOAT (d = log (d), "log", arg);
545 else
547 double b = extract_float (base);
549 #ifdef FLOAT_CHECK_DOMAIN
550 if (b <= 0.0 || b == 1.0)
551 domain_error2 ("log", arg, base);
552 #endif
553 if (b == 10.0)
554 IN_FLOAT2 (d = log10 (d), "log", arg, base);
555 else
556 IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
558 return make_float (d);
561 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
562 doc: /* Return the logarithm base 10 of ARG. */)
563 (register Lisp_Object arg)
565 double d = extract_float (arg);
566 #ifdef FLOAT_CHECK_DOMAIN
567 if (d <= 0.0)
568 domain_error ("log10", arg);
569 #endif
570 IN_FLOAT (d = log10 (d), "log10", arg);
571 return make_float (d);
574 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
575 doc: /* Return the square root of ARG. */)
576 (register Lisp_Object arg)
578 double d = extract_float (arg);
579 #ifdef FLOAT_CHECK_DOMAIN
580 if (d < 0.0)
581 domain_error ("sqrt", arg);
582 #endif
583 IN_FLOAT (d = sqrt (d), "sqrt", arg);
584 return make_float (d);
587 #if 0 /* Not clearly worth adding. */
589 DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
590 doc: /* Return the inverse hyperbolic cosine of ARG. */)
591 (register Lisp_Object arg)
593 double d = extract_float (arg);
594 #ifdef FLOAT_CHECK_DOMAIN
595 if (d < 1.0)
596 domain_error ("acosh", arg);
597 #endif
598 #ifdef HAVE_INVERSE_HYPERBOLIC
599 IN_FLOAT (d = acosh (d), "acosh", arg);
600 #else
601 IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
602 #endif
603 return make_float (d);
606 DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
607 doc: /* Return the inverse hyperbolic sine of ARG. */)
608 (register Lisp_Object arg)
610 double d = extract_float (arg);
611 #ifdef HAVE_INVERSE_HYPERBOLIC
612 IN_FLOAT (d = asinh (d), "asinh", arg);
613 #else
614 IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
615 #endif
616 return make_float (d);
619 DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
620 doc: /* Return the inverse hyperbolic tangent of ARG. */)
621 (register Lisp_Object arg)
623 double d = extract_float (arg);
624 #ifdef FLOAT_CHECK_DOMAIN
625 if (d >= 1.0 || d <= -1.0)
626 domain_error ("atanh", arg);
627 #endif
628 #ifdef HAVE_INVERSE_HYPERBOLIC
629 IN_FLOAT (d = atanh (d), "atanh", arg);
630 #else
631 IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
632 #endif
633 return make_float (d);
636 DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
637 doc: /* Return the hyperbolic cosine of ARG. */)
638 (register Lisp_Object arg)
640 double d = extract_float (arg);
641 #ifdef FLOAT_CHECK_DOMAIN
642 if (d > 710.0 || d < -710.0)
643 range_error ("cosh", arg);
644 #endif
645 IN_FLOAT (d = cosh (d), "cosh", arg);
646 return make_float (d);
649 DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
650 doc: /* Return the hyperbolic sine of ARG. */)
651 (register Lisp_Object arg)
653 double d = extract_float (arg);
654 #ifdef FLOAT_CHECK_DOMAIN
655 if (d > 710.0 || d < -710.0)
656 range_error ("sinh", arg);
657 #endif
658 IN_FLOAT (d = sinh (d), "sinh", arg);
659 return make_float (d);
662 DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
663 doc: /* Return the hyperbolic tangent of ARG. */)
664 (register Lisp_Object arg)
666 double d = extract_float (arg);
667 IN_FLOAT (d = tanh (d), "tanh", arg);
668 return make_float (d);
670 #endif
672 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
673 doc: /* Return the absolute value of ARG. */)
674 (register Lisp_Object arg)
676 CHECK_NUMBER_OR_FLOAT (arg);
678 if (FLOATP (arg))
679 IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
680 else if (XINT (arg) < 0)
681 XSETINT (arg, - XINT (arg));
683 return arg;
686 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
687 doc: /* Return the floating point number equal to ARG. */)
688 (register Lisp_Object arg)
690 CHECK_NUMBER_OR_FLOAT (arg);
692 if (INTEGERP (arg))
693 return make_float ((double) XINT (arg));
694 else /* give 'em the same float back */
695 return arg;
698 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
699 doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
700 This is the same as the exponent of a float. */)
701 (Lisp_Object arg)
703 Lisp_Object val;
704 EMACS_INT value;
705 double f = extract_float (arg);
707 if (f == 0.0)
708 value = MOST_NEGATIVE_FIXNUM;
709 else
711 #ifdef HAVE_LOGB
712 IN_FLOAT (value = logb (f), "logb", arg);
713 #else
714 #ifdef HAVE_FREXP
715 int ivalue;
716 IN_FLOAT (frexp (f, &ivalue), "logb", arg);
717 value = ivalue - 1;
718 #else
719 int i;
720 double d;
721 if (f < 0.0)
722 f = -f;
723 value = -1;
724 while (f < 0.5)
726 for (i = 1, d = 0.5; d * d >= f; i += i)
727 d *= d;
728 f /= d;
729 value -= i;
731 while (f >= 1.0)
733 for (i = 1, d = 2.0; d * d <= f; i += i)
734 d *= d;
735 f /= d;
736 value += i;
738 #endif
739 #endif
741 XSETINT (val, value);
742 return val;
746 /* the rounding functions */
748 static Lisp_Object
749 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
750 double (*double_round) (double),
751 EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
752 const char *name)
754 CHECK_NUMBER_OR_FLOAT (arg);
756 if (! NILP (divisor))
758 EMACS_INT i1, i2;
760 CHECK_NUMBER_OR_FLOAT (divisor);
762 if (FLOATP (arg) || FLOATP (divisor))
764 double f1, f2;
766 f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
767 f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
768 if (! IEEE_FLOATING_POINT && f2 == 0)
769 xsignal0 (Qarith_error);
771 IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
772 FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
773 return arg;
776 i1 = XINT (arg);
777 i2 = XINT (divisor);
779 if (i2 == 0)
780 xsignal0 (Qarith_error);
782 XSETINT (arg, (*int_round2) (i1, i2));
783 return arg;
786 if (FLOATP (arg))
788 double d;
790 IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
791 FLOAT_TO_INT (d, arg, name, arg);
794 return arg;
797 /* With C's /, the result is implementation-defined if either operand
798 is negative, so take care with negative operands in the following
799 integer functions. */
801 static EMACS_INT
802 ceiling2 (EMACS_INT i1, EMACS_INT i2)
804 return (i2 < 0
805 ? (i1 < 0 ? ((-1 - i1) / -i2) + 1 : - (i1 / -i2))
806 : (i1 <= 0 ? - (-i1 / i2) : ((i1 - 1) / i2) + 1));
809 static EMACS_INT
810 floor2 (EMACS_INT i1, EMACS_INT i2)
812 return (i2 < 0
813 ? (i1 <= 0 ? -i1 / -i2 : -1 - ((i1 - 1) / -i2))
814 : (i1 < 0 ? -1 - ((-1 - i1) / i2) : i1 / i2));
817 static EMACS_INT
818 truncate2 (EMACS_INT i1, EMACS_INT i2)
820 return (i2 < 0
821 ? (i1 < 0 ? -i1 / -i2 : - (i1 / -i2))
822 : (i1 < 0 ? - (-i1 / i2) : i1 / i2));
825 static EMACS_INT
826 round2 (EMACS_INT i1, EMACS_INT i2)
828 /* The C language's division operator gives us one remainder R, but
829 we want the remainder R1 on the other side of 0 if R1 is closer
830 to 0 than R is; because we want to round to even, we also want R1
831 if R and R1 are the same distance from 0 and if C's quotient is
832 odd. */
833 EMACS_INT q = i1 / i2;
834 EMACS_INT r = i1 % i2;
835 EMACS_INT abs_r = r < 0 ? -r : r;
836 EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
837 return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
840 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
841 if `rint' exists but does not work right. */
842 #ifdef HAVE_RINT
843 #define emacs_rint rint
844 #else
845 static double
846 emacs_rint (double d)
848 return floor (d + 0.5);
850 #endif
852 static double
853 double_identity (double d)
855 return d;
858 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
859 doc: /* Return the smallest integer no less than ARG.
860 This rounds the value towards +inf.
861 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR. */)
862 (Lisp_Object arg, Lisp_Object divisor)
864 return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
867 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
868 doc: /* Return the largest integer no greater than ARG.
869 This rounds the value towards -inf.
870 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR. */)
871 (Lisp_Object arg, Lisp_Object divisor)
873 return rounding_driver (arg, divisor, floor, floor2, "floor");
876 DEFUN ("round", Fround, Sround, 1, 2, 0,
877 doc: /* Return the nearest integer to ARG.
878 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
880 Rounding a value equidistant between two integers may choose the
881 integer closer to zero, or it may prefer an even integer, depending on
882 your machine. For example, \(round 2.5\) can return 3 on some
883 systems, but 2 on others. */)
884 (Lisp_Object arg, Lisp_Object divisor)
886 return rounding_driver (arg, divisor, emacs_rint, round2, "round");
889 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
890 doc: /* Truncate a floating point number to an int.
891 Rounds ARG toward zero.
892 With optional DIVISOR, truncate ARG/DIVISOR. */)
893 (Lisp_Object arg, Lisp_Object divisor)
895 return rounding_driver (arg, divisor, double_identity, truncate2,
896 "truncate");
900 Lisp_Object
901 fmod_float (Lisp_Object x, Lisp_Object y)
903 double f1, f2;
905 f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
906 f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
908 if (! IEEE_FLOATING_POINT && f2 == 0)
909 xsignal0 (Qarith_error);
911 /* If the "remainder" comes out with the wrong sign, fix it. */
912 IN_FLOAT2 ((f1 = fmod (f1, f2),
913 f1 = (f2 < 0 ? f1 > 0 : f1 < 0) ? f1 + f2 : f1),
914 "mod", x, y);
915 return make_float (f1);
918 /* It's not clear these are worth adding. */
920 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
921 doc: /* Return the smallest integer no less than ARG, as a float.
922 \(Round toward +inf.\) */)
923 (register Lisp_Object arg)
925 double d = extract_float (arg);
926 IN_FLOAT (d = ceil (d), "fceiling", arg);
927 return make_float (d);
930 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
931 doc: /* Return the largest integer no greater than ARG, as a float.
932 \(Round towards -inf.\) */)
933 (register Lisp_Object arg)
935 double d = extract_float (arg);
936 IN_FLOAT (d = floor (d), "ffloor", arg);
937 return make_float (d);
940 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
941 doc: /* Return the nearest integer to ARG, as a float. */)
942 (register Lisp_Object arg)
944 double d = extract_float (arg);
945 IN_FLOAT (d = emacs_rint (d), "fround", arg);
946 return make_float (d);
949 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
950 doc: /* Truncate a floating point number to an integral float value.
951 Rounds the value toward zero. */)
952 (register Lisp_Object arg)
954 double d = extract_float (arg);
955 if (d >= 0.0)
956 IN_FLOAT (d = floor (d), "ftruncate", arg);
957 else
958 IN_FLOAT (d = ceil (d), "ftruncate", arg);
959 return make_float (d);
962 #ifdef FLOAT_CATCH_SIGILL
963 static void
964 float_error (int signo)
966 if (! in_float)
967 fatal_error_signal (signo);
969 #ifdef BSD_SYSTEM
970 sigsetmask (SIGEMPTYMASK);
971 #else
972 /* Must reestablish handler each time it is called. */
973 signal (SIGILL, float_error);
974 #endif /* BSD_SYSTEM */
976 SIGNAL_THREAD_CHECK (signo);
977 in_float = 0;
979 xsignal1 (Qarith_error, float_error_arg);
982 /* Another idea was to replace the library function `infnan'
983 where SIGILL is signaled. */
985 #endif /* FLOAT_CATCH_SIGILL */
987 #ifdef HAVE_MATHERR
989 matherr (struct exception *x)
991 Lisp_Object args;
992 const char *name = x->name;
994 if (! in_float)
995 /* Not called from emacs-lisp float routines; do the default thing. */
996 return 0;
997 if (!strcmp (x->name, "pow"))
998 name = "expt";
1000 args
1001 = Fcons (build_string (name),
1002 Fcons (make_float (x->arg1),
1003 ((!strcmp (name, "log") || !strcmp (name, "pow"))
1004 ? Fcons (make_float (x->arg2), Qnil)
1005 : Qnil)));
1006 switch (x->type)
1008 case DOMAIN: xsignal (Qdomain_error, args); break;
1009 case SING: xsignal (Qsingularity_error, args); break;
1010 case OVERFLOW: xsignal (Qoverflow_error, args); break;
1011 case UNDERFLOW: xsignal (Qunderflow_error, args); break;
1012 default: xsignal (Qarith_error, args); break;
1014 return (1); /* don't set errno or print a message */
1016 #endif /* HAVE_MATHERR */
1018 void
1019 init_floatfns (void)
1021 #ifdef FLOAT_CATCH_SIGILL
1022 signal (SIGILL, float_error);
1023 #endif
1024 in_float = 0;
1027 void
1028 syms_of_floatfns (void)
1030 defsubr (&Sacos);
1031 defsubr (&Sasin);
1032 defsubr (&Satan);
1033 defsubr (&Scos);
1034 defsubr (&Ssin);
1035 defsubr (&Stan);
1036 defsubr (&Sisnan);
1037 #ifdef HAVE_COPYSIGN
1038 defsubr (&Scopysign);
1039 defsubr (&Sfrexp);
1040 defsubr (&Sldexp);
1041 #endif
1042 #if 0
1043 defsubr (&Sacosh);
1044 defsubr (&Sasinh);
1045 defsubr (&Satanh);
1046 defsubr (&Scosh);
1047 defsubr (&Ssinh);
1048 defsubr (&Stanh);
1049 defsubr (&Sbessel_y0);
1050 defsubr (&Sbessel_y1);
1051 defsubr (&Sbessel_yn);
1052 defsubr (&Sbessel_j0);
1053 defsubr (&Sbessel_j1);
1054 defsubr (&Sbessel_jn);
1055 defsubr (&Serf);
1056 defsubr (&Serfc);
1057 defsubr (&Slog_gamma);
1058 defsubr (&Scube_root);
1059 #endif
1060 defsubr (&Sfceiling);
1061 defsubr (&Sffloor);
1062 defsubr (&Sfround);
1063 defsubr (&Sftruncate);
1064 defsubr (&Sexp);
1065 defsubr (&Sexpt);
1066 defsubr (&Slog);
1067 defsubr (&Slog10);
1068 defsubr (&Ssqrt);
1070 defsubr (&Sabs);
1071 defsubr (&Sfloat);
1072 defsubr (&Slogb);
1073 defsubr (&Sceiling);
1074 defsubr (&Sfloor);
1075 defsubr (&Sround);
1076 defsubr (&Struncate);