1 ;;; calc-math.el --- mathematical functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Colin Walters <walters@debian.org>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is distributed in the hope that it will be useful,
11 ;; but WITHOUT ANY WARRANTY. No author or distributor
12 ;; accepts responsibility to anyone for the consequences of using it
13 ;; or for whether it serves any particular purpose or works at all,
14 ;; unless he says so in writing. Refer to the GNU Emacs General Public
15 ;; License for full details.
17 ;; Everyone is granted permission to copy, modify and redistribute
18 ;; GNU Emacs, but only under the conditions described in the
19 ;; GNU Emacs General Public License. A copy of this license is
20 ;; supposed to have been given to you along with GNU Emacs so you
21 ;; can know your rights and responsibilities. It should be in a
22 ;; file named COPYING. Among other things, the copyright notice
23 ;; and this notice must be preserved on all copies.
29 ;; This file is autoloaded from calc-ext.el.
34 (defun calc-Need-calc-math () nil
)
37 (defun calc-sqrt (arg)
41 (calc-unary-op "^2" 'calcFunc-sqr arg
)
42 (calc-unary-op "sqrt" 'calcFunc-sqrt arg
))))
44 (defun calc-isqrt (arg)
48 (calc-unary-op "^2" 'calcFunc-sqr arg
)
49 (calc-unary-op "isqt" 'calcFunc-isqrt arg
))))
52 (defun calc-hypot (arg)
55 (calc-binary-op "hypt" 'calcFunc-hypot arg
)))
62 (defun calc-log10 (arg)
64 (calc-hyperbolic-func)
71 (calc-binary-op "alog" 'calcFunc-alog arg
)
72 (calc-binary-op "log" 'calcFunc-log arg
))))
74 (defun calc-ilog (arg)
78 (calc-binary-op "alog" 'calcFunc-alog arg
)
79 (calc-binary-op "ilog" 'calcFunc-ilog arg
))))
81 (defun calc-lnp1 (arg)
89 (if (calc-is-hyperbolic)
91 (calc-unary-op "lg10" 'calcFunc-log10 arg
)
92 (calc-unary-op "10^" 'calcFunc-exp10 arg
))
94 (calc-unary-op "ln" 'calcFunc-ln arg
)
95 (calc-unary-op "exp" 'calcFunc-exp arg
)))))
97 (defun calc-expm1 (arg)
100 (if (calc-is-inverse)
101 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg
)
102 (calc-unary-op "ex-1" 'calcFunc-expm1 arg
))))
107 (if (calc-is-inverse)
108 (if (calc-is-hyperbolic)
109 (if calc-symbolic-mode
110 (calc-pop-push-record 0 "phi" '(var phi var-phi
))
111 (calc-pop-push-record 0 "phi" (math-phi)))
112 (if calc-symbolic-mode
113 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma
))
114 (calc-pop-push-record 0 "gmma" (math-gamma-const))))
115 (if (calc-is-hyperbolic)
116 (if calc-symbolic-mode
117 (calc-pop-push-record 0 "e" '(var e var-e
))
118 (calc-pop-push-record 0 "e" (math-e)))
119 (if calc-symbolic-mode
120 (calc-pop-push-record 0 "pi" '(var pi var-pi
))
121 (calc-pop-push-record 0 "pi" (math-pi)))))))
123 (defun calc-sin (arg)
126 (if (calc-is-hyperbolic)
127 (if (calc-is-inverse)
128 (calc-unary-op "asnh" 'calcFunc-arcsinh arg
)
129 (calc-unary-op "sinh" 'calcFunc-sinh arg
))
130 (if (calc-is-inverse)
131 (calc-unary-op "asin" 'calcFunc-arcsin arg
)
132 (calc-unary-op "sin" 'calcFunc-sin arg
)))))
134 (defun calc-arcsin (arg)
139 (defun calc-sinh (arg)
141 (calc-hyperbolic-func)
144 (defun calc-arcsinh (arg)
147 (calc-hyperbolic-func)
150 (defun calc-cos (arg)
153 (if (calc-is-hyperbolic)
154 (if (calc-is-inverse)
155 (calc-unary-op "acsh" 'calcFunc-arccosh arg
)
156 (calc-unary-op "cosh" 'calcFunc-cosh arg
))
157 (if (calc-is-inverse)
158 (calc-unary-op "acos" 'calcFunc-arccos arg
)
159 (calc-unary-op "cos" 'calcFunc-cos arg
)))))
161 (defun calc-arccos (arg)
166 (defun calc-cosh (arg)
168 (calc-hyperbolic-func)
171 (defun calc-arccosh (arg)
174 (calc-hyperbolic-func)
177 (defun calc-sincos ()
180 (if (calc-is-inverse)
181 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos
(calc-top-n 1)))
182 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos
(calc-top-n 1))))))
184 (defun calc-tan (arg)
187 (if (calc-is-hyperbolic)
188 (if (calc-is-inverse)
189 (calc-unary-op "atnh" 'calcFunc-arctanh arg
)
190 (calc-unary-op "tanh" 'calcFunc-tanh arg
))
191 (if (calc-is-inverse)
192 (calc-unary-op "atan" 'calcFunc-arctan arg
)
193 (calc-unary-op "tan" 'calcFunc-tan arg
)))))
195 (defun calc-arctan (arg)
200 (defun calc-tanh (arg)
202 (calc-hyperbolic-func)
205 (defun calc-arctanh (arg)
208 (calc-hyperbolic-func)
211 (defun calc-arctan2 ()
214 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2
(calc-top-list-n 2)))))
216 (defun calc-conj (arg)
219 (calc-unary-op "conj" 'calcFunc-conj arg
)))
221 (defun calc-imaginary ()
224 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
228 (defun calc-to-degrees (arg)
231 (calc-unary-op ">deg" 'calcFunc-deg arg
)))
233 (defun calc-to-radians (arg)
236 (calc-unary-op ">rad" 'calcFunc-rad arg
)))
239 (defun calc-degrees-mode (arg)
243 (calc-change-mode 'calc-angle-mode
'deg
)
244 (message "Angles measured in degrees")))
245 ((= arg
2) (calc-radians-mode))
246 ((= arg
3) (calc-hms-mode))
247 (t (error "Prefix argument out of range"))))
249 (defun calc-radians-mode ()
252 (calc-change-mode 'calc-angle-mode
'rad
)
253 (message "Angles measured in radians")))
256 ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public]
257 ;;; This method takes advantage of the fact that Newton's method starting
258 ;;; with an overestimate always works, even using truncating integer division!
259 (defun math-isqrt (a)
260 (cond ((Math-zerop a
) a
)
261 ((not (math-natnump a
))
262 (math-reject-arg a
'natnump
))
264 (math-isqrt-small a
))
266 (math-normalize (cons 'bigpos
(cdr (math-isqrt-bignum (cdr a
))))))))
268 (defun calcFunc-isqrt (a)
270 (math-isqrt (math-floor a
))
271 (math-floor (math-sqrt a
))))
274 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
275 (defun math-isqrt-bignum (a) ; [P.l L]
276 (let ((len (length a
)))
278 (let* ((top (nthcdr (- len
2) a
)))
279 (math-isqrt-bignum-iter
283 (1+ (math-isqrt-small
284 (+ (* (nth 1 top
) 1000) (car top
)))))
286 (let* ((top (nth (1- len
) a
)))
287 (math-isqrt-bignum-iter
290 (list (1+ (math-isqrt-small top
)))
293 (defun math-isqrt-bignum-iter (a guess
) ; [l L l]
294 (math-working "isqrt" (cons 'bigpos guess
))
295 (let* ((q (math-div-bignum a guess
))
296 (s (math-add-bignum (car q
) guess
))
297 (g2 (math-div2-bignum s
))
298 (comp (math-compare-bignum g2 guess
)))
300 (math-isqrt-bignum-iter a g2
)
301 (cons (and (= comp
0)
302 (math-zerop-bignum (cdr q
))
306 (defun math-zerop-bignum (a)
309 (while (eq (car (setq a
(cdr a
))) 0))
312 (defun math-scale-bignum-3 (a n
) ; [L L S]
318 (defun math-isqrt-small (a) ; A > 0. [S S]
319 (let ((g (cond ((>= a
10000) 1000)
323 (while (< (setq g2
(/ (+ g
(/ a g
)) 2)) g
)
330 ;;; Compute the square root of a number.
331 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
334 (and (Math-zerop a
) a
)
335 (and (math-known-nonposp a
)
336 (math-imaginary (math-sqrt (math-neg a
))))
338 (let ((sqrt (math-isqrt-small a
)))
339 (if (= (* sqrt sqrt
) a
)
341 (if calc-symbolic-mode
342 (list 'calcFunc-sqrt a
)
343 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
344 (and (eq (car-safe a
) 'bigpos
)
345 (let* ((res (math-isqrt-bignum (cdr a
)))
346 (sqrt (math-normalize (cons 'bigpos
(cdr res
)))))
349 (if calc-symbolic-mode
350 (list 'calcFunc-sqrt a
)
351 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
352 (and (eq (car-safe a
) 'frac
)
353 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a
)))))
354 (num-sqrt (math-normalize (cons 'bigpos
(cdr num-res
))))
355 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a
)))))
356 (den-sqrt (math-normalize (cons 'bigpos
(cdr den-res
)))))
357 (if (and (car num-res
) (car den-res
))
358 (list 'frac num-sqrt den-sqrt
)
359 (if calc-symbolic-mode
360 (if (or (car num-res
) (car den-res
))
361 (math-div (if (car num-res
)
362 num-sqrt
(list 'calcFunc-sqrt
(nth 1 a
)))
364 den-sqrt
(list 'calcFunc-sqrt
(nth 2 a
))))
365 (list 'calcFunc-sqrt a
))
366 (math-sqrt-float (math-float a
)
367 (math-div (math-float num-sqrt
) den-sqrt
))))))
368 (and (eq (car-safe a
) 'float
)
369 (if calc-symbolic-mode
370 (if (= (%
(nth 2 a
) 2) 0)
371 (let ((res (math-isqrt-bignum
372 (cdr (Math-bignum-test (nth 1 a
))))))
374 (math-make-float (math-normalize
375 (cons 'bigpos
(cdr res
)))
377 (signal 'inexact-result nil
)))
378 (signal 'inexact-result nil
))
379 (math-sqrt-float a
)))
380 (and (eq (car-safe a
) 'cplx
)
381 (math-with-extra-prec 2
382 (let* ((d (math-abs a
))
383 (imag (math-sqrt (math-mul (math-sub d
(nth 1 a
))
386 (math-sqrt (math-mul (math-add d
(nth 1 a
)) '(float 5 -
1)))
387 (if (math-negp (nth 2 a
)) (math-neg imag
) imag
)))))
388 (and (eq (car-safe a
) 'polar
)
390 (math-sqrt (nth 1 a
))
391 (math-mul (nth 2 a
) '(float 5 -
1))))
392 (and (eq (car-safe a
) 'sdev
)
393 (let ((sqrt (math-sqrt (nth 1 a
))))
395 (math-div (nth 2 a
) (math-mul sqrt
2)))))
396 (and (eq (car-safe a
) 'intv
)
397 (not (math-negp (nth 2 a
)))
398 (math-make-intv (nth 1 a
) (math-sqrt (nth 2 a
)) (math-sqrt (nth 3 a
))))
399 (and (eq (car-safe a
) '*)
400 (or (math-known-nonnegp (nth 1 a
))
401 (math-known-nonnegp (nth 2 a
)))
402 (math-mul (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
403 (and (eq (car-safe a
) '/)
404 (or (and (math-known-nonnegp (nth 2 a
))
405 (math-div (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
406 (and (math-known-nonnegp (nth 1 a
))
407 (not (math-equal-int (nth 1 a
) 1))
408 (math-mul (math-sqrt (nth 1 a
))
409 (math-sqrt (math-div 1 (nth 2 a
)))))))
410 (and (eq (car-safe a
) '^
)
411 (math-known-evenp (nth 2 a
))
412 (math-known-realp (nth 1 a
))
413 (math-abs (math-pow (nth 1 a
) (math-div (nth 2 a
) 2))))
414 (let ((inf (math-infinitep a
)))
416 (math-mul (math-sqrt (math-infinite-dir a inf
)) inf
)))
418 (calc-record-why 'numberp a
)
419 (list 'calcFunc-sqrt a
))))
420 (defalias 'calcFunc-sqrt
'math-sqrt
)
422 (defun math-infinite-dir (a &optional inf
)
423 (or inf
(setq inf
(math-infinitep a
)))
424 (math-normalize (math-expr-subst a inf
1)))
426 (defun math-sqrt-float (a &optional guess
) ; [F F F]
427 (if calc-symbolic-mode
428 (signal 'inexact-result nil
)
429 (math-with-extra-prec 1 (math-sqrt-raw a guess
))))
431 (defun math-sqrt-raw (a &optional guess
) ; [F F F]
432 (if (not (Math-posp a
))
435 (let ((ldiff (- (math-numdigs (nth 1 a
)) 6)))
436 (or (= (%
(+ (nth 2 a
) ldiff
) 2) 0) (setq ldiff
(1+ ldiff
)))
437 (setq guess
(math-make-float (math-isqrt-small
438 (math-scale-int (nth 1 a
) (- ldiff
)))
439 (/ (+ (nth 2 a
) ldiff
) 2)))))
440 (math-sqrt-float-iter a guess
)))
442 (defun math-sqrt-float-iter (a guess
) ; [F F F]
443 (math-working "sqrt" guess
)
444 (let ((g2 (math-mul-float (math-add-float guess
(math-div-float a guess
))
446 (if (math-nearly-equal-float g2 guess
)
448 (math-sqrt-float-iter a g2
))))
450 ;;; True if A and B differ only in the last digit of precision. [P F F]
451 (defun math-nearly-equal-float (a b
)
452 (let ((ediff (- (nth 2 a
) (nth 2 b
))))
453 (cond ((= ediff
0) ;; Expanded out for speed
454 (setq ediff
(math-add (Math-integer-neg (nth 1 a
)) (nth 1 b
)))
456 (and (not (consp ediff
))
459 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))
461 (setq ediff
(math-add (Math-integer-neg (nth 1 b
))
462 (math-scale-int (nth 1 a
) 1)))
463 (and (not (consp ediff
))
466 (= (math-numdigs (nth 1 b
)) calc-internal-prec
)))
468 (setq ediff
(math-add (Math-integer-neg (nth 1 a
))
469 (math-scale-int (nth 1 b
) 1)))
470 (and (not (consp ediff
))
473 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))))
475 (defun math-nearly-equal (a b
) ; [P N N] [Public]
476 (setq a
(math-float a
))
477 (setq b
(math-float b
))
478 (if (eq (car a
) 'polar
) (setq a
(math-complex a
)))
479 (if (eq (car b
) 'polar
) (setq b
(math-complex b
)))
480 (if (eq (car a
) 'cplx
)
481 (if (eq (car b
) 'cplx
)
482 (and (or (math-nearly-equal-float (nth 1 a
) (nth 1 b
))
483 (and (math-nearly-zerop-float (nth 1 a
) (nth 2 a
))
484 (math-nearly-zerop-float (nth 1 b
) (nth 2 b
))))
485 (or (math-nearly-equal-float (nth 2 a
) (nth 2 b
))
486 (and (math-nearly-zerop-float (nth 2 a
) (nth 1 a
))
487 (math-nearly-zerop-float (nth 2 b
) (nth 1 b
)))))
488 (and (math-nearly-equal-float (nth 1 a
) b
)
489 (math-nearly-zerop-float (nth 2 a
) b
)))
490 (if (eq (car b
) 'cplx
)
491 (and (math-nearly-equal-float a
(nth 1 b
))
492 (math-nearly-zerop-float a
(nth 2 b
)))
493 (math-nearly-equal-float a b
))))
495 ;;; True if A is nearly zero compared to B. [P F F]
496 (defun math-nearly-zerop-float (a b
)
498 (<= (+ (math-numdigs (nth 1 a
)) (nth 2 a
))
499 (1+ (- (+ (math-numdigs (nth 1 b
)) (nth 2 b
)) calc-internal-prec
)))))
501 (defun math-nearly-zerop (a b
) ; [P N R] [Public]
502 (setq a
(math-float a
))
503 (setq b
(math-float b
))
504 (if (eq (car a
) 'cplx
)
505 (and (math-nearly-zerop-float (nth 1 a
) b
)
506 (math-nearly-zerop-float (nth 2 a
) b
))
507 (if (eq (car a
) 'polar
)
508 (math-nearly-zerop-float (nth 1 a
) b
)
509 (math-nearly-zerop-float a b
))))
511 ;;; This implementation could be improved, accuracy-wise.
512 (defun math-hypot (a b
)
513 (cond ((Math-zerop a
) (math-abs b
))
514 ((Math-zerop b
) (math-abs a
))
515 ((not (Math-scalarp a
))
516 (if (math-infinitep a
)
517 (if (math-infinitep b
)
522 (calc-record-why 'scalarp a
)
523 (list 'calcFunc-hypot a b
)))
524 ((not (Math-scalarp b
))
525 (if (math-infinitep b
)
527 (calc-record-why 'scalarp b
)
528 (list 'calcFunc-hypot a b
)))
529 ((and (Math-numberp a
) (Math-numberp b
))
530 (math-with-extra-prec 1
531 (math-sqrt (math-add (calcFunc-abssqr a
) (calcFunc-abssqr b
)))))
532 ((eq (car-safe a
) 'hms
)
533 (if (eq (car-safe b
) 'hms
) ; this helps sdev's of hms forms
534 (math-to-hms (math-hypot (math-from-hms a
'deg
)
535 (math-from-hms b
'deg
)))
536 (math-to-hms (math-hypot (math-from-hms a
'deg
) b
))))
537 ((eq (car-safe b
) 'hms
)
538 (math-to-hms (math-hypot a
(math-from-hms b
'deg
))))
540 (defalias 'calcFunc-hypot
'math-hypot
)
542 (defun calcFunc-sqr (x)
547 (defun math-nth-root (a n
)
548 (cond ((= n
2) (math-sqrt a
))
552 (let ((root (math-nth-root-integer a n
)))
555 (and (not calc-symbolic-mode
)
556 (math-nth-root-float (math-float a
) n
557 (math-float (cdr root
)))))))
558 ((eq (car-safe a
) 'frac
)
559 (let* ((num-root (math-nth-root-integer (nth 1 a
) n
))
560 (den-root (math-nth-root-integer (nth 2 a
) n
)))
561 (if (and (car num-root
) (car den-root
))
562 (list 'frac
(cdr num-root
) (cdr den-root
))
563 (and (not calc-symbolic-mode
)
566 (math-div-float (math-float (cdr num-root
))
567 (math-float (cdr den-root
))))))))
568 ((eq (car-safe a
) 'float
)
569 (and (not calc-symbolic-mode
)
570 (math-nth-root-float a n
)))
571 ((eq (car-safe a
) 'polar
)
572 (let ((root (math-nth-root (nth 1 a
) n
)))
573 (and root
(list 'polar root
(math-div (nth 2 a
) n
)))))
576 (defun math-nth-root-float (a n
&optional guess
)
577 (math-inexact-result)
578 (math-with-extra-prec 1
579 (let ((nf (math-float n
))
580 (nfm1 (math-float (1- n
))))
581 (math-nth-root-float-iter a
(or guess
583 1 (/ (+ (math-numdigs (nth 1 a
))
588 (defun math-nth-root-float-iter (a guess
) ; uses "n", "nf", "nfm1"
589 (math-working "root" guess
)
590 (let ((g2 (math-div-float (math-add-float (math-mul nfm1 guess
)
592 a
(math-ipow guess
(1- n
))))
594 (if (math-nearly-equal-float g2 guess
)
596 (math-nth-root-float-iter a g2
))))
598 (defun math-nth-root-integer (a n
&optional guess
) ; [I I S]
599 (math-nth-root-int-iter a
(or guess
600 (math-scale-int 1 (/ (+ (math-numdigs a
)
604 (defun math-nth-root-int-iter (a guess
) ; uses "n"
605 (math-working "root" guess
)
606 (let* ((q (math-idivmod a
(math-ipow guess
(1- n
))))
607 (s (math-add (car q
) (math-mul (1- n
) guess
)))
608 (g2 (math-idivmod s n
)))
609 (if (Math-natnum-lessp (car g2
) guess
)
610 (math-nth-root-int-iter a
(car g2
))
611 (cons (and (equal (car g2
) guess
)
616 (defun calcFunc-nroot (x n
)
617 (calcFunc-pow x
(if (integerp n
)
624 ;;;; Transcendental functions.
626 ;;; All of these functions are defined on the complex plane.
627 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
629 ;;; Most functions increase calc-internal-prec by 2 digits, then round
630 ;;; down afterward. "-raw" functions use the current precision, require
631 ;;; their arguments to be in float (or complex float) format, and always
632 ;;; work in radians (where applicable).
634 (defun math-to-radians (a) ; [N N]
635 (cond ((eq (car-safe a
) 'hms
)
636 (math-from-hms a
'rad
))
637 ((memq calc-angle-mode
'(deg hms
))
638 (math-mul a
(math-pi-over-180)))
641 (defun math-from-radians (a) ; [N N]
642 (cond ((eq calc-angle-mode
'deg
)
644 (math-div a
(math-pi-over-180))
645 (list 'calcFunc-deg a
)))
646 ((eq calc-angle-mode
'hms
)
647 (math-to-hms a
'rad
))
650 (defun math-to-radians-2 (a) ; [N N]
651 (cond ((eq (car-safe a
) 'hms
)
652 (math-from-hms a
'rad
))
653 ((memq calc-angle-mode
'(deg hms
))
654 (if calc-symbolic-mode
655 (math-div (math-mul a
'(var pi var-pi
)) 180)
656 (math-mul a
(math-pi-over-180))))
659 (defun math-from-radians-2 (a) ; [N N]
660 (cond ((memq calc-angle-mode
'(deg hms
))
661 (if calc-symbolic-mode
662 (math-div (math-mul 180 a
) '(var pi var-pi
))
663 (math-div a
(math-pi-over-180))))
668 ;;; Sine, cosine, and tangent.
670 (defun calcFunc-sin (x) ; [N N] [Public]
671 (cond ((and (integerp x
)
672 (if (eq calc-angle-mode
'deg
)
675 (aref [0 1 0 -
1] (math-mod (/ x
90) 4)))
677 (math-with-extra-prec 2
678 (math-sin-raw (math-to-radians (math-float x
)))))
681 (math-with-extra-prec 2
682 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
683 (xs (math-to-radians (math-float (nth 2 x
))))
684 (sc (math-sin-cos-raw xx
)))
685 (math-make-sdev (car sc
) (math-mul xs
(cdr sc
)))))
686 (math-make-sdev (calcFunc-sin (nth 1 x
))
687 (math-mul (nth 2 x
) (calcFunc-cos (nth 1 x
))))))
688 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
689 (calcFunc-cos (math-sub x
(math-quarter-circle nil
))))
690 ((equal x
'(var nan var-nan
))
692 (t (calc-record-why 'scalarp x
)
693 (list 'calcFunc-sin x
))))
695 (defun calcFunc-cos (x) ; [N N] [Public]
696 (cond ((and (integerp x
)
697 (if (eq calc-angle-mode
'deg
)
700 (aref [1 0 -
1 0] (math-mod (/ x
90) 4)))
702 (math-with-extra-prec 2
703 (math-cos-raw (math-to-radians (math-float x
)))))
706 (math-with-extra-prec 2
707 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
708 (xs (math-to-radians (math-float (nth 2 x
))))
709 (sc (math-sin-cos-raw xx
)))
710 (math-make-sdev (cdr sc
) (math-mul xs
(car sc
)))))
711 (math-make-sdev (calcFunc-cos (nth 1 x
))
712 (math-mul (nth 2 x
) (calcFunc-sin (nth 1 x
))))))
713 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
714 (math-with-extra-prec 2
715 (let* ((xx (math-to-radians (math-float x
)))
716 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
717 (nb (math-floor (math-div (nth 3 xx
) (math-pi))))
718 (span (math-sub nb na
)))
719 (if (memq span
'(0 1))
720 (let ((int (math-sort-intv (nth 1 x
)
721 (math-cos-raw (nth 2 xx
))
722 (math-cos-raw (nth 3 xx
)))))
725 (math-make-intv (logior (nth 1 x
) 2)
728 (math-make-intv (logior (nth 1 x
) 1)
732 (list 'intv
3 -
1 1)))))
733 ((equal x
'(var nan var-nan
))
735 (t (calc-record-why 'scalarp x
)
736 (list 'calcFunc-cos x
))))
738 (defun calcFunc-sincos (x) ; [V N] [Public]
740 (math-with-extra-prec 2
741 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x
)))))
742 (list 'vec
(cdr sc
) (car sc
)))) ; the vector [cos, sin]
743 (list 'vec
(calcFunc-sin x
) (calcFunc-cos x
))))
745 (defun calcFunc-tan (x) ; [N N] [Public]
746 (cond ((and (integerp x
)
747 (if (eq calc-angle-mode
'deg
)
752 (math-with-extra-prec 2
753 (math-tan-raw (math-to-radians (math-float x
)))))
756 (math-with-extra-prec 2
757 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
758 (xs (math-to-radians (math-float (nth 2 x
))))
759 (sc (math-sin-cos-raw xx
)))
760 (if (and (math-zerop (cdr sc
)) (not calc-infinite-mode
))
762 (calc-record-why "*Division by zero")
763 (list 'calcFunc-tan x
))
764 (math-make-sdev (math-div-float (car sc
) (cdr sc
))
765 (math-div-float xs
(math-sqr (cdr sc
)))))))
766 (math-make-sdev (calcFunc-tan (nth 1 x
))
768 (math-sqr (calcFunc-cos (nth 1 x
)))))))
769 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
770 (or (math-with-extra-prec 2
771 (let* ((xx (math-to-radians (math-float x
)))
772 (na (math-floor (math-div (math-sub (nth 2 xx
)
775 (nb (math-floor (math-div (math-sub (nth 3 xx
)
779 (math-sort-intv (nth 1 x
)
780 (math-tan-raw (nth 2 xx
))
781 (math-tan-raw (nth 3 xx
))))))
782 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))))
783 ((equal x
'(var nan var-nan
))
785 (t (calc-record-why 'scalarp x
)
786 (list 'calcFunc-tan x
))))
788 (defun math-sin-raw (x) ; [N N]
789 (cond ((eq (car x
) 'cplx
)
790 (let* ((expx (math-exp-raw (nth 2 x
)))
791 (expmx (math-div-float '(float 1 0) expx
))
792 (sc (math-sin-cos-raw (nth 1 x
))))
794 (math-mul-float (car sc
)
795 (math-mul-float (math-add-float expx expmx
)
797 (math-mul-float (cdr sc
)
798 (math-mul-float (math-sub-float expx expmx
)
801 (math-polar (math-sin-raw (math-complex x
))))
802 ((Math-integer-negp (nth 1 x
))
803 (math-neg-float (math-sin-raw (math-neg-float x
))))
804 ((math-lessp-float '(float 7 0) x
) ; avoid inf loops due to roundoff
805 (math-sin-raw (math-mod x
(math-two-pi))))
806 (t (math-sin-raw-2 x x
))))
808 (defun math-cos-raw (x) ; [N N]
809 (if (eq (car-safe x
) 'polar
)
810 (math-polar (math-cos-raw (math-complex x
)))
811 (math-sin-raw (math-sub (math-pi-over-2) x
))))
813 ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
814 ;;; compute either sin(x) or cos(x), whichever is smaller, and compute
815 ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
816 (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
817 (cons (math-sin-raw x
) (math-cos-raw x
)))
819 (defun math-tan-raw (x) ; [N N]
820 (cond ((eq (car x
) 'cplx
)
821 (let* ((x (math-mul x
'(float 2 0)))
822 (expx (math-exp-raw (nth 2 x
)))
823 (expmx (math-div-float '(float 1 0) expx
))
824 (sc (math-sin-cos-raw (nth 1 x
)))
825 (d (math-add-float (cdr sc
)
826 (math-mul-float (math-add-float expx expmx
)
828 (and (not (eq (nth 1 d
) 0))
830 (math-div-float (car sc
) d
)
831 (math-div-float (math-mul-float (math-sub-float expx
833 '(float 5 -
1)) d
)))))
835 (math-polar (math-tan-raw (math-complex x
))))
837 (let ((sc (math-sin-cos-raw x
)))
838 (if (eq (nth 1 (cdr sc
)) 0)
839 (math-div (car sc
) 0)
840 (math-div-float (car sc
) (cdr sc
)))))))
842 (defun math-sin-raw-2 (x orgx
) ; This avoids poss of inf recursion. [F F]
843 (let ((xmpo2 (math-sub-float (math-pi-over-2) x
)))
844 (cond ((Math-integer-negp (nth 1 xmpo2
))
845 (math-neg-float (math-sin-raw-2 (math-sub-float x
(math-pi))
847 ((math-lessp-float (math-pi-over-4) x
)
848 (math-cos-raw-2 xmpo2 orgx
))
849 ((math-lessp-float x
(math-neg (math-pi-over-4)))
850 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x
) orgx
)))
851 ((math-nearly-zerop-float x orgx
) '(float 0 0))
852 (calc-symbolic-mode (signal 'inexact-result nil
))
853 (t (math-sin-series x
6 4 x
(math-neg-float (math-sqr-float x
)))))))
855 (defun math-cos-raw-2 (x orgx
) ; [F F]
856 (cond ((math-nearly-zerop-float x orgx
) '(float 1 0))
857 (calc-symbolic-mode (signal 'inexact-result nil
))
858 (t (let ((xnegsqr (math-neg-float (math-sqr-float x
))))
860 (math-add-float '(float 1 0)
861 (math-mul-float xnegsqr
'(float 5 -
1)))
862 24 5 xnegsqr xnegsqr
)))))
864 (defun math-sin-series (sum nfac n x xnegsqr
)
865 (math-working "sin" sum
)
866 (let* ((nextx (math-mul-float x xnegsqr
))
867 (nextsum (math-add-float sum
(math-div-float nextx
868 (math-float nfac
)))))
869 (if (math-nearly-equal-float sum nextsum
)
871 (math-sin-series nextsum
(math-mul nfac
(* n
(1+ n
)))
872 (+ n
2) nextx xnegsqr
))))
875 ;;; Inverse sine, cosine, tangent.
877 (defun calcFunc-arcsin (x) ; [N N] [Public]
879 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 90)
880 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
90)
881 (calc-symbolic-mode (signal 'inexact-result nil
))
883 (math-with-extra-prec 2
884 (math-from-radians (math-arcsin-raw (math-float x
)))))
886 (math-make-sdev (calcFunc-arcsin (nth 1 x
))
890 (math-sub 1 (math-sqr (nth 1 x
))))))))
892 (math-sort-intv (nth 1 x
)
893 (calcFunc-arcsin (nth 2 x
))
894 (calcFunc-arcsin (nth 3 x
))))
895 ((equal x
'(var nan var-nan
))
897 (t (calc-record-why 'numberp x
)
898 (list 'calcFunc-arcsin x
))))
900 (defun calcFunc-arccos (x) ; [N N] [Public]
902 ((and (eq x
0) (eq calc-angle-mode
'deg
)) 90)
903 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) 180)
904 (calc-symbolic-mode (signal 'inexact-result nil
))
906 (math-with-extra-prec 2
907 (math-from-radians (math-arccos-raw (math-float x
)))))
909 (math-make-sdev (calcFunc-arccos (nth 1 x
))
913 (math-sub 1 (math-sqr (nth 1 x
))))))))
915 (math-sort-intv (nth 1 x
)
916 (calcFunc-arccos (nth 2 x
))
917 (calcFunc-arccos (nth 3 x
))))
918 ((equal x
'(var nan var-nan
))
920 (t (calc-record-why 'numberp x
)
921 (list 'calcFunc-arccos x
))))
923 (defun calcFunc-arctan (x) ; [N N] [Public]
925 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 45)
926 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
45)
928 (math-with-extra-prec 2
929 (math-from-radians (math-arctan-raw (math-float x
)))))
931 (math-make-sdev (calcFunc-arctan (nth 1 x
))
934 (math-add 1 (math-sqr (nth 1 x
)))))))
936 (math-sort-intv (nth 1 x
)
937 (calcFunc-arctan (nth 2 x
))
938 (calcFunc-arctan (nth 3 x
))))
939 ((equal x
'(var inf var-inf
))
940 (math-quarter-circle t
))
941 ((equal x
'(neg (var inf var-inf
)))
942 (math-neg (math-quarter-circle t
)))
943 ((equal x
'(var nan var-nan
))
945 (t (calc-record-why 'numberp x
)
946 (list 'calcFunc-arctan x
))))
948 (defun math-arcsin-raw (x) ; [N N]
949 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x
)))))
950 (if (or (memq (car x
) '(cplx polar
))
951 (memq (car a
) '(cplx polar
)))
952 (math-with-extra-prec 2 ; use extra precision for difficult case
953 (math-mul '(cplx 0 -
1)
954 (math-ln-raw (math-add (math-mul '(cplx 0 1) x
) a
))))
955 (math-arctan2-raw x a
))))
957 (defun math-arccos-raw (x) ; [N N]
958 (math-sub (math-pi-over-2) (math-arcsin-raw x
)))
960 (defun math-arctan-raw (x) ; [N N]
961 (cond ((memq (car x
) '(cplx polar
))
962 (math-with-extra-prec 2 ; extra-extra
964 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x
)))
965 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -
1) x
))))
967 ((Math-integer-negp (nth 1 x
))
968 (math-neg-float (math-arctan-raw (math-neg-float x
))))
970 (calc-symbolic-mode (signal 'inexact-result nil
))
971 ((math-equal-int x
1) (math-pi-over-4))
972 ((math-equal-int x -
1) (math-neg (math-pi-over-4)))
973 ((math-lessp-float '(float 414214 -
6) x
) ; if x > sqrt(2) - 1, reduce
974 (if (math-lessp-float '(float 1 0) x
)
975 (math-sub-float (math-mul-float (math-pi) '(float 5 -
1))
976 (math-arctan-raw (math-div-float '(float 1 0) x
)))
977 (math-sub-float (math-mul-float (math-pi) '(float 25 -
2))
978 (math-arctan-raw (math-div-float
979 (math-sub-float '(float 1 0) x
)
980 (math-add-float '(float 1 0)
982 (t (math-arctan-series x
3 x
(math-neg-float (math-sqr-float x
))))))
984 (defun math-arctan-series (sum n x xnegsqr
)
985 (math-working "arctan" sum
)
986 (let* ((nextx (math-mul-float x xnegsqr
))
987 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
988 (if (math-nearly-equal-float sum nextsum
)
990 (math-arctan-series nextsum
(+ n
2) nextx xnegsqr
))))
992 (defun calcFunc-arctan2 (y x
) ; [F R R] [Public]
995 (math-with-extra-prec 2
996 (math-from-radians (math-arctan2-raw (math-float y
)
998 (calc-record-why 'anglep x
)
999 (list 'calcFunc-arctan2 y x
))
1000 (if (and (or (math-infinitep x
) (math-anglep x
))
1001 (or (math-infinitep y
) (math-anglep y
)))
1016 (calcFunc-arctan2 y x
)
1017 '(var nan var-nan
)))
1018 (calc-record-why 'anglep y
)
1019 (list 'calcFunc-arctan2 y x
))))
1021 (defun math-arctan2-raw (y x
) ; [F R R]
1022 (cond ((math-zerop y
)
1023 (if (math-negp x
) (math-pi)
1024 (if (or (math-floatp x
) (math-floatp y
)) '(float 0 0) 0)))
1028 (math-neg (math-pi-over-2))))
1030 (math-arctan-raw (math-div-float y x
)))
1032 (math-add-float (math-arctan-raw (math-div-float y x
))
1035 (math-sub-float (math-arctan-raw (math-div-float y x
))
1038 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1039 (if (and (Math-vectorp x
)
1041 (calcFunc-arctan2 (nth 2 x
) (nth 1 x
))
1042 (math-reject-arg x
"*Two-element vector expected")))
1046 ;;; Exponential function.
1048 (defun calcFunc-exp (x) ; [N N] [Public]
1050 ((and (memq x
'(1 -
1)) calc-symbolic-mode
)
1051 (if (eq x
1) '(var e var-e
) (math-div 1 '(var e var-e
))))
1053 (math-with-extra-prec 2 (math-exp-raw (math-float x
))))
1054 ((eq (car-safe x
) 'sdev
)
1055 (let ((ex (calcFunc-exp (nth 1 x
))))
1056 (math-make-sdev ex
(math-mul (nth 2 x
) ex
))))
1057 ((eq (car-safe x
) 'intv
)
1058 (math-make-intv (nth 1 x
) (calcFunc-exp (nth 2 x
))
1059 (calcFunc-exp (nth 3 x
))))
1060 ((equal x
'(var inf var-inf
))
1062 ((equal x
'(neg (var inf var-inf
)))
1064 ((equal x
'(var nan var-nan
))
1066 (t (calc-record-why 'numberp x
)
1067 (list 'calcFunc-exp x
))))
1069 (defun calcFunc-expm1 (x) ; [N N] [Public]
1071 ((math-zerop x
) '(float 0 0))
1072 (calc-symbolic-mode (signal 'inexact-result nil
))
1074 (math-with-extra-prec 2
1075 (let ((x (math-float x
)))
1076 (if (and (eq (car x
) 'float
)
1077 (math-lessp-float x
'(float 1 0))
1078 (math-lessp-float '(float -
1 0) x
))
1079 (math-exp-minus-1-raw x
)
1080 (math-add (math-exp-raw x
) -
1)))))
1081 ((eq (car-safe x
) 'sdev
)
1083 (let ((ex (calcFunc-expm1 (nth 1 x
))))
1084 (math-make-sdev ex
(math-mul (nth 2 x
) (math-add ex
1))))
1085 (math-make-sdev (calcFunc-expm1 (nth 1 x
))
1086 (math-mul (nth 2 x
) (calcFunc-exp (nth 1 x
))))))
1087 ((eq (car-safe x
) 'intv
)
1088 (math-make-intv (nth 1 x
)
1089 (calcFunc-expm1 (nth 2 x
))
1090 (calcFunc-expm1 (nth 3 x
))))
1091 ((equal x
'(var inf var-inf
))
1093 ((equal x
'(neg (var inf var-inf
)))
1095 ((equal x
'(var nan var-nan
))
1097 (t (calc-record-why 'numberp x
)
1098 (list 'calcFunc-expm1 x
))))
1100 (defun calcFunc-exp10 (x) ; [N N] [Public]
1103 (math-pow '(float 1 1) x
)))
1105 (defun math-exp-raw (x) ; [N N]
1106 (cond ((math-zerop x
) '(float 1 0))
1107 (calc-symbolic-mode (signal 'inexact-result nil
))
1109 (let ((expx (math-exp-raw (nth 1 x
)))
1110 (sc (math-sin-cos-raw (nth 2 x
))))
1112 (math-mul-float expx
(cdr sc
))
1113 (math-mul-float expx
(car sc
)))))
1114 ((eq (car x
) 'polar
)
1115 (let ((xc (math-complex x
)))
1117 (math-exp-raw (nth 1 xc
))
1118 (math-from-radians (nth 2 xc
)))))
1119 ((or (math-lessp-float '(float 5 -
1) x
)
1120 (math-lessp-float x
'(float -
5 -
1)))
1121 (if (math-lessp-float '(float 921035 1) x
)
1123 (if (math-lessp-float x
'(float -
921035 1))
1125 (let* ((two-x (math-mul-float x
'(float 2 0)))
1126 (hint (math-scale-int (nth 1 two-x
) (nth 2 two-x
)))
1127 (hfrac (math-sub-float x
(math-mul-float (math-float hint
)
1129 (math-mul-float (math-ipow (math-sqrt-e) hint
)
1130 (math-add-float '(float 1 0)
1131 (math-exp-minus-1-raw hfrac
)))))
1132 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x
)))))
1134 (defun math-exp-minus-1-raw (x) ; [F F]
1135 (math-exp-series x
2 3 x x
))
1137 (defun math-exp-series (sum nfac n xpow x
)
1138 (math-working "exp" sum
)
1139 (let* ((nextx (math-mul-float xpow x
))
1140 (nextsum (math-add-float sum
(math-div-float nextx
1141 (math-float nfac
)))))
1142 (if (math-nearly-equal-float sum nextsum
)
1144 (math-exp-series nextsum
(math-mul nfac n
) (1+ n
) nextx x
))))
1150 (defun calcFunc-ln (x) ; [N N] [Public]
1151 (cond ((math-zerop x
)
1152 (if calc-infinite-mode
1153 '(neg (var inf var-inf
))
1154 (math-reject-arg x
"*Logarithm of zero")))
1157 (math-with-extra-prec 2 (math-ln-raw (math-float x
))))
1158 ((eq (car-safe x
) 'sdev
)
1159 (math-make-sdev (calcFunc-ln (nth 1 x
))
1160 (math-div (nth 2 x
) (nth 1 x
))))
1161 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1162 (Math-zerop (nth 2 x
))
1163 (not (math-intv-constp x
))))
1164 (let ((calc-infinite-mode t
))
1165 (math-make-intv (nth 1 x
) (calcFunc-ln (nth 2 x
))
1166 (calcFunc-ln (nth 3 x
)))))
1167 ((equal x
'(var e var-e
))
1169 ((and (eq (car-safe x
) '^
)
1170 (equal (nth 1 x
) '(var e var-e
))
1171 (math-known-realp (nth 2 x
)))
1174 (if (equal x
'(var nan var-nan
))
1176 '(var inf var-inf
)))
1177 (t (calc-record-why 'numberp x
)
1178 (list 'calcFunc-ln x
))))
1180 (defun calcFunc-log10 (x) ; [N N] [Public]
1181 (cond ((math-equal-int x
1)
1182 (if (math-floatp x
) '(float 0 0) 0))
1183 ((and (Math-integerp x
)
1185 (let ((res (math-integer-log x
10)))
1187 (setq x
(cdr res
)))))
1189 ((and (eq (car-safe x
) 'frac
)
1191 (let ((res (math-integer-log (nth 2 x
) 10)))
1193 (setq x
(- (cdr res
))))))
1196 (if calc-infinite-mode
1197 '(neg (var inf var-inf
))
1198 (math-reject-arg x
"*Logarithm of zero")))
1199 (calc-symbolic-mode (signal 'inexact-result nil
))
1201 (math-with-extra-prec 2
1202 (let ((xf (math-float x
)))
1203 (if (eq (nth 1 xf
) 0)
1204 (math-reject-arg x
"*Logarithm of zero"))
1205 (if (Math-integer-posp (nth 1 xf
))
1206 (if (eq (nth 1 xf
) 1) ; log10(1*10^n) = n
1207 (math-float (nth 2 xf
))
1208 (let ((xdigs (1- (math-numdigs (nth 1 xf
)))))
1210 (math-div-float (math-ln-raw-2
1211 (list 'float
(nth 1 xf
) (- xdigs
)))
1213 (math-float (+ (nth 2 xf
) xdigs
)))))
1214 (math-div (calcFunc-ln xf
) (math-ln-10))))))
1215 ((eq (car-safe x
) 'sdev
)
1216 (math-make-sdev (calcFunc-log10 (nth 1 x
))
1218 (math-mul (nth 1 x
) (math-ln-10)))))
1219 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1220 (not (math-intv-constp x
))))
1221 (math-make-intv (nth 1 x
)
1222 (calcFunc-log10 (nth 2 x
))
1223 (calcFunc-log10 (nth 3 x
))))
1225 (if (equal x
'(var nan var-nan
))
1227 '(var inf var-inf
)))
1228 (t (calc-record-why 'numberp x
)
1229 (list 'calcFunc-log10 x
))))
1231 (defun calcFunc-log (x &optional b
) ; [N N N] [Public]
1232 (cond ((or (null b
) (equal b
'(var e var-e
)))
1233 (math-normalize (list 'calcFunc-ln x
)))
1234 ((or (eq b
10) (equal b
'(float 1 1)))
1235 (math-normalize (list 'calcFunc-log10 x
)))
1237 (if calc-infinite-mode
1238 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1239 (math-reject-arg x
"*Logarithm of zero")))
1241 (if calc-infinite-mode
1242 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1243 (math-reject-arg b
"*Logarithm of zero")))
1244 ((math-equal-int b
1)
1245 (if calc-infinite-mode
1246 (math-div (calcFunc-ln x
) 0)
1247 (math-reject-arg b
"*Logarithm base one")))
1248 ((math-equal-int x
1)
1249 (if (or (math-floatp a
) (math-floatp b
)) '(float 0 0) 0))
1250 ((and (Math-ratp x
) (Math-ratp b
)
1251 (math-posp x
) (math-posp b
)
1252 (let* ((sign 1) (inv nil
)
1253 (xx (if (Math-lessp 1 x
)
1257 (bb (if (Math-lessp 1 b
)
1259 (setq sign
(- sign
))
1261 (res (if (Math-lessp xx bb
)
1262 (setq inv
(math-integer-log bb xx
))
1263 (math-integer-log xx bb
))))
1266 (math-div 1 (* sign
(cdr res
)))
1267 (* sign
(cdr res
)))))))
1269 (calc-symbolic-mode (signal 'inexact-result nil
))
1270 ((and (Math-numberp x
) (Math-numberp b
))
1271 (math-with-extra-prec 2
1272 (math-div (math-ln-raw (math-float x
))
1273 (math-log-base-raw b
))))
1274 ((and (eq (car-safe x
) 'sdev
)
1276 (math-make-sdev (calcFunc-log (nth 1 x
) b
)
1279 (math-log-base-raw b
)))))
1280 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1281 (not (math-intv-constp x
)))
1283 (math-make-intv (nth 1 x
)
1284 (calcFunc-log (nth 2 x
) b
)
1285 (calcFunc-log (nth 3 x
) b
)))
1286 ((or (eq (car-safe x
) 'intv
) (eq (car-safe b
) 'intv
))
1287 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1288 ((or (math-infinitep x
)
1290 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1291 (t (if (Math-numberp b
)
1292 (calc-record-why 'numberp x
)
1293 (calc-record-why 'numberp b
))
1294 (list 'calcFunc-log x b
))))
1296 (defun calcFunc-alog (x &optional b
)
1297 (cond ((or (null b
) (equal b
'(var e var-e
)))
1298 (math-normalize (list 'calcFunc-exp x
)))
1299 (t (math-pow b x
))))
1301 (defun calcFunc-ilog (x b
)
1302 (if (and (math-natnump x
) (not (eq x
0))
1303 (math-natnump b
) (not (eq b
0)))
1305 (math-reject-arg x
"*Logarithm base one")
1306 (if (Math-natnum-lessp x b
)
1308 (cdr (math-integer-log x b
))))
1309 (math-floor (calcFunc-log x b
))))
1311 (defun math-integer-log (x b
)
1312 (let ((pows (list b
))
1316 (while (not (Math-lessp x pow
))
1317 (setq pows
(cons pow pows
)
1318 pow
(math-sqr pow
)))
1319 (setq n
(lsh 1 (1- (length pows
)))
1322 (while (and (setq pows
(cdr pows
))
1325 next
(math-mul pow
(car pows
)))
1326 (or (Math-lessp x next
)
1329 (cons (equal pow x
) sum
)))
1332 (defvar math-log-base-cache nil
)
1333 (defun math-log-base-raw (b) ; [N N]
1334 (if (not (and (equal (car math-log-base-cache
) b
)
1335 (eq (nth 1 math-log-base-cache
) calc-internal-prec
)))
1336 (setq math-log-base-cache
(list b calc-internal-prec
1337 (math-ln-raw (math-float b
)))))
1338 (nth 2 math-log-base-cache
))
1340 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1341 (cond ((Math-equal-int x -
1)
1342 (if calc-infinite-mode
1343 '(neg (var inf var-inf
))
1344 (math-reject-arg x
"*Logarithm of zero")))
1346 ((math-zerop x
) '(float 0 0))
1347 (calc-symbolic-mode (signal 'inexact-result nil
))
1349 (math-with-extra-prec 2
1350 (let ((x (math-float x
)))
1351 (if (and (eq (car x
) 'float
)
1352 (math-lessp-float x
'(float 5 -
1))
1353 (math-lessp-float '(float -
5 -
1) x
))
1354 (math-ln-plus-1-raw x
)
1355 (math-ln-raw (math-add-float x
'(float 1 0)))))))
1356 ((eq (car-safe x
) 'sdev
)
1357 (math-make-sdev (calcFunc-lnp1 (nth 1 x
))
1358 (math-div (nth 2 x
) (math-add (nth 1 x
) 1))))
1359 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1360 (not (math-intv-constp x
))))
1361 (math-make-intv (nth 1 x
)
1362 (calcFunc-lnp1 (nth 2 x
))
1363 (calcFunc-lnp1 (nth 3 x
))))
1365 (if (equal x
'(var nan var-nan
))
1367 '(var inf var-inf
)))
1368 (t (calc-record-why 'numberp x
)
1369 (list 'calcFunc-lnp1 x
))))
1371 (defun math-ln-raw (x) ; [N N] --- must be float format!
1372 (cond ((eq (car-safe x
) 'cplx
)
1374 (math-mul-float (math-ln-raw
1375 (math-add-float (math-sqr-float (nth 1 x
))
1376 (math-sqr-float (nth 2 x
))))
1378 (math-arctan2-raw (nth 2 x
) (nth 1 x
))))
1379 ((eq (car x
) 'polar
)
1380 (math-polar (list 'cplx
1381 (math-ln-raw (nth 1 x
))
1382 (math-to-radians (nth 2 x
)))))
1383 ((Math-equal-int x
1)
1385 (calc-symbolic-mode (signal 'inexact-result nil
))
1386 ((math-posp (nth 1 x
)) ; positive and real
1387 (let ((xdigs (1- (math-numdigs (nth 1 x
)))))
1388 (math-add-float (math-ln-raw-2 (list 'float
(nth 1 x
) (- xdigs
)))
1389 (math-mul-float (math-float (+ (nth 2 x
) xdigs
))
1392 (math-reject-arg x
"*Logarithm of zero"))
1393 ((eq calc-complex-mode
'polar
) ; negative and real
1395 (list 'cplx
; negative and real
1396 (math-ln-raw (math-neg-float x
))
1398 (t (list 'cplx
; negative and real
1399 (math-ln-raw (math-neg-float x
))
1402 (defun math-ln-raw-2 (x) ; [F F]
1403 (cond ((math-lessp-float '(float 14 -
1) x
)
1404 (math-add-float (math-ln-raw-2 (math-mul-float x
'(float 5 -
1)))
1406 (t ; now .7 < x <= 1.4
1407 (math-ln-raw-3 (math-div-float (math-sub-float x
'(float 1 0))
1408 (math-add-float x
'(float 1 0)))))))
1410 (defun math-ln-raw-3 (x) ; [F F]
1411 (math-mul-float (math-ln-raw-series x
3 x
(math-sqr-float x
))
1414 ;;; Compute ln((1+x)/(1-x))
1415 (defun math-ln-raw-series (sum n x xsqr
)
1416 (math-working "log" sum
)
1417 (let* ((nextx (math-mul-float x xsqr
))
1418 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1419 (if (math-nearly-equal-float sum nextsum
)
1421 (math-ln-raw-series nextsum
(+ n
2) nextx xsqr
))))
1423 (defun math-ln-plus-1-raw (x)
1424 (math-lnp1-series x
2 x
(math-neg x
)))
1426 (defun math-lnp1-series (sum n xpow x
)
1427 (math-working "lnp1" sum
)
1428 (let* ((nextx (math-mul-float xpow x
))
1429 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1430 (if (math-nearly-equal-float sum nextsum
)
1432 (math-lnp1-series nextsum
(1+ n
) nextx x
))))
1434 (math-defcache math-ln-10
(float (bigpos 018 684 045 994 092 585 302 2) -
21)
1435 (math-ln-raw-2 '(float 1 1)))
1437 (math-defcache math-ln-2
(float (bigpos 417 309 945 559 180 147 693) -
21)
1438 (math-ln-raw-3 (math-float '(frac 1 3))))
1442 ;;; Hyperbolic functions.
1444 (defun calcFunc-sinh (x) ; [N N] [Public]
1446 (math-expand-formulas
1448 (list '/ (list '-
(list 'calcFunc-exp x
)
1449 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1451 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1452 (math-with-extra-prec 2
1453 (let ((expx (math-exp-raw (math-float x
))))
1454 (math-mul (math-add expx
(math-div -
1 expx
)) '(float 5 -
1)))))
1455 ((eq (car-safe x
) 'sdev
)
1456 (math-make-sdev (calcFunc-sinh (nth 1 x
))
1457 (math-mul (nth 2 x
) (calcFunc-cosh (nth 1 x
)))))
1459 (math-sort-intv (nth 1 x
)
1460 (calcFunc-sinh (nth 2 x
))
1461 (calcFunc-sinh (nth 3 x
))))
1462 ((or (equal x
'(var inf var-inf
))
1463 (equal x
'(neg (var inf var-inf
)))
1464 (equal x
'(var nan var-nan
)))
1466 (t (calc-record-why 'numberp x
)
1467 (list 'calcFunc-sinh x
))))
1468 (put 'calcFunc-sinh
'math-expandable t
)
1470 (defun calcFunc-cosh (x) ; [N N] [Public]
1472 (math-expand-formulas
1474 (list '/ (list '+ (list 'calcFunc-exp x
)
1475 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1477 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1478 (math-with-extra-prec 2
1479 (let ((expx (math-exp-raw (math-float x
))))
1480 (math-mul (math-add expx
(math-div 1 expx
)) '(float 5 -
1)))))
1481 ((eq (car-safe x
) 'sdev
)
1482 (math-make-sdev (calcFunc-cosh (nth 1 x
))
1484 (calcFunc-sinh (nth 1 x
)))))
1485 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1486 (setq x
(math-abs x
))
1487 (math-sort-intv (nth 1 x
)
1488 (calcFunc-cosh (nth 2 x
))
1489 (calcFunc-cosh (nth 3 x
))))
1490 ((or (equal x
'(var inf var-inf
))
1491 (equal x
'(neg (var inf var-inf
)))
1492 (equal x
'(var nan var-nan
)))
1494 (t (calc-record-why 'numberp x
)
1495 (list 'calcFunc-cosh x
))))
1496 (put 'calcFunc-cosh
'math-expandable t
)
1498 (defun calcFunc-tanh (x) ; [N N] [Public]
1500 (math-expand-formulas
1502 (let ((expx (list 'calcFunc-exp x
))
1503 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1505 (list '/ (list '- expx expmx
) (list '+ expx expmx
))))))
1507 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1508 (math-with-extra-prec 2
1509 (let* ((expx (calcFunc-exp (math-float x
)))
1510 (expmx (math-div 1 expx
)))
1511 (math-div (math-sub expx expmx
)
1512 (math-add expx expmx
)))))
1513 ((eq (car-safe x
) 'sdev
)
1514 (math-make-sdev (calcFunc-tanh (nth 1 x
))
1516 (math-sqr (calcFunc-cosh (nth 1 x
))))))
1518 (math-sort-intv (nth 1 x
)
1519 (calcFunc-tanh (nth 2 x
))
1520 (calcFunc-tanh (nth 3 x
))))
1521 ((equal x
'(var inf var-inf
))
1523 ((equal x
'(neg (var inf var-inf
)))
1525 ((equal x
'(var nan var-nan
))
1527 (t (calc-record-why 'numberp x
)
1528 (list 'calcFunc-tanh x
))))
1529 (put 'calcFunc-tanh
'math-expandable t
)
1531 (defun calcFunc-arcsinh (x) ; [N N] [Public]
1533 (math-expand-formulas
1535 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
1536 (list '+ (list '^ x
2) 1))))))
1538 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1539 (math-with-extra-prec 2
1540 (math-ln-raw (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
1542 ((eq (car-safe x
) 'sdev
)
1543 (math-make-sdev (calcFunc-arcsinh (nth 1 x
))
1546 (math-add (math-sqr (nth 1 x
)) 1)))))
1548 (math-sort-intv (nth 1 x
)
1549 (calcFunc-arcsinh (nth 2 x
))
1550 (calcFunc-arcsinh (nth 3 x
))))
1551 ((or (equal x
'(var inf var-inf
))
1552 (equal x
'(neg (var inf var-inf
)))
1553 (equal x
'(var nan var-nan
)))
1555 (t (calc-record-why 'numberp x
)
1556 (list 'calcFunc-arcsinh x
))))
1557 (put 'calcFunc-arcsinh
'math-expandable t
)
1559 (defun calcFunc-arccosh (x) ; [N N] [Public]
1561 ((and (eq x -
1) calc-symbolic-mode
)
1563 ((and (eq x
0) calc-symbolic-mode
)
1564 (math-div (math-mul '(var pi var-pi
) '(var i var-i
)) 2))
1565 (math-expand-formulas
1567 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
1568 (list '-
(list '^ x
2) 1))))))
1570 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1571 (if (Math-equal-int x -
1)
1572 (math-imaginary (math-pi))
1573 (math-with-extra-prec 2
1574 (if (or t
; need to do this even in the real case!
1575 (memq (car-safe x
) '(cplx polar
)))
1576 (let ((xp1 (math-add 1 x
))) ; this gets the branch cuts right
1578 (math-add x
(math-mul xp1
1585 (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
1586 '(float -
1 0)))))))))
1587 ((eq (car-safe x
) 'sdev
)
1588 (math-make-sdev (calcFunc-arccosh (nth 1 x
))
1591 (math-add (math-sqr (nth 1 x
)) -
1)))))
1593 (math-sort-intv (nth 1 x
)
1594 (calcFunc-arccosh (nth 2 x
))
1595 (calcFunc-arccosh (nth 3 x
))))
1596 ((or (equal x
'(var inf var-inf
))
1597 (equal x
'(neg (var inf var-inf
)))
1598 (equal x
'(var nan var-nan
)))
1600 (t (calc-record-why 'numberp x
)
1601 (list 'calcFunc-arccosh x
))))
1602 (put 'calcFunc-arccosh
'math-expandable t
)
1604 (defun calcFunc-arctanh (x) ; [N N] [Public]
1606 ((and (Math-equal-int x
1) calc-infinite-mode
)
1608 ((and (Math-equal-int x -
1) calc-infinite-mode
)
1609 '(neg (var inf var-inf
)))
1610 (math-expand-formulas
1612 (list 'calcFunc-ln
(list '+ 1 x
))
1613 (list 'calcFunc-ln
(list '-
1 x
))) 2))
1615 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1616 (math-with-extra-prec 2
1617 (if (or (memq (car-safe x
) '(cplx polar
))
1619 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x
))
1620 (math-ln-raw (math-sub '(float 1 0) x
)))
1622 (if (and (math-equal-int x
1) calc-infinite-mode
)
1624 (if (and (math-equal-int x -
1) calc-infinite-mode
)
1625 '(neg (var inf var-inf
))
1626 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x
)
1629 ((eq (car-safe x
) 'sdev
)
1630 (math-make-sdev (calcFunc-arctanh (nth 1 x
))
1632 (math-sub 1 (math-sqr (nth 1 x
))))))
1634 (math-sort-intv (nth 1 x
)
1635 (calcFunc-arctanh (nth 2 x
))
1636 (calcFunc-arctanh (nth 3 x
))))
1637 ((equal x
'(var nan var-nan
))
1639 (t (calc-record-why 'numberp x
)
1640 (list 'calcFunc-arctanh x
))))
1641 (put 'calcFunc-arctanh
'math-expandable t
)
1644 ;;; Convert A from HMS or degrees to radians.
1645 (defun calcFunc-rad (a) ; [R R] [Public]
1646 (cond ((or (Math-numberp a
)
1648 (math-with-extra-prec 2
1649 (math-mul a
(math-pi-over-180))))
1651 (math-from-hms a
'rad
))
1653 (math-make-sdev (calcFunc-rad (nth 1 a
))
1654 (calcFunc-rad (nth 2 a
))))
1655 (math-expand-formulas
1656 (math-div (math-mul a
'(var pi var-pi
)) 180))
1657 ((math-infinitep a
) a
)
1658 (t (list 'calcFunc-rad a
))))
1659 (put 'calcFunc-rad
'math-expandable t
)
1661 ;;; Convert A from HMS or radians to degrees.
1662 (defun calcFunc-deg (a) ; [R R] [Public]
1663 (cond ((or (Math-numberp a
)
1665 (math-with-extra-prec 2
1666 (math-div a
(math-pi-over-180))))
1668 (math-from-hms a
'deg
))
1670 (math-make-sdev (calcFunc-deg (nth 1 a
))
1671 (calcFunc-deg (nth 2 a
))))
1672 (math-expand-formulas
1673 (math-div (math-mul 180 a
) '(var pi var-pi
)))
1674 ((math-infinitep a
) a
)
1675 (t (list 'calcFunc-deg a
))))
1676 (put 'calcFunc-deg
'math-expandable t
)
1678 ;;; calc-math.el ends here