*** empty log message ***
[emacs.git] / src / floatfns.c
blobe0ed9b625cd8d39382c27d3295a56bfc1734f9e4
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 Free Software Foundation, Inc.
5 This file is part of GNU Emacs.
7 GNU Emacs is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2, or (at your option)
10 any later version.
12 GNU Emacs is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
17 You should have received a copy of the GNU General Public License
18 along with GNU Emacs; see the file COPYING. If not, write to
19 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
20 Boston, MA 02110-1301, USA. */
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.
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 a working rint.
30 If you don't define these, then the appropriate routines will be simulated.
32 Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
33 (This should happen automatically.)
35 Define FLOAT_CHECK_ERRNO if the float library routines set errno.
36 This has no effect if HAVE_MATHERR is defined.
38 Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
39 (What systems actually do this? Please let us know.)
41 Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
42 either setting errno, or signaling 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.)
48 #include <config.h>
49 #include <signal.h>
50 #include "lisp.h"
51 #include "syssignal.h"
53 #if STDC_HEADERS
54 #include <float.h>
55 #endif
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 /* Work around a problem that happens because math.h on hpux 7
68 defines two static variables--which, in Emacs, are not really static,
69 because `static' is defined as nothing. The problem is that they are
70 defined both here and in lread.c.
71 These macros prevent the name conflict. */
72 #if defined (HPUX) && !defined (HPUX8)
73 #define _MAXLDBL floatfns_maxldbl
74 #define _NMAXLDBL floatfns_nmaxldbl
75 #endif
77 #include <math.h>
79 /* This declaration is omitted on some systems, like Ultrix. */
80 #if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
81 extern double logb ();
82 #endif /* not HPUX and HAVE_LOGB and no logb macro */
84 #if defined(DOMAIN) && defined(SING) && defined(OVERFLOW)
85 /* If those are defined, then this is probably a `matherr' machine. */
86 # ifndef HAVE_MATHERR
87 # define HAVE_MATHERR
88 # endif
89 #endif
91 #ifdef NO_MATHERR
92 #undef HAVE_MATHERR
93 #endif
95 #ifdef HAVE_MATHERR
96 # ifdef FLOAT_CHECK_ERRNO
97 # undef FLOAT_CHECK_ERRNO
98 # endif
99 # ifdef FLOAT_CHECK_DOMAIN
100 # undef FLOAT_CHECK_DOMAIN
101 # endif
102 #endif
104 #ifndef NO_FLOAT_CHECK_ERRNO
105 #define FLOAT_CHECK_ERRNO
106 #endif
108 #ifdef FLOAT_CHECK_ERRNO
109 # include <errno.h>
111 #ifndef USE_CRT_DLL
112 extern int errno;
113 #endif
114 #endif
116 /* Avoid traps on VMS from sinh and cosh.
117 All the other functions set errno instead. */
119 #ifdef VMS
120 #undef cosh
121 #undef sinh
122 #define cosh(x) ((exp(x)+exp(-x))*0.5)
123 #define sinh(x) ((exp(x)-exp(-x))*0.5)
124 #endif /* VMS */
126 #ifdef FLOAT_CATCH_SIGILL
127 static SIGTYPE float_error ();
128 #endif
130 /* Nonzero while executing in floating point.
131 This tells float_error what to do. */
133 static int in_float;
135 /* If an argument is out of range for a mathematical function,
136 here is the actual argument value to use in the error message.
137 These variables are used only across the floating point library call
138 so there is no need to staticpro them. */
140 static Lisp_Object float_error_arg, float_error_arg2;
142 static char *float_error_fn_name;
144 /* Evaluate the floating point expression D, recording NUM
145 as the original argument for error messages.
146 D is normally an assignment expression.
147 Handle errors which may result in signals or may set errno.
149 Note that float_error may be declared to return void, so you can't
150 just cast the zero after the colon to (SIGTYPE) to make the types
151 check properly. */
153 #ifdef FLOAT_CHECK_ERRNO
154 #define IN_FLOAT(d, name, num) \
155 do { \
156 float_error_arg = num; \
157 float_error_fn_name = name; \
158 in_float = 1; errno = 0; (d); in_float = 0; \
159 switch (errno) { \
160 case 0: break; \
161 case EDOM: domain_error (float_error_fn_name, float_error_arg); \
162 case ERANGE: range_error (float_error_fn_name, float_error_arg); \
163 default: arith_error (float_error_fn_name, float_error_arg); \
165 } while (0)
166 #define IN_FLOAT2(d, name, num, num2) \
167 do { \
168 float_error_arg = num; \
169 float_error_arg2 = num2; \
170 float_error_fn_name = name; \
171 in_float = 1; errno = 0; (d); in_float = 0; \
172 switch (errno) { \
173 case 0: break; \
174 case EDOM: domain_error (float_error_fn_name, float_error_arg); \
175 case ERANGE: range_error (float_error_fn_name, float_error_arg); \
176 default: arith_error (float_error_fn_name, float_error_arg); \
178 } while (0)
179 #else
180 #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
181 #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
182 #endif
184 /* Convert float to Lisp_Int if it fits, else signal a range error
185 using the given arguments. */
186 #define FLOAT_TO_INT(x, i, name, num) \
187 do \
189 if (FIXNUM_OVERFLOW_P (x)) \
190 range_error (name, num); \
191 XSETINT (i, (EMACS_INT)(x)); \
193 while (0)
194 #define FLOAT_TO_INT2(x, i, name, num1, num2) \
195 do \
197 if (FIXNUM_OVERFLOW_P (x)) \
198 range_error2 (name, num1, num2); \
199 XSETINT (i, (EMACS_INT)(x)); \
201 while (0)
203 #define arith_error(op,arg) \
204 xsignal2 (Qarith_error, build_string ((op)), (arg))
205 #define range_error(op,arg) \
206 xsignal2 (Qrange_error, build_string ((op)), (arg))
207 #define range_error2(op,a1,a2) \
208 xsignal3 (Qrange_error, build_string ((op)), (a1), (a2))
209 #define domain_error(op,arg) \
210 xsignal2 (Qdomain_error, build_string ((op)), (arg))
211 #define domain_error2(op,a1,a2) \
212 xsignal3 (Qdomain_error, build_string ((op)), (a1), (a2))
214 /* Extract a Lisp number as a `double', or signal an error. */
216 double
217 extract_float (num)
218 Lisp_Object num;
220 CHECK_NUMBER_OR_FLOAT (num);
222 if (FLOATP (num))
223 return XFLOAT_DATA (num);
224 return (double) XINT (num);
227 /* Trig functions. */
229 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
230 doc: /* Return the inverse cosine of ARG. */)
231 (arg)
232 register Lisp_Object arg;
234 double d = extract_float (arg);
235 #ifdef FLOAT_CHECK_DOMAIN
236 if (d > 1.0 || d < -1.0)
237 domain_error ("acos", arg);
238 #endif
239 IN_FLOAT (d = acos (d), "acos", arg);
240 return make_float (d);
243 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
244 doc: /* Return the inverse sine of ARG. */)
245 (arg)
246 register Lisp_Object arg;
248 double d = extract_float (arg);
249 #ifdef FLOAT_CHECK_DOMAIN
250 if (d > 1.0 || d < -1.0)
251 domain_error ("asin", arg);
252 #endif
253 IN_FLOAT (d = asin (d), "asin", arg);
254 return make_float (d);
257 DEFUN ("atan", Fatan, Satan, 1, 2, 0,
258 doc: /* Return the inverse tangent of the arguments.
259 If only one argument Y is given, return the inverse tangent of Y.
260 If two arguments Y and X are given, return the inverse tangent of Y
261 divided by X, i.e. the angle in radians between the vector (X, Y)
262 and the x-axis. */)
263 (y, x)
264 register Lisp_Object y, x;
266 double d = extract_float (y);
268 if (NILP (x))
269 IN_FLOAT (d = atan (d), "atan", y);
270 else
272 double d2 = extract_float (x);
274 IN_FLOAT2 (d = atan2 (d, d2), "atan", y, x);
276 return make_float (d);
279 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
280 doc: /* Return the cosine of ARG. */)
281 (arg)
282 register Lisp_Object arg;
284 double d = extract_float (arg);
285 IN_FLOAT (d = cos (d), "cos", arg);
286 return make_float (d);
289 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
290 doc: /* Return the sine of ARG. */)
291 (arg)
292 register Lisp_Object arg;
294 double d = extract_float (arg);
295 IN_FLOAT (d = sin (d), "sin", arg);
296 return make_float (d);
299 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
300 doc: /* Return the tangent of ARG. */)
301 (arg)
302 register Lisp_Object arg;
304 double d = extract_float (arg);
305 double c = cos (d);
306 #ifdef FLOAT_CHECK_DOMAIN
307 if (c == 0.0)
308 domain_error ("tan", arg);
309 #endif
310 IN_FLOAT (d = sin (d) / c, "tan", arg);
311 return make_float (d);
314 #if 0 /* Leave these out unless we find there's a reason for them. */
316 DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
317 doc: /* Return the bessel function j0 of ARG. */)
318 (arg)
319 register Lisp_Object arg;
321 double d = extract_float (arg);
322 IN_FLOAT (d = j0 (d), "bessel-j0", arg);
323 return make_float (d);
326 DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
327 doc: /* Return the bessel function j1 of ARG. */)
328 (arg)
329 register Lisp_Object arg;
331 double d = extract_float (arg);
332 IN_FLOAT (d = j1 (d), "bessel-j1", arg);
333 return make_float (d);
336 DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
337 doc: /* Return the order N bessel function output jn of ARG.
338 The first arg (the order) is truncated to an integer. */)
339 (n, arg)
340 register Lisp_Object n, arg;
342 int i1 = extract_float (n);
343 double f2 = extract_float (arg);
345 IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
346 return make_float (f2);
349 DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
350 doc: /* Return the bessel function y0 of ARG. */)
351 (arg)
352 register Lisp_Object arg;
354 double d = extract_float (arg);
355 IN_FLOAT (d = y0 (d), "bessel-y0", arg);
356 return make_float (d);
359 DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
360 doc: /* Return the bessel function y1 of ARG. */)
361 (arg)
362 register Lisp_Object arg;
364 double d = extract_float (arg);
365 IN_FLOAT (d = y1 (d), "bessel-y0", arg);
366 return make_float (d);
369 DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
370 doc: /* Return the order N bessel function output yn of ARG.
371 The first arg (the order) is truncated to an integer. */)
372 (n, arg)
373 register Lisp_Object n, arg;
375 int i1 = extract_float (n);
376 double f2 = extract_float (arg);
378 IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
379 return make_float (f2);
382 #endif
384 #if 0 /* Leave these out unless we see they are worth having. */
386 DEFUN ("erf", Ferf, Serf, 1, 1, 0,
387 doc: /* Return the mathematical error function of ARG. */)
388 (arg)
389 register Lisp_Object arg;
391 double d = extract_float (arg);
392 IN_FLOAT (d = erf (d), "erf", arg);
393 return make_float (d);
396 DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
397 doc: /* Return the complementary error function of ARG. */)
398 (arg)
399 register Lisp_Object arg;
401 double d = extract_float (arg);
402 IN_FLOAT (d = erfc (d), "erfc", arg);
403 return make_float (d);
406 DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
407 doc: /* Return the log gamma of ARG. */)
408 (arg)
409 register Lisp_Object arg;
411 double d = extract_float (arg);
412 IN_FLOAT (d = lgamma (d), "log-gamma", arg);
413 return make_float (d);
416 DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
417 doc: /* Return the cube root of ARG. */)
418 (arg)
419 register Lisp_Object arg;
421 double d = extract_float (arg);
422 #ifdef HAVE_CBRT
423 IN_FLOAT (d = cbrt (d), "cube-root", arg);
424 #else
425 if (d >= 0.0)
426 IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
427 else
428 IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
429 #endif
430 return make_float (d);
433 #endif
435 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
436 doc: /* Return the exponential base e of ARG. */)
437 (arg)
438 register Lisp_Object arg;
440 double d = extract_float (arg);
441 #ifdef FLOAT_CHECK_DOMAIN
442 if (d > 709.7827) /* Assume IEEE doubles here */
443 range_error ("exp", arg);
444 else if (d < -709.0)
445 return make_float (0.0);
446 else
447 #endif
448 IN_FLOAT (d = exp (d), "exp", arg);
449 return make_float (d);
452 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
453 doc: /* Return the exponential ARG1 ** ARG2. */)
454 (arg1, arg2)
455 register Lisp_Object arg1, arg2;
457 double f1, f2;
459 CHECK_NUMBER_OR_FLOAT (arg1);
460 CHECK_NUMBER_OR_FLOAT (arg2);
461 if (INTEGERP (arg1) /* common lisp spec */
462 && INTEGERP (arg2) /* don't promote, if both are ints, and */
463 && 0 <= XINT (arg2)) /* we are sure the result is not fractional */
464 { /* this can be improved by pre-calculating */
465 EMACS_INT acc, x, y; /* some binary powers of x then accumulating */
466 Lisp_Object val;
468 x = XINT (arg1);
469 y = XINT (arg2);
470 acc = 1;
472 if (y < 0)
474 if (x == 1)
475 acc = 1;
476 else if (x == -1)
477 acc = (y & 1) ? -1 : 1;
478 else
479 acc = 0;
481 else
483 while (y > 0)
485 if (y & 1)
486 acc *= x;
487 x *= x;
488 y = (unsigned)y >> 1;
491 XSETINT (val, acc);
492 return val;
494 f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
495 f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
496 /* Really should check for overflow, too */
497 if (f1 == 0.0 && f2 == 0.0)
498 f1 = 1.0;
499 #ifdef FLOAT_CHECK_DOMAIN
500 else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor(f2)))
501 domain_error2 ("expt", arg1, arg2);
502 #endif
503 IN_FLOAT2 (f1 = pow (f1, f2), "expt", arg1, arg2);
504 return make_float (f1);
507 DEFUN ("log", Flog, Slog, 1, 2, 0,
508 doc: /* Return the natural logarithm of ARG.
509 If the optional argument BASE is given, return log ARG using that base. */)
510 (arg, base)
511 register Lisp_Object arg, base;
513 double d = extract_float (arg);
515 #ifdef FLOAT_CHECK_DOMAIN
516 if (d <= 0.0)
517 domain_error2 ("log", arg, base);
518 #endif
519 if (NILP (base))
520 IN_FLOAT (d = log (d), "log", arg);
521 else
523 double b = extract_float (base);
525 #ifdef FLOAT_CHECK_DOMAIN
526 if (b <= 0.0 || b == 1.0)
527 domain_error2 ("log", arg, base);
528 #endif
529 if (b == 10.0)
530 IN_FLOAT2 (d = log10 (d), "log", arg, base);
531 else
532 IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
534 return make_float (d);
537 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
538 doc: /* Return the logarithm base 10 of ARG. */)
539 (arg)
540 register Lisp_Object arg;
542 double d = extract_float (arg);
543 #ifdef FLOAT_CHECK_DOMAIN
544 if (d <= 0.0)
545 domain_error ("log10", arg);
546 #endif
547 IN_FLOAT (d = log10 (d), "log10", arg);
548 return make_float (d);
551 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
552 doc: /* Return the square root of ARG. */)
553 (arg)
554 register Lisp_Object arg;
556 double d = extract_float (arg);
557 #ifdef FLOAT_CHECK_DOMAIN
558 if (d < 0.0)
559 domain_error ("sqrt", arg);
560 #endif
561 IN_FLOAT (d = sqrt (d), "sqrt", arg);
562 return make_float (d);
565 #if 0 /* Not clearly worth adding. */
567 DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
568 doc: /* Return the inverse hyperbolic cosine of ARG. */)
569 (arg)
570 register Lisp_Object arg;
572 double d = extract_float (arg);
573 #ifdef FLOAT_CHECK_DOMAIN
574 if (d < 1.0)
575 domain_error ("acosh", arg);
576 #endif
577 #ifdef HAVE_INVERSE_HYPERBOLIC
578 IN_FLOAT (d = acosh (d), "acosh", arg);
579 #else
580 IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
581 #endif
582 return make_float (d);
585 DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
586 doc: /* Return the inverse hyperbolic sine of ARG. */)
587 (arg)
588 register Lisp_Object arg;
590 double d = extract_float (arg);
591 #ifdef HAVE_INVERSE_HYPERBOLIC
592 IN_FLOAT (d = asinh (d), "asinh", arg);
593 #else
594 IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
595 #endif
596 return make_float (d);
599 DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
600 doc: /* Return the inverse hyperbolic tangent of ARG. */)
601 (arg)
602 register Lisp_Object arg;
604 double d = extract_float (arg);
605 #ifdef FLOAT_CHECK_DOMAIN
606 if (d >= 1.0 || d <= -1.0)
607 domain_error ("atanh", arg);
608 #endif
609 #ifdef HAVE_INVERSE_HYPERBOLIC
610 IN_FLOAT (d = atanh (d), "atanh", arg);
611 #else
612 IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
613 #endif
614 return make_float (d);
617 DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
618 doc: /* Return the hyperbolic cosine of ARG. */)
619 (arg)
620 register Lisp_Object arg;
622 double d = extract_float (arg);
623 #ifdef FLOAT_CHECK_DOMAIN
624 if (d > 710.0 || d < -710.0)
625 range_error ("cosh", arg);
626 #endif
627 IN_FLOAT (d = cosh (d), "cosh", arg);
628 return make_float (d);
631 DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
632 doc: /* Return the hyperbolic sine of ARG. */)
633 (arg)
634 register Lisp_Object arg;
636 double d = extract_float (arg);
637 #ifdef FLOAT_CHECK_DOMAIN
638 if (d > 710.0 || d < -710.0)
639 range_error ("sinh", arg);
640 #endif
641 IN_FLOAT (d = sinh (d), "sinh", arg);
642 return make_float (d);
645 DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
646 doc: /* Return the hyperbolic tangent of ARG. */)
647 (arg)
648 register Lisp_Object arg;
650 double d = extract_float (arg);
651 IN_FLOAT (d = tanh (d), "tanh", arg);
652 return make_float (d);
654 #endif
656 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
657 doc: /* Return the absolute value of ARG. */)
658 (arg)
659 register Lisp_Object arg;
661 CHECK_NUMBER_OR_FLOAT (arg);
663 if (FLOATP (arg))
664 IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
665 else if (XINT (arg) < 0)
666 XSETINT (arg, - XINT (arg));
668 return arg;
671 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
672 doc: /* Return the floating point number equal to ARG. */)
673 (arg)
674 register Lisp_Object arg;
676 CHECK_NUMBER_OR_FLOAT (arg);
678 if (INTEGERP (arg))
679 return make_float ((double) XINT (arg));
680 else /* give 'em the same float back */
681 return arg;
684 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
685 doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
686 This is the same as the exponent of a float. */)
687 (arg)
688 Lisp_Object arg;
690 Lisp_Object val;
691 EMACS_INT value;
692 double f = extract_float (arg);
694 if (f == 0.0)
695 value = MOST_NEGATIVE_FIXNUM;
696 else
698 #ifdef HAVE_LOGB
699 IN_FLOAT (value = logb (f), "logb", arg);
700 #else
701 #ifdef HAVE_FREXP
702 int ivalue;
703 IN_FLOAT (frexp (f, &ivalue), "logb", arg);
704 value = ivalue - 1;
705 #else
706 int i;
707 double d;
708 if (f < 0.0)
709 f = -f;
710 value = -1;
711 while (f < 0.5)
713 for (i = 1, d = 0.5; d * d >= f; i += i)
714 d *= d;
715 f /= d;
716 value -= i;
718 while (f >= 1.0)
720 for (i = 1, d = 2.0; d * d <= f; i += i)
721 d *= d;
722 f /= d;
723 value += i;
725 #endif
726 #endif
728 XSETINT (val, value);
729 return val;
733 /* the rounding functions */
735 static Lisp_Object
736 rounding_driver (arg, divisor, double_round, int_round2, name)
737 register Lisp_Object arg, divisor;
738 double (*double_round) ();
739 EMACS_INT (*int_round2) ();
740 char *name;
742 CHECK_NUMBER_OR_FLOAT (arg);
744 if (! NILP (divisor))
746 EMACS_INT i1, i2;
748 CHECK_NUMBER_OR_FLOAT (divisor);
750 if (FLOATP (arg) || FLOATP (divisor))
752 double f1, f2;
754 f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
755 f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
756 if (! IEEE_FLOATING_POINT && f2 == 0)
757 xsignal0 (Qarith_error);
759 IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
760 FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
761 return arg;
764 i1 = XINT (arg);
765 i2 = XINT (divisor);
767 if (i2 == 0)
768 xsignal0 (Qarith_error);
770 XSETINT (arg, (*int_round2) (i1, i2));
771 return arg;
774 if (FLOATP (arg))
776 double d;
778 IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
779 FLOAT_TO_INT (d, arg, name, arg);
782 return arg;
785 /* With C's /, the result is implementation-defined if either operand
786 is negative, so take care with negative operands in the following
787 integer functions. */
789 static EMACS_INT
790 ceiling2 (i1, i2)
791 EMACS_INT i1, i2;
793 return (i2 < 0
794 ? (i1 < 0 ? ((-1 - i1) / -i2) + 1 : - (i1 / -i2))
795 : (i1 <= 0 ? - (-i1 / i2) : ((i1 - 1) / i2) + 1));
798 static EMACS_INT
799 floor2 (i1, i2)
800 EMACS_INT i1, i2;
802 return (i2 < 0
803 ? (i1 <= 0 ? -i1 / -i2 : -1 - ((i1 - 1) / -i2))
804 : (i1 < 0 ? -1 - ((-1 - i1) / i2) : i1 / i2));
807 static EMACS_INT
808 truncate2 (i1, i2)
809 EMACS_INT i1, i2;
811 return (i2 < 0
812 ? (i1 < 0 ? -i1 / -i2 : - (i1 / -i2))
813 : (i1 < 0 ? - (-i1 / i2) : i1 / i2));
816 static EMACS_INT
817 round2 (i1, i2)
818 EMACS_INT i1, i2;
820 /* The C language's division operator gives us one remainder R, but
821 we want the remainder R1 on the other side of 0 if R1 is closer
822 to 0 than R is; because we want to round to even, we also want R1
823 if R and R1 are the same distance from 0 and if C's quotient is
824 odd. */
825 EMACS_INT q = i1 / i2;
826 EMACS_INT r = i1 % i2;
827 EMACS_INT abs_r = r < 0 ? -r : r;
828 EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
829 return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
832 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
833 if `rint' exists but does not work right. */
834 #ifdef HAVE_RINT
835 #define emacs_rint rint
836 #else
837 static double
838 emacs_rint (d)
839 double d;
841 return floor (d + 0.5);
843 #endif
845 static double
846 double_identity (d)
847 double d;
849 return d;
852 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
853 doc: /* Return the smallest integer no less than ARG.
854 This rounds the value towards +inf.
855 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR. */)
856 (arg, divisor)
857 Lisp_Object arg, divisor;
859 return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
862 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
863 doc: /* Return the largest integer no greater than ARG.
864 This rounds the value towards -inf.
865 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR. */)
866 (arg, divisor)
867 Lisp_Object arg, divisor;
869 return rounding_driver (arg, divisor, floor, floor2, "floor");
872 DEFUN ("round", Fround, Sround, 1, 2, 0,
873 doc: /* Return the nearest integer to ARG.
874 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
876 Rounding a value equidistant between two integers may choose the
877 integer closer to zero, or it may prefer an even integer, depending on
878 your machine. For example, \(round 2.5\) can return 3 on some
879 systems, but 2 on others. */)
880 (arg, divisor)
881 Lisp_Object arg, divisor;
883 return rounding_driver (arg, divisor, emacs_rint, round2, "round");
886 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
887 doc: /* Truncate a floating point number to an int.
888 Rounds ARG toward zero.
889 With optional DIVISOR, truncate ARG/DIVISOR. */)
890 (arg, divisor)
891 Lisp_Object arg, divisor;
893 return rounding_driver (arg, divisor, double_identity, truncate2,
894 "truncate");
898 Lisp_Object
899 fmod_float (x, y)
900 register Lisp_Object x, y;
902 double f1, f2;
904 f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
905 f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
907 if (! IEEE_FLOATING_POINT && f2 == 0)
908 xsignal0 (Qarith_error);
910 /* If the "remainder" comes out with the wrong sign, fix it. */
911 IN_FLOAT2 ((f1 = fmod (f1, f2),
912 f1 = (f2 < 0 ? f1 > 0 : f1 < 0) ? f1 + f2 : f1),
913 "mod", x, y);
914 return make_float (f1);
917 /* It's not clear these are worth adding. */
919 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
920 doc: /* Return the smallest integer no less than ARG, as a float.
921 \(Round toward +inf.\) */)
922 (arg)
923 register Lisp_Object arg;
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 (arg)
934 register Lisp_Object arg;
936 double d = extract_float (arg);
937 IN_FLOAT (d = floor (d), "ffloor", arg);
938 return make_float (d);
941 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
942 doc: /* Return the nearest integer to ARG, as a float. */)
943 (arg)
944 register Lisp_Object arg;
946 double d = extract_float (arg);
947 IN_FLOAT (d = emacs_rint (d), "fround", arg);
948 return make_float (d);
951 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
952 doc: /* Truncate a floating point number to an integral float value.
953 Rounds the value toward zero. */)
954 (arg)
955 register Lisp_Object arg;
957 double d = extract_float (arg);
958 if (d >= 0.0)
959 IN_FLOAT (d = floor (d), "ftruncate", arg);
960 else
961 IN_FLOAT (d = ceil (d), "ftruncate", arg);
962 return make_float (d);
965 #ifdef FLOAT_CATCH_SIGILL
966 static SIGTYPE
967 float_error (signo)
968 int signo;
970 if (! in_float)
971 fatal_error_signal (signo);
973 #ifdef BSD_SYSTEM
974 #ifdef BSD4_1
975 sigrelse (SIGILL);
976 #else /* not BSD4_1 */
977 sigsetmask (SIGEMPTYMASK);
978 #endif /* not BSD4_1 */
979 #else
980 /* Must reestablish handler each time it is called. */
981 signal (SIGILL, float_error);
982 #endif /* BSD_SYSTEM */
984 SIGNAL_THREAD_CHECK (signo);
985 in_float = 0;
987 xsignal1 (Qarith_error, float_error_arg);
990 /* Another idea was to replace the library function `infnan'
991 where SIGILL is signaled. */
993 #endif /* FLOAT_CATCH_SIGILL */
995 #ifdef HAVE_MATHERR
997 matherr (x)
998 struct exception *x;
1000 Lisp_Object args;
1001 if (! in_float)
1002 /* Not called from emacs-lisp float routines; do the default thing. */
1003 return 0;
1004 if (!strcmp (x->name, "pow"))
1005 x->name = "expt";
1007 args
1008 = Fcons (build_string (x->name),
1009 Fcons (make_float (x->arg1),
1010 ((!strcmp (x->name, "log") || !strcmp (x->name, "pow"))
1011 ? Fcons (make_float (x->arg2), Qnil)
1012 : Qnil)));
1013 switch (x->type)
1015 case DOMAIN: xsignal (Qdomain_error, args); break;
1016 case SING: xsignal (Qsingularity_error, args); break;
1017 case OVERFLOW: xsignal (Qoverflow_error, args); break;
1018 case UNDERFLOW: xsignal (Qunderflow_error, args); break;
1019 default: xsignal (Qarith_error, args); break;
1021 return (1); /* don't set errno or print a message */
1023 #endif /* HAVE_MATHERR */
1025 void
1026 init_floatfns ()
1028 #ifdef FLOAT_CATCH_SIGILL
1029 signal (SIGILL, float_error);
1030 #endif
1031 in_float = 0;
1034 void
1035 syms_of_floatfns ()
1037 defsubr (&Sacos);
1038 defsubr (&Sasin);
1039 defsubr (&Satan);
1040 defsubr (&Scos);
1041 defsubr (&Ssin);
1042 defsubr (&Stan);
1043 #if 0
1044 defsubr (&Sacosh);
1045 defsubr (&Sasinh);
1046 defsubr (&Satanh);
1047 defsubr (&Scosh);
1048 defsubr (&Ssinh);
1049 defsubr (&Stanh);
1050 defsubr (&Sbessel_y0);
1051 defsubr (&Sbessel_y1);
1052 defsubr (&Sbessel_yn);
1053 defsubr (&Sbessel_j0);
1054 defsubr (&Sbessel_j1);
1055 defsubr (&Sbessel_jn);
1056 defsubr (&Serf);
1057 defsubr (&Serfc);
1058 defsubr (&Slog_gamma);
1059 defsubr (&Scube_root);
1060 #endif
1061 defsubr (&Sfceiling);
1062 defsubr (&Sffloor);
1063 defsubr (&Sfround);
1064 defsubr (&Sftruncate);
1065 defsubr (&Sexp);
1066 defsubr (&Sexpt);
1067 defsubr (&Slog);
1068 defsubr (&Slog10);
1069 defsubr (&Ssqrt);
1071 defsubr (&Sabs);
1072 defsubr (&Sfloat);
1073 defsubr (&Slogb);
1074 defsubr (&Sceiling);
1075 defsubr (&Sfloor);
1076 defsubr (&Sround);
1077 defsubr (&Struncate);
1080 /* arch-tag: be05bf9d-049e-4e31-91b9-e6153d483ae7
1081 (do not change this comment) */