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