Merge from trunk.
[bpt/emacs.git] / src / floatfns.c
1 /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
2
3 Copyright (C) 1988, 1993-1994, 1999, 2001-2011
4 Free Software Foundation, Inc.
5
6 Author: Wolfgang Rupprecht
7 (according to ack.texi)
8
9 This file is part of GNU Emacs.
10
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.
15
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.
20
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/>. */
23
24
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.
28
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.
33
34 Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
35 (This should happen automatically.)
36
37 Define FLOAT_CHECK_ERRNO if the float library routines set errno.
38 This has no effect if HAVE_MATHERR is defined.
39
40 Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
41 (What systems actually do this? Please let us know.)
42
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.)
48 */
49
50 #include <config.h>
51 #include <signal.h>
52 #include <setjmp.h>
53 #include "lisp.h"
54 #include "syssignal.h"
55
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
66
67 #include <math.h>
68
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 */
73
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
80
81 #ifdef NO_MATHERR
82 #undef HAVE_MATHERR
83 #endif
84
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
93
94 #ifndef NO_FLOAT_CHECK_ERRNO
95 #define FLOAT_CHECK_ERRNO
96 #endif
97
98 #ifdef FLOAT_CHECK_ERRNO
99 # include <errno.h>
100 #endif
101
102 #ifdef FLOAT_CATCH_SIGILL
103 static void float_error ();
104 #endif
105
106 /* Nonzero while executing in floating point.
107 This tells float_error what to do. */
108
109 static int in_float;
110
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. */
115
116 static Lisp_Object float_error_arg, float_error_arg2;
117
118 static const char *float_error_fn_name;
119
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.
124
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. */
128
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); \
140 } \
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); \
153 } \
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
159
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 \
164 { \
165 if (FIXNUM_OVERFLOW_P (x)) \
166 range_error (name, num); \
167 XSETINT (i, (EMACS_INT)(x)); \
168 } \
169 while (0)
170 #define FLOAT_TO_INT2(x, i, name, num1, num2) \
171 do \
172 { \
173 if (FIXNUM_OVERFLOW_P (x)) \
174 range_error2 (name, num1, num2); \
175 XSETINT (i, (EMACS_INT)(x)); \
176 } \
177 while (0)
178
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
191
192 /* Extract a Lisp number as a `double', or signal an error. */
193
194 double
195 extract_float (Lisp_Object num)
196 {
197 CHECK_NUMBER_OR_FLOAT (num);
198
199 if (FLOATP (num))
200 return XFLOAT_DATA (num);
201 return (double) XINT (num);
202 }
203 \f
204 /* Trig functions. */
205
206 DEFUN ("acos", Facos, Sacos, 1, 1, 0,
207 doc: /* Return the inverse cosine of ARG. */)
208 (register Lisp_Object arg)
209 {
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);
217 }
218
219 DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
220 doc: /* Return the inverse sine of ARG. */)
221 (register Lisp_Object arg)
222 {
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);
230 }
231
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)
239 {
240 double d = extract_float (y);
241
242 if (NILP (x))
243 IN_FLOAT (d = atan (d), "atan", y);
244 else
245 {
246 double d2 = extract_float (x);
247
248 IN_FLOAT2 (d = atan2 (d, d2), "atan", y, x);
249 }
250 return make_float (d);
251 }
252
253 DEFUN ("cos", Fcos, Scos, 1, 1, 0,
254 doc: /* Return the cosine of ARG. */)
255 (register Lisp_Object arg)
256 {
257 double d = extract_float (arg);
258 IN_FLOAT (d = cos (d), "cos", arg);
259 return make_float (d);
260 }
261
262 DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
263 doc: /* Return the sine of ARG. */)
264 (register Lisp_Object arg)
265 {
266 double d = extract_float (arg);
267 IN_FLOAT (d = sin (d), "sin", arg);
268 return make_float (d);
269 }
270
271 DEFUN ("tan", Ftan, Stan, 1, 1, 0,
272 doc: /* Return the tangent of ARG. */)
273 (register Lisp_Object arg)
274 {
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);
283 }
284
285 #undef isnan
286 #define isnan(x) ((x) != (x))
287
288 DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
289 doc: /* Return non nil iff argument X is a NaN. */)
290 (Lisp_Object x)
291 {
292 CHECK_FLOAT (x);
293 return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
294 }
295
296 #ifdef HAVE_COPYSIGN
297 DEFUN ("copysign", Fcopysign, Scopysign, 1, 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)
301 {
302 double f1, f2;
303
304 CHECK_FLOAT (x1);
305 CHECK_FLOAT (x2);
306
307 f1 = XFLOAT_DATA (x1);
308 f2 = XFLOAT_DATA (x2);
309
310 return make_float (copysign (f1, f2));
311 }
312
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:
318
319 X = SGNFCAND * 2^EXP
320
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)
324 {
325 double f = XFLOATINT (x);
326
327 if (f == 0.0)
328 return Fcons (make_float (0.0), make_number (0));
329 else
330 {
331 int exponent;
332 double sgnfcand = frexp (f, &exponent);
333 return Fcons (make_float (sgnfcand), make_number (exponent));
334 }
335 }
336
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)
342 {
343 CHECK_NUMBER (exponent);
344 return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exponent)));
345 }
346 #endif
347 \f
348 #if 0 /* Leave these out unless we find there's a reason for them. */
349
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)
353 {
354 double d = extract_float (arg);
355 IN_FLOAT (d = j0 (d), "bessel-j0", arg);
356 return make_float (d);
357 }
358
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)
362 {
363 double d = extract_float (arg);
364 IN_FLOAT (d = j1 (d), "bessel-j1", arg);
365 return make_float (d);
366 }
367
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)
372 {
373 int i1 = extract_float (n);
374 double f2 = extract_float (arg);
375
376 IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
377 return make_float (f2);
378 }
379
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)
383 {
384 double d = extract_float (arg);
385 IN_FLOAT (d = y0 (d), "bessel-y0", arg);
386 return make_float (d);
387 }
388
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)
392 {
393 double d = extract_float (arg);
394 IN_FLOAT (d = y1 (d), "bessel-y0", arg);
395 return make_float (d);
396 }
397
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)
402 {
403 int i1 = extract_float (n);
404 double f2 = extract_float (arg);
405
406 IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
407 return make_float (f2);
408 }
409
410 #endif
411 \f
412 #if 0 /* Leave these out unless we see they are worth having. */
413
414 DEFUN ("erf", Ferf, Serf, 1, 1, 0,
415 doc: /* Return the mathematical error function of ARG. */)
416 (register Lisp_Object arg)
417 {
418 double d = extract_float (arg);
419 IN_FLOAT (d = erf (d), "erf", arg);
420 return make_float (d);
421 }
422
423 DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
424 doc: /* Return the complementary error function of ARG. */)
425 (register Lisp_Object arg)
426 {
427 double d = extract_float (arg);
428 IN_FLOAT (d = erfc (d), "erfc", arg);
429 return make_float (d);
430 }
431
432 DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
433 doc: /* Return the log gamma of ARG. */)
434 (register Lisp_Object arg)
435 {
436 double d = extract_float (arg);
437 IN_FLOAT (d = lgamma (d), "log-gamma", arg);
438 return make_float (d);
439 }
440
441 DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
442 doc: /* Return the cube root of ARG. */)
443 (register Lisp_Object arg)
444 {
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);
455 }
456
457 #endif
458 \f
459 DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
460 doc: /* Return the exponential base e of ARG. */)
461 (register Lisp_Object arg)
462 {
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);
473 }
474
475 DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
476 doc: /* Return the exponential ARG1 ** ARG2. */)
477 (register Lisp_Object arg1, Lisp_Object arg2)
478 {
479 double f1, f2, f3;
480
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 y; /* some binary powers of x then accumulating */
488 EMACS_UINT acc, x; /* Unsigned so that overflow is well defined. */
489 Lisp_Object val;
490
491 x = XINT (arg1);
492 y = XINT (arg2);
493 acc = 1;
494
495 if (y < 0)
496 {
497 if (x == 1)
498 acc = 1;
499 else if (x == -1)
500 acc = (y & 1) ? -1 : 1;
501 else
502 acc = 0;
503 }
504 else
505 {
506 while (y > 0)
507 {
508 if (y & 1)
509 acc *= x;
510 x *= x;
511 y >>= 1;
512 }
513 }
514 XSETINT (val, acc);
515 return val;
516 }
517 f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
518 f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
519 /* Really should check for overflow, too */
520 if (f1 == 0.0 && f2 == 0.0)
521 f1 = 1.0;
522 #ifdef FLOAT_CHECK_DOMAIN
523 else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor (f2)))
524 domain_error2 ("expt", arg1, arg2);
525 #endif
526 IN_FLOAT2 (f3 = pow (f1, f2), "expt", arg1, arg2);
527 /* Check for overflow in the result. */
528 if (f1 != 0.0 && f3 == 0.0)
529 range_error ("expt", arg1);
530 return make_float (f3);
531 }
532
533 DEFUN ("log", Flog, Slog, 1, 2, 0,
534 doc: /* Return the natural logarithm of ARG.
535 If the optional argument BASE is given, return log ARG using that base. */)
536 (register Lisp_Object arg, Lisp_Object base)
537 {
538 double d = extract_float (arg);
539
540 #ifdef FLOAT_CHECK_DOMAIN
541 if (d <= 0.0)
542 domain_error2 ("log", arg, base);
543 #endif
544 if (NILP (base))
545 IN_FLOAT (d = log (d), "log", arg);
546 else
547 {
548 double b = extract_float (base);
549
550 #ifdef FLOAT_CHECK_DOMAIN
551 if (b <= 0.0 || b == 1.0)
552 domain_error2 ("log", arg, base);
553 #endif
554 if (b == 10.0)
555 IN_FLOAT2 (d = log10 (d), "log", arg, base);
556 else
557 IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
558 }
559 return make_float (d);
560 }
561
562 DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
563 doc: /* Return the logarithm base 10 of ARG. */)
564 (register Lisp_Object arg)
565 {
566 double d = extract_float (arg);
567 #ifdef FLOAT_CHECK_DOMAIN
568 if (d <= 0.0)
569 domain_error ("log10", arg);
570 #endif
571 IN_FLOAT (d = log10 (d), "log10", arg);
572 return make_float (d);
573 }
574
575 DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
576 doc: /* Return the square root of ARG. */)
577 (register Lisp_Object arg)
578 {
579 double d = extract_float (arg);
580 #ifdef FLOAT_CHECK_DOMAIN
581 if (d < 0.0)
582 domain_error ("sqrt", arg);
583 #endif
584 IN_FLOAT (d = sqrt (d), "sqrt", arg);
585 return make_float (d);
586 }
587 \f
588 #if 0 /* Not clearly worth adding. */
589
590 DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
591 doc: /* Return the inverse hyperbolic cosine of ARG. */)
592 (register Lisp_Object arg)
593 {
594 double d = extract_float (arg);
595 #ifdef FLOAT_CHECK_DOMAIN
596 if (d < 1.0)
597 domain_error ("acosh", arg);
598 #endif
599 #ifdef HAVE_INVERSE_HYPERBOLIC
600 IN_FLOAT (d = acosh (d), "acosh", arg);
601 #else
602 IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
603 #endif
604 return make_float (d);
605 }
606
607 DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
608 doc: /* Return the inverse hyperbolic sine of ARG. */)
609 (register Lisp_Object arg)
610 {
611 double d = extract_float (arg);
612 #ifdef HAVE_INVERSE_HYPERBOLIC
613 IN_FLOAT (d = asinh (d), "asinh", arg);
614 #else
615 IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
616 #endif
617 return make_float (d);
618 }
619
620 DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
621 doc: /* Return the inverse hyperbolic tangent of ARG. */)
622 (register Lisp_Object arg)
623 {
624 double d = extract_float (arg);
625 #ifdef FLOAT_CHECK_DOMAIN
626 if (d >= 1.0 || d <= -1.0)
627 domain_error ("atanh", arg);
628 #endif
629 #ifdef HAVE_INVERSE_HYPERBOLIC
630 IN_FLOAT (d = atanh (d), "atanh", arg);
631 #else
632 IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
633 #endif
634 return make_float (d);
635 }
636
637 DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
638 doc: /* Return the hyperbolic cosine of ARG. */)
639 (register Lisp_Object arg)
640 {
641 double d = extract_float (arg);
642 #ifdef FLOAT_CHECK_DOMAIN
643 if (d > 710.0 || d < -710.0)
644 range_error ("cosh", arg);
645 #endif
646 IN_FLOAT (d = cosh (d), "cosh", arg);
647 return make_float (d);
648 }
649
650 DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
651 doc: /* Return the hyperbolic sine of ARG. */)
652 (register Lisp_Object arg)
653 {
654 double d = extract_float (arg);
655 #ifdef FLOAT_CHECK_DOMAIN
656 if (d > 710.0 || d < -710.0)
657 range_error ("sinh", arg);
658 #endif
659 IN_FLOAT (d = sinh (d), "sinh", arg);
660 return make_float (d);
661 }
662
663 DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
664 doc: /* Return the hyperbolic tangent of ARG. */)
665 (register Lisp_Object arg)
666 {
667 double d = extract_float (arg);
668 IN_FLOAT (d = tanh (d), "tanh", arg);
669 return make_float (d);
670 }
671 #endif
672 \f
673 DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
674 doc: /* Return the absolute value of ARG. */)
675 (register Lisp_Object arg)
676 {
677 CHECK_NUMBER_OR_FLOAT (arg);
678
679 if (FLOATP (arg))
680 IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
681 else if (XINT (arg) < 0)
682 XSETINT (arg, - XINT (arg));
683
684 return arg;
685 }
686
687 DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
688 doc: /* Return the floating point number equal to ARG. */)
689 (register Lisp_Object arg)
690 {
691 CHECK_NUMBER_OR_FLOAT (arg);
692
693 if (INTEGERP (arg))
694 return make_float ((double) XINT (arg));
695 else /* give 'em the same float back */
696 return arg;
697 }
698
699 DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
700 doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
701 This is the same as the exponent of a float. */)
702 (Lisp_Object arg)
703 {
704 Lisp_Object val;
705 EMACS_INT value;
706 double f = extract_float (arg);
707
708 if (f == 0.0)
709 value = MOST_NEGATIVE_FIXNUM;
710 else
711 {
712 #ifdef HAVE_LOGB
713 IN_FLOAT (value = logb (f), "logb", arg);
714 #else
715 #ifdef HAVE_FREXP
716 int ivalue;
717 IN_FLOAT (frexp (f, &ivalue), "logb", arg);
718 value = ivalue - 1;
719 #else
720 int i;
721 double d;
722 if (f < 0.0)
723 f = -f;
724 value = -1;
725 while (f < 0.5)
726 {
727 for (i = 1, d = 0.5; d * d >= f; i += i)
728 d *= d;
729 f /= d;
730 value -= i;
731 }
732 while (f >= 1.0)
733 {
734 for (i = 1, d = 2.0; d * d <= f; i += i)
735 d *= d;
736 f /= d;
737 value += i;
738 }
739 #endif
740 #endif
741 }
742 XSETINT (val, value);
743 return val;
744 }
745
746
747 /* the rounding functions */
748
749 static Lisp_Object
750 rounding_driver (Lisp_Object arg, Lisp_Object divisor,
751 double (*double_round) (double),
752 EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
753 const char *name)
754 {
755 CHECK_NUMBER_OR_FLOAT (arg);
756
757 if (! NILP (divisor))
758 {
759 EMACS_INT i1, i2;
760
761 CHECK_NUMBER_OR_FLOAT (divisor);
762
763 if (FLOATP (arg) || FLOATP (divisor))
764 {
765 double f1, f2;
766
767 f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
768 f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
769 if (! IEEE_FLOATING_POINT && f2 == 0)
770 xsignal0 (Qarith_error);
771
772 IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
773 FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
774 return arg;
775 }
776
777 i1 = XINT (arg);
778 i2 = XINT (divisor);
779
780 if (i2 == 0)
781 xsignal0 (Qarith_error);
782
783 XSETINT (arg, (*int_round2) (i1, i2));
784 return arg;
785 }
786
787 if (FLOATP (arg))
788 {
789 double d;
790
791 IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
792 FLOAT_TO_INT (d, arg, name, arg);
793 }
794
795 return arg;
796 }
797
798 /* With C's /, the result is implementation-defined if either operand
799 is negative, so take care with negative operands in the following
800 integer functions. */
801
802 static EMACS_INT
803 ceiling2 (EMACS_INT i1, EMACS_INT i2)
804 {
805 return (i2 < 0
806 ? (i1 < 0 ? ((-1 - i1) / -i2) + 1 : - (i1 / -i2))
807 : (i1 <= 0 ? - (-i1 / i2) : ((i1 - 1) / i2) + 1));
808 }
809
810 static EMACS_INT
811 floor2 (EMACS_INT i1, EMACS_INT i2)
812 {
813 return (i2 < 0
814 ? (i1 <= 0 ? -i1 / -i2 : -1 - ((i1 - 1) / -i2))
815 : (i1 < 0 ? -1 - ((-1 - i1) / i2) : i1 / i2));
816 }
817
818 static EMACS_INT
819 truncate2 (EMACS_INT i1, EMACS_INT i2)
820 {
821 return (i2 < 0
822 ? (i1 < 0 ? -i1 / -i2 : - (i1 / -i2))
823 : (i1 < 0 ? - (-i1 / i2) : i1 / i2));
824 }
825
826 static EMACS_INT
827 round2 (EMACS_INT i1, EMACS_INT i2)
828 {
829 /* The C language's division operator gives us one remainder R, but
830 we want the remainder R1 on the other side of 0 if R1 is closer
831 to 0 than R is; because we want to round to even, we also want R1
832 if R and R1 are the same distance from 0 and if C's quotient is
833 odd. */
834 EMACS_INT q = i1 / i2;
835 EMACS_INT r = i1 % i2;
836 EMACS_INT abs_r = r < 0 ? -r : r;
837 EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
838 return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
839 }
840
841 /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
842 if `rint' exists but does not work right. */
843 #ifdef HAVE_RINT
844 #define emacs_rint rint
845 #else
846 static double
847 emacs_rint (double d)
848 {
849 return floor (d + 0.5);
850 }
851 #endif
852
853 static double
854 double_identity (double d)
855 {
856 return d;
857 }
858
859 DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
860 doc: /* Return the smallest integer no less than ARG.
861 This rounds the value towards +inf.
862 With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR. */)
863 (Lisp_Object arg, Lisp_Object divisor)
864 {
865 return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
866 }
867
868 DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
869 doc: /* Return the largest integer no greater than ARG.
870 This rounds the value towards -inf.
871 With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR. */)
872 (Lisp_Object arg, Lisp_Object divisor)
873 {
874 return rounding_driver (arg, divisor, floor, floor2, "floor");
875 }
876
877 DEFUN ("round", Fround, Sround, 1, 2, 0,
878 doc: /* Return the nearest integer to ARG.
879 With optional DIVISOR, return the nearest integer to ARG/DIVISOR.
880
881 Rounding a value equidistant between two integers may choose the
882 integer closer to zero, or it may prefer an even integer, depending on
883 your machine. For example, \(round 2.5\) can return 3 on some
884 systems, but 2 on others. */)
885 (Lisp_Object arg, Lisp_Object divisor)
886 {
887 return rounding_driver (arg, divisor, emacs_rint, round2, "round");
888 }
889
890 DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
891 doc: /* Truncate a floating point number to an int.
892 Rounds ARG toward zero.
893 With optional DIVISOR, truncate ARG/DIVISOR. */)
894 (Lisp_Object arg, Lisp_Object divisor)
895 {
896 return rounding_driver (arg, divisor, double_identity, truncate2,
897 "truncate");
898 }
899
900
901 Lisp_Object
902 fmod_float (Lisp_Object x, Lisp_Object y)
903 {
904 double f1, f2;
905
906 f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
907 f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
908
909 if (! IEEE_FLOATING_POINT && f2 == 0)
910 xsignal0 (Qarith_error);
911
912 /* If the "remainder" comes out with the wrong sign, fix it. */
913 IN_FLOAT2 ((f1 = fmod (f1, f2),
914 f1 = (f2 < 0 ? f1 > 0 : f1 < 0) ? f1 + f2 : f1),
915 "mod", x, y);
916 return make_float (f1);
917 }
918 \f
919 /* It's not clear these are worth adding. */
920
921 DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
922 doc: /* Return the smallest integer no less than ARG, as a float.
923 \(Round toward +inf.\) */)
924 (register Lisp_Object arg)
925 {
926 double d = extract_float (arg);
927 IN_FLOAT (d = ceil (d), "fceiling", arg);
928 return make_float (d);
929 }
930
931 DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
932 doc: /* Return the largest integer no greater than ARG, as a float.
933 \(Round towards -inf.\) */)
934 (register Lisp_Object arg)
935 {
936 double d = extract_float (arg);
937 IN_FLOAT (d = floor (d), "ffloor", arg);
938 return make_float (d);
939 }
940
941 DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
942 doc: /* Return the nearest integer to ARG, as a float. */)
943 (register Lisp_Object arg)
944 {
945 double d = extract_float (arg);
946 IN_FLOAT (d = emacs_rint (d), "fround", arg);
947 return make_float (d);
948 }
949
950 DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
951 doc: /* Truncate a floating point number to an integral float value.
952 Rounds the value toward zero. */)
953 (register Lisp_Object arg)
954 {
955 double d = extract_float (arg);
956 if (d >= 0.0)
957 IN_FLOAT (d = floor (d), "ftruncate", arg);
958 else
959 IN_FLOAT (d = ceil (d), "ftruncate", arg);
960 return make_float (d);
961 }
962 \f
963 #ifdef FLOAT_CATCH_SIGILL
964 static void
965 float_error (int signo)
966 {
967 if (! in_float)
968 fatal_error_signal (signo);
969
970 #ifdef BSD_SYSTEM
971 sigsetmask (SIGEMPTYMASK);
972 #else
973 /* Must reestablish handler each time it is called. */
974 signal (SIGILL, float_error);
975 #endif /* BSD_SYSTEM */
976
977 SIGNAL_THREAD_CHECK (signo);
978 in_float = 0;
979
980 xsignal1 (Qarith_error, float_error_arg);
981 }
982
983 /* Another idea was to replace the library function `infnan'
984 where SIGILL is signaled. */
985
986 #endif /* FLOAT_CATCH_SIGILL */
987
988 #ifdef HAVE_MATHERR
989 int
990 matherr (struct exception *x)
991 {
992 Lisp_Object args;
993 const char *name = x->name;
994
995 if (! in_float)
996 /* Not called from emacs-lisp float routines; do the default thing. */
997 return 0;
998 if (!strcmp (x->name, "pow"))
999 name = "expt";
1000
1001 args
1002 = Fcons (build_string (name),
1003 Fcons (make_float (x->arg1),
1004 ((!strcmp (name, "log") || !strcmp (name, "pow"))
1005 ? Fcons (make_float (x->arg2), Qnil)
1006 : Qnil)));
1007 switch (x->type)
1008 {
1009 case DOMAIN: xsignal (Qdomain_error, args); break;
1010 case SING: xsignal (Qsingularity_error, args); break;
1011 case OVERFLOW: xsignal (Qoverflow_error, args); break;
1012 case UNDERFLOW: xsignal (Qunderflow_error, args); break;
1013 default: xsignal (Qarith_error, args); break;
1014 }
1015 return (1); /* don't set errno or print a message */
1016 }
1017 #endif /* HAVE_MATHERR */
1018
1019 void
1020 init_floatfns (void)
1021 {
1022 #ifdef FLOAT_CATCH_SIGILL
1023 signal (SIGILL, float_error);
1024 #endif
1025 in_float = 0;
1026 }
1027
1028 void
1029 syms_of_floatfns (void)
1030 {
1031 defsubr (&Sacos);
1032 defsubr (&Sasin);
1033 defsubr (&Satan);
1034 defsubr (&Scos);
1035 defsubr (&Ssin);
1036 defsubr (&Stan);
1037 defsubr (&Sisnan);
1038 #ifdef HAVE_COPYSIGN
1039 defsubr (&Scopysign);
1040 defsubr (&Sfrexp);
1041 defsubr (&Sldexp);
1042 #endif
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);
1070
1071 defsubr (&Sabs);
1072 defsubr (&Sfloat);
1073 defsubr (&Slogb);
1074 defsubr (&Sceiling);
1075 defsubr (&Sfloor);
1076 defsubr (&Sround);
1077 defsubr (&Struncate);
1078 }