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: Jay Belanger <belanger@truman.edu>
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-sqrt (arg)
38 (calc-unary-op "^2" 'calcFunc-sqr arg)
39 (calc-unary-op "sqrt" 'calcFunc-sqrt arg))))
41 (defun calc-isqrt (arg)
45 (calc-unary-op "^2" 'calcFunc-sqr arg)
46 (calc-unary-op "isqt" 'calcFunc-isqrt arg))))
49 (defun calc-hypot (arg)
52 (calc-binary-op "hypt" 'calcFunc-hypot arg)))
59 (defun calc-log10 (arg)
61 (calc-hyperbolic-func)
68 (calc-binary-op "alog" 'calcFunc-alog arg)
69 (calc-binary-op "log" 'calcFunc-log arg))))
71 (defun calc-ilog (arg)
75 (calc-binary-op "alog" 'calcFunc-alog arg)
76 (calc-binary-op "ilog" 'calcFunc-ilog arg))))
78 (defun calc-lnp1 (arg)
86 (if (calc-is-hyperbolic)
88 (calc-unary-op "lg10" 'calcFunc-log10 arg)
89 (calc-unary-op "10^" 'calcFunc-exp10 arg))
91 (calc-unary-op "ln" 'calcFunc-ln arg)
92 (calc-unary-op "exp" 'calcFunc-exp arg)))))
94 (defun calc-expm1 (arg)
98 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
99 (calc-unary-op "ex-1" 'calcFunc-expm1 arg))))
104 (if (calc-is-inverse)
105 (if (calc-is-hyperbolic)
106 (if calc-symbolic-mode
107 (calc-pop-push-record 0 "phi" '(var phi var-phi))
108 (calc-pop-push-record 0 "phi" (math-phi)))
109 (if calc-symbolic-mode
110 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma))
111 (calc-pop-push-record 0 "gmma" (math-gamma-const))))
112 (if (calc-is-hyperbolic)
113 (if calc-symbolic-mode
114 (calc-pop-push-record 0 "e" '(var e var-e))
115 (calc-pop-push-record 0 "e" (math-e)))
116 (if calc-symbolic-mode
117 (calc-pop-push-record 0 "pi" '(var pi var-pi))
118 (calc-pop-push-record 0 "pi" (math-pi)))))))
120 (defun calc-sin (arg)
123 (if (calc-is-hyperbolic)
124 (if (calc-is-inverse)
125 (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
126 (calc-unary-op "sinh" 'calcFunc-sinh arg))
127 (if (calc-is-inverse)
128 (calc-unary-op "asin" 'calcFunc-arcsin arg)
129 (calc-unary-op "sin" 'calcFunc-sin arg)))))
131 (defun calc-arcsin (arg)
136 (defun calc-sinh (arg)
138 (calc-hyperbolic-func)
141 (defun calc-arcsinh (arg)
144 (calc-hyperbolic-func)
147 (defun calc-cos (arg)
150 (if (calc-is-hyperbolic)
151 (if (calc-is-inverse)
152 (calc-unary-op "acsh" 'calcFunc-arccosh arg)
153 (calc-unary-op "cosh" 'calcFunc-cosh arg))
154 (if (calc-is-inverse)
155 (calc-unary-op "acos" 'calcFunc-arccos arg)
156 (calc-unary-op "cos" 'calcFunc-cos arg)))))
158 (defun calc-arccos (arg)
163 (defun calc-cosh (arg)
165 (calc-hyperbolic-func)
168 (defun calc-arccosh (arg)
171 (calc-hyperbolic-func)
174 (defun calc-sincos ()
177 (if (calc-is-inverse)
178 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
179 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1))))))
181 (defun calc-tan (arg)
184 (if (calc-is-hyperbolic)
185 (if (calc-is-inverse)
186 (calc-unary-op "atnh" 'calcFunc-arctanh arg)
187 (calc-unary-op "tanh" 'calcFunc-tanh arg))
188 (if (calc-is-inverse)
189 (calc-unary-op "atan" 'calcFunc-arctan arg)
190 (calc-unary-op "tan" 'calcFunc-tan arg)))))
192 (defun calc-arctan (arg)
197 (defun calc-tanh (arg)
199 (calc-hyperbolic-func)
202 (defun calc-arctanh (arg)
205 (calc-hyperbolic-func)
208 (defun calc-arctan2 ()
211 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2)))))
213 (defun calc-conj (arg)
216 (calc-unary-op "conj" 'calcFunc-conj arg)))
218 (defun calc-imaginary ()
221 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
225 (defun calc-to-degrees (arg)
228 (calc-unary-op ">deg" 'calcFunc-deg arg)))
230 (defun calc-to-radians (arg)
233 (calc-unary-op ">rad" 'calcFunc-rad arg)))
236 (defun calc-degrees-mode (arg)
240 (calc-change-mode 'calc-angle-mode 'deg)
241 (message "Angles measured in degrees")))
242 ((= arg 2) (calc-radians-mode))
243 ((= arg 3) (calc-hms-mode))
244 (t (error "Prefix argument out of range"))))
246 (defun calc-radians-mode ()
249 (calc-change-mode 'calc-angle-mode 'rad)
250 (message "Angles measured in radians")))
253 ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public]
254 ;;; This method takes advantage of the fact that Newton's method starting
255 ;;; with an overestimate always works, even using truncating integer division!
256 (defun math-isqrt (a)
257 (cond ((Math-zerop a) a)
258 ((not (math-natnump a))
259 (math-reject-arg a 'natnump))
261 (math-isqrt-small a))
263 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))))
265 (defun calcFunc-isqrt (a)
267 (math-isqrt (math-floor a))
268 (math-floor (math-sqrt a))))
271 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
272 (defun math-isqrt-bignum (a) ; [P.l L]
273 (let ((len (length a)))
275 (let* ((top (nthcdr (- len 2) a)))
276 (math-isqrt-bignum-iter
280 (1+ (math-isqrt-small
281 (+ (* (nth 1 top) 1000) (car top)))))
283 (let* ((top (nth (1- len) a)))
284 (math-isqrt-bignum-iter
287 (list (1+ (math-isqrt-small top)))
290 (defun math-isqrt-bignum-iter (a guess) ; [l L l]
291 (math-working "isqrt" (cons 'bigpos guess))
292 (let* ((q (math-div-bignum a guess))
293 (s (math-add-bignum (car q) guess))
294 (g2 (math-div2-bignum s))
295 (comp (math-compare-bignum g2 guess)))
297 (math-isqrt-bignum-iter a g2)
298 (cons (and (= comp 0)
299 (math-zerop-bignum (cdr q))
303 (defun math-zerop-bignum (a)
306 (while (eq (car (setq a (cdr a))) 0))
309 (defun math-scale-bignum-3 (a n) ; [L L S]
315 (defun math-isqrt-small (a) ; A > 0. [S S]
316 (let ((g (cond ((>= a 10000) 1000)
320 (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
327 ;;; Compute the square root of a number.
328 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
331 (and (Math-zerop a) a)
332 (and (math-known-nonposp a)
333 (math-imaginary (math-sqrt (math-neg a))))
335 (let ((sqrt (math-isqrt-small a)))
336 (if (= (* sqrt sqrt) a)
338 (if calc-symbolic-mode
339 (list 'calcFunc-sqrt a)
340 (math-sqrt-float (math-float a) (math-float sqrt))))))
341 (and (eq (car-safe a) 'bigpos)
342 (let* ((res (math-isqrt-bignum (cdr a)))
343 (sqrt (math-normalize (cons 'bigpos (cdr res)))))
346 (if calc-symbolic-mode
347 (list 'calcFunc-sqrt a)
348 (math-sqrt-float (math-float a) (math-float sqrt))))))
349 (and (eq (car-safe a) 'frac)
350 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a)))))
351 (num-sqrt (math-normalize (cons 'bigpos (cdr num-res))))
352 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a)))))
353 (den-sqrt (math-normalize (cons 'bigpos (cdr den-res)))))
354 (if (and (car num-res) (car den-res))
355 (list 'frac num-sqrt den-sqrt)
356 (if calc-symbolic-mode
357 (if (or (car num-res) (car den-res))
358 (math-div (if (car num-res)
359 num-sqrt (list 'calcFunc-sqrt (nth 1 a)))
361 den-sqrt (list 'calcFunc-sqrt (nth 2 a))))
362 (list 'calcFunc-sqrt a))
363 (math-sqrt-float (math-float a)
364 (math-div (math-float num-sqrt) den-sqrt))))))
365 (and (eq (car-safe a) 'float)
366 (if calc-symbolic-mode
367 (if (= (% (nth 2 a) 2) 0)
368 (let ((res (math-isqrt-bignum
369 (cdr (Math-bignum-test (nth 1 a))))))
371 (math-make-float (math-normalize
372 (cons 'bigpos (cdr res)))
374 (signal 'inexact-result nil)))
375 (signal 'inexact-result nil))
376 (math-sqrt-float a)))
377 (and (eq (car-safe a) 'cplx)
378 (math-with-extra-prec 2
379 (let* ((d (math-abs a))
380 (imag (math-sqrt (math-mul (math-sub d (nth 1 a))
383 (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1)))
384 (if (math-negp (nth 2 a)) (math-neg imag) imag)))))
385 (and (eq (car-safe a) 'polar)
387 (math-sqrt (nth 1 a))
388 (math-mul (nth 2 a) '(float 5 -1))))
389 (and (eq (car-safe a) 'sdev)
390 (let ((sqrt (math-sqrt (nth 1 a))))
392 (math-div (nth 2 a) (math-mul sqrt 2)))))
393 (and (eq (car-safe a) 'intv)
394 (not (math-negp (nth 2 a)))
395 (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a))))
396 (and (eq (car-safe a) '*)
397 (or (math-known-nonnegp (nth 1 a))
398 (math-known-nonnegp (nth 2 a)))
399 (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
400 (and (eq (car-safe a) '/)
401 (or (and (math-known-nonnegp (nth 2 a))
402 (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
403 (and (math-known-nonnegp (nth 1 a))
404 (not (math-equal-int (nth 1 a) 1))
405 (math-mul (math-sqrt (nth 1 a))
406 (math-sqrt (math-div 1 (nth 2 a)))))))
407 (and (eq (car-safe a) '^)
408 (math-known-evenp (nth 2 a))
409 (math-known-realp (nth 1 a))
410 (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2))))
411 (let ((inf (math-infinitep a)))
413 (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
415 (calc-record-why 'numberp a)
416 (list 'calcFunc-sqrt a))))
417 (defalias 'calcFunc-sqrt 'math-sqrt)
419 (defun math-infinite-dir (a &optional inf)
420 (or inf (setq inf (math-infinitep a)))
421 (math-normalize (math-expr-subst a inf 1)))
423 (defun math-sqrt-float (a &optional guess) ; [F F F]
424 (if calc-symbolic-mode
425 (signal 'inexact-result nil)
426 (math-with-extra-prec 1 (math-sqrt-raw a guess))))
428 (defun math-sqrt-raw (a &optional guess) ; [F F F]
429 (if (not (Math-posp a))
432 (let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
433 (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
434 (setq guess (math-make-float (math-isqrt-small
435 (math-scale-int (nth 1 a) (- ldiff)))
436 (/ (+ (nth 2 a) ldiff) 2)))))
437 (math-sqrt-float-iter a guess)))
439 (defun math-sqrt-float-iter (a guess) ; [F F F]
440 (math-working "sqrt" guess)
441 (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
443 (if (math-nearly-equal-float g2 guess)
445 (math-sqrt-float-iter a g2))))
447 ;;; True if A and B differ only in the last digit of precision. [P F F]
448 (defun math-nearly-equal-float (a b)
449 (let ((ediff (- (nth 2 a) (nth 2 b))))
450 (cond ((= ediff 0) ;; Expanded out for speed
451 (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b)))
453 (and (not (consp ediff))
456 (= (math-numdigs (nth 1 a)) calc-internal-prec))))
458 (setq ediff (math-add (Math-integer-neg (nth 1 b))
459 (math-scale-int (nth 1 a) 1)))
460 (and (not (consp ediff))
463 (= (math-numdigs (nth 1 b)) calc-internal-prec)))
465 (setq ediff (math-add (Math-integer-neg (nth 1 a))
466 (math-scale-int (nth 1 b) 1)))
467 (and (not (consp ediff))
470 (= (math-numdigs (nth 1 a)) calc-internal-prec))))))
472 (defun math-nearly-equal (a b) ; [P N N] [Public]
473 (setq a (math-float a))
474 (setq b (math-float b))
475 (if (eq (car a) 'polar) (setq a (math-complex a)))
476 (if (eq (car b) 'polar) (setq b (math-complex b)))
477 (if (eq (car a) 'cplx)
478 (if (eq (car b) 'cplx)
479 (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b))
480 (and (math-nearly-zerop-float (nth 1 a) (nth 2 a))
481 (math-nearly-zerop-float (nth 1 b) (nth 2 b))))
482 (or (math-nearly-equal-float (nth 2 a) (nth 2 b))
483 (and (math-nearly-zerop-float (nth 2 a) (nth 1 a))
484 (math-nearly-zerop-float (nth 2 b) (nth 1 b)))))
485 (and (math-nearly-equal-float (nth 1 a) b)
486 (math-nearly-zerop-float (nth 2 a) b)))
487 (if (eq (car b) 'cplx)
488 (and (math-nearly-equal-float a (nth 1 b))
489 (math-nearly-zerop-float a (nth 2 b)))
490 (math-nearly-equal-float a b))))
492 ;;; True if A is nearly zero compared to B. [P F F]
493 (defun math-nearly-zerop-float (a b)
495 (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
496 (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec)))))
498 (defun math-nearly-zerop (a b) ; [P N R] [Public]
499 (setq a (math-float a))
500 (setq b (math-float b))
501 (if (eq (car a) 'cplx)
502 (and (math-nearly-zerop-float (nth 1 a) b)
503 (math-nearly-zerop-float (nth 2 a) b))
504 (if (eq (car a) 'polar)
505 (math-nearly-zerop-float (nth 1 a) b)
506 (math-nearly-zerop-float a b))))
508 ;;; This implementation could be improved, accuracy-wise.
509 (defun math-hypot (a b)
510 (cond ((Math-zerop a) (math-abs b))
511 ((Math-zerop b) (math-abs a))
512 ((not (Math-scalarp a))
513 (if (math-infinitep a)
514 (if (math-infinitep b)
519 (calc-record-why 'scalarp a)
520 (list 'calcFunc-hypot a b)))
521 ((not (Math-scalarp b))
522 (if (math-infinitep b)
524 (calc-record-why 'scalarp b)
525 (list 'calcFunc-hypot a b)))
526 ((and (Math-numberp a) (Math-numberp b))
527 (math-with-extra-prec 1
528 (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b)))))
529 ((eq (car-safe a) 'hms)
530 (if (eq (car-safe b) 'hms) ; this helps sdev's of hms forms
531 (math-to-hms (math-hypot (math-from-hms a 'deg)
532 (math-from-hms b 'deg)))
533 (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
534 ((eq (car-safe b) 'hms)
535 (math-to-hms (math-hypot a (math-from-hms b 'deg))))
537 (defalias 'calcFunc-hypot 'math-hypot)
539 (defun calcFunc-sqr (x)
544 (defun math-nth-root (a n)
545 (cond ((= n 2) (math-sqrt a))
549 (let ((root (math-nth-root-integer a n)))
552 (and (not calc-symbolic-mode)
553 (math-nth-root-float (math-float a) n
554 (math-float (cdr root)))))))
555 ((eq (car-safe a) 'frac)
556 (let* ((num-root (math-nth-root-integer (nth 1 a) n))
557 (den-root (math-nth-root-integer (nth 2 a) n)))
558 (if (and (car num-root) (car den-root))
559 (list 'frac (cdr num-root) (cdr den-root))
560 (and (not calc-symbolic-mode)
563 (math-div-float (math-float (cdr num-root))
564 (math-float (cdr den-root))))))))
565 ((eq (car-safe a) 'float)
566 (and (not calc-symbolic-mode)
567 (math-nth-root-float a n)))
568 ((eq (car-safe a) 'polar)
569 (let ((root (math-nth-root (nth 1 a) n)))
570 (and root (list 'polar root (math-div (nth 2 a) n)))))
573 ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
574 ;; to math-nth-root-float, but are used by math-nth-root-float-iter,
575 ;; which is called by math-nth-root-float.
578 (defvar math-nrf-nfm1)
580 (defun math-nth-root-float (a math-nrf-n &optional guess)
581 (math-inexact-result)
582 (math-with-extra-prec 1
583 (let ((math-nrf-nf (math-float math-nrf-n))
584 (math-nrf-nfm1 (math-float (1- math-nrf-n))))
585 (math-nth-root-float-iter a (or guess
587 1 (/ (+ (math-numdigs (nth 1 a))
592 (defun math-nth-root-float-iter (a guess)
593 (math-working "root" guess)
594 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess)
596 a (math-ipow guess (1- math-nrf-n))))
598 (if (math-nearly-equal-float g2 guess)
600 (math-nth-root-float-iter a g2))))
602 ;; The variable math-nri-n is local to math-nth-root-integer, but
603 ;; is used by math-nth-root-int-iter, which is called by
604 ;; math-nth-root-int.
607 (defun math-nth-root-integer (a math-nri-n &optional guess) ; [I I S]
608 (math-nth-root-int-iter a (or guess
609 (math-scale-int 1 (/ (+ (math-numdigs a)
613 (defun math-nth-root-int-iter (a guess)
614 (math-working "root" guess)
615 (let* ((q (math-idivmod a (math-ipow guess (1- math-nri-n))))
616 (s (math-add (car q) (math-mul (1- math-nri-n) guess)))
617 (g2 (math-idivmod s math-nri-n)))
618 (if (Math-natnum-lessp (car g2) guess)
619 (math-nth-root-int-iter a (car g2))
620 (cons (and (equal (car g2) guess)
625 (defun calcFunc-nroot (x n)
626 (calcFunc-pow x (if (integerp n)
633 ;;;; Transcendental functions.
635 ;;; All of these functions are defined on the complex plane.
636 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
638 ;;; Most functions increase calc-internal-prec by 2 digits, then round
639 ;;; down afterward. "-raw" functions use the current precision, require
640 ;;; their arguments to be in float (or complex float) format, and always
641 ;;; work in radians (where applicable).
643 (defun math-to-radians (a) ; [N N]
644 (cond ((eq (car-safe a) 'hms)
645 (math-from-hms a 'rad))
646 ((memq calc-angle-mode '(deg hms))
647 (math-mul a (math-pi-over-180)))
650 (defun math-from-radians (a) ; [N N]
651 (cond ((eq calc-angle-mode 'deg)
653 (math-div a (math-pi-over-180))
654 (list 'calcFunc-deg a)))
655 ((eq calc-angle-mode 'hms)
656 (math-to-hms a 'rad))
659 (defun math-to-radians-2 (a) ; [N N]
660 (cond ((eq (car-safe a) 'hms)
661 (math-from-hms a 'rad))
662 ((memq calc-angle-mode '(deg hms))
663 (if calc-symbolic-mode
664 (math-div (math-mul a '(var pi var-pi)) 180)
665 (math-mul a (math-pi-over-180))))
668 (defun math-from-radians-2 (a) ; [N N]
669 (cond ((memq calc-angle-mode '(deg hms))
670 (if calc-symbolic-mode
671 (math-div (math-mul 180 a) '(var pi var-pi))
672 (math-div a (math-pi-over-180))))
677 ;;; Sine, cosine, and tangent.
679 (defun calcFunc-sin (x) ; [N N] [Public]
680 (cond ((and (integerp x)
681 (if (eq calc-angle-mode 'deg)
684 (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
686 (math-with-extra-prec 2
687 (math-sin-raw (math-to-radians (math-float x)))))
690 (math-with-extra-prec 2
691 (let* ((xx (math-to-radians (math-float (nth 1 x))))
692 (xs (math-to-radians (math-float (nth 2 x))))
693 (sc (math-sin-cos-raw xx)))
694 (math-make-sdev (car sc) (math-mul xs (cdr sc)))))
695 (math-make-sdev (calcFunc-sin (nth 1 x))
696 (math-mul (nth 2 x) (calcFunc-cos (nth 1 x))))))
697 ((and (eq (car x) 'intv) (math-intv-constp x))
698 (calcFunc-cos (math-sub x (math-quarter-circle nil))))
699 ((equal x '(var nan var-nan))
701 (t (calc-record-why 'scalarp x)
702 (list 'calcFunc-sin x))))
704 (defun calcFunc-cos (x) ; [N N] [Public]
705 (cond ((and (integerp x)
706 (if (eq calc-angle-mode 'deg)
709 (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
711 (math-with-extra-prec 2
712 (math-cos-raw (math-to-radians (math-float x)))))
715 (math-with-extra-prec 2
716 (let* ((xx (math-to-radians (math-float (nth 1 x))))
717 (xs (math-to-radians (math-float (nth 2 x))))
718 (sc (math-sin-cos-raw xx)))
719 (math-make-sdev (cdr sc) (math-mul xs (car sc)))))
720 (math-make-sdev (calcFunc-cos (nth 1 x))
721 (math-mul (nth 2 x) (calcFunc-sin (nth 1 x))))))
722 ((and (eq (car x) 'intv) (math-intv-constp x))
723 (math-with-extra-prec 2
724 (let* ((xx (math-to-radians (math-float x)))
725 (na (math-floor (math-div (nth 2 xx) (math-pi))))
726 (nb (math-floor (math-div (nth 3 xx) (math-pi))))
727 (span (math-sub nb na)))
728 (if (memq span '(0 1))
729 (let ((int (math-sort-intv (nth 1 x)
730 (math-cos-raw (nth 2 xx))
731 (math-cos-raw (nth 3 xx)))))
734 (math-make-intv (logior (nth 1 x) 2)
737 (math-make-intv (logior (nth 1 x) 1)
741 (list 'intv 3 -1 1)))))
742 ((equal x '(var nan var-nan))
744 (t (calc-record-why 'scalarp x)
745 (list 'calcFunc-cos x))))
747 (defun calcFunc-sincos (x) ; [V N] [Public]
749 (math-with-extra-prec 2
750 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
751 (list 'vec (cdr sc) (car sc)))) ; the vector [cos, sin]
752 (list 'vec (calcFunc-sin x) (calcFunc-cos x))))
754 (defun calcFunc-tan (x) ; [N N] [Public]
755 (cond ((and (integerp x)
756 (if (eq calc-angle-mode 'deg)
761 (math-with-extra-prec 2
762 (math-tan-raw (math-to-radians (math-float x)))))
765 (math-with-extra-prec 2
766 (let* ((xx (math-to-radians (math-float (nth 1 x))))
767 (xs (math-to-radians (math-float (nth 2 x))))
768 (sc (math-sin-cos-raw xx)))
769 (if (and (math-zerop (cdr sc)) (not calc-infinite-mode))
771 (calc-record-why "*Division by zero")
772 (list 'calcFunc-tan x))
773 (math-make-sdev (math-div-float (car sc) (cdr sc))
774 (math-div-float xs (math-sqr (cdr sc)))))))
775 (math-make-sdev (calcFunc-tan (nth 1 x))
777 (math-sqr (calcFunc-cos (nth 1 x)))))))
778 ((and (eq (car x) 'intv) (math-intv-constp x))
779 (or (math-with-extra-prec 2
780 (let* ((xx (math-to-radians (math-float x)))
781 (na (math-floor (math-div (math-sub (nth 2 xx)
784 (nb (math-floor (math-div (math-sub (nth 3 xx)
788 (math-sort-intv (nth 1 x)
789 (math-tan-raw (nth 2 xx))
790 (math-tan-raw (nth 3 xx))))))
791 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
792 ((equal x '(var nan var-nan))
794 (t (calc-record-why 'scalarp x)
795 (list 'calcFunc-tan x))))
797 (defun math-sin-raw (x) ; [N N]
798 (cond ((eq (car x) 'cplx)
799 (let* ((expx (math-exp-raw (nth 2 x)))
800 (expmx (math-div-float '(float 1 0) expx))
801 (sc (math-sin-cos-raw (nth 1 x))))
803 (math-mul-float (car sc)
804 (math-mul-float (math-add-float expx expmx)
806 (math-mul-float (cdr sc)
807 (math-mul-float (math-sub-float expx expmx)
810 (math-polar (math-sin-raw (math-complex x))))
811 ((Math-integer-negp (nth 1 x))
812 (math-neg-float (math-sin-raw (math-neg-float x))))
813 ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff
814 (math-sin-raw (math-mod x (math-two-pi))))
815 (t (math-sin-raw-2 x x))))
817 (defun math-cos-raw (x) ; [N N]
818 (if (eq (car-safe x) 'polar)
819 (math-polar (math-cos-raw (math-complex x)))
820 (math-sin-raw (math-sub (math-pi-over-2) x))))
822 ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
823 ;;; compute either sin(x) or cos(x), whichever is smaller, and compute
824 ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
825 (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
826 (cons (math-sin-raw x) (math-cos-raw x)))
828 (defun math-tan-raw (x) ; [N N]
829 (cond ((eq (car x) 'cplx)
830 (let* ((x (math-mul x '(float 2 0)))
831 (expx (math-exp-raw (nth 2 x)))
832 (expmx (math-div-float '(float 1 0) expx))
833 (sc (math-sin-cos-raw (nth 1 x)))
834 (d (math-add-float (cdr sc)
835 (math-mul-float (math-add-float expx expmx)
837 (and (not (eq (nth 1 d) 0))
839 (math-div-float (car sc) d)
840 (math-div-float (math-mul-float (math-sub-float expx
842 '(float 5 -1)) d)))))
844 (math-polar (math-tan-raw (math-complex x))))
846 (let ((sc (math-sin-cos-raw x)))
847 (if (eq (nth 1 (cdr sc)) 0)
848 (math-div (car sc) 0)
849 (math-div-float (car sc) (cdr sc)))))))
851 (defun math-sin-raw-2 (x orgx) ; This avoids poss of inf recursion. [F F]
852 (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
853 (cond ((Math-integer-negp (nth 1 xmpo2))
854 (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
856 ((math-lessp-float (math-pi-over-4) x)
857 (math-cos-raw-2 xmpo2 orgx))
858 ((math-lessp-float x (math-neg (math-pi-over-4)))
859 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
860 ((math-nearly-zerop-float x orgx) '(float 0 0))
861 (calc-symbolic-mode (signal 'inexact-result nil))
862 (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x)))))))
864 (defun math-cos-raw-2 (x orgx) ; [F F]
865 (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
866 (calc-symbolic-mode (signal 'inexact-result nil))
867 (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
869 (math-add-float '(float 1 0)
870 (math-mul-float xnegsqr '(float 5 -1)))
871 24 5 xnegsqr xnegsqr)))))
873 (defun math-sin-series (sum nfac n x xnegsqr)
874 (math-working "sin" sum)
875 (let* ((nextx (math-mul-float x xnegsqr))
876 (nextsum (math-add-float sum (math-div-float nextx
877 (math-float nfac)))))
878 (if (math-nearly-equal-float sum nextsum)
880 (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
881 (+ n 2) nextx xnegsqr))))
884 ;;; Inverse sine, cosine, tangent.
886 (defun calcFunc-arcsin (x) ; [N N] [Public]
888 ((and (eq x 1) (eq calc-angle-mode 'deg)) 90)
889 ((and (eq x -1) (eq calc-angle-mode 'deg)) -90)
890 (calc-symbolic-mode (signal 'inexact-result nil))
892 (math-with-extra-prec 2
893 (math-from-radians (math-arcsin-raw (math-float x)))))
895 (math-make-sdev (calcFunc-arcsin (nth 1 x))
899 (math-sub 1 (math-sqr (nth 1 x))))))))
901 (math-sort-intv (nth 1 x)
902 (calcFunc-arcsin (nth 2 x))
903 (calcFunc-arcsin (nth 3 x))))
904 ((equal x '(var nan var-nan))
906 (t (calc-record-why 'numberp x)
907 (list 'calcFunc-arcsin x))))
909 (defun calcFunc-arccos (x) ; [N N] [Public]
911 ((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
912 ((and (eq x -1) (eq calc-angle-mode 'deg)) 180)
913 (calc-symbolic-mode (signal 'inexact-result nil))
915 (math-with-extra-prec 2
916 (math-from-radians (math-arccos-raw (math-float x)))))
918 (math-make-sdev (calcFunc-arccos (nth 1 x))
922 (math-sub 1 (math-sqr (nth 1 x))))))))
924 (math-sort-intv (nth 1 x)
925 (calcFunc-arccos (nth 2 x))
926 (calcFunc-arccos (nth 3 x))))
927 ((equal x '(var nan var-nan))
929 (t (calc-record-why 'numberp x)
930 (list 'calcFunc-arccos x))))
932 (defun calcFunc-arctan (x) ; [N N] [Public]
934 ((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
935 ((and (eq x -1) (eq calc-angle-mode 'deg)) -45)
937 (math-with-extra-prec 2
938 (math-from-radians (math-arctan-raw (math-float x)))))
940 (math-make-sdev (calcFunc-arctan (nth 1 x))
943 (math-add 1 (math-sqr (nth 1 x)))))))
945 (math-sort-intv (nth 1 x)
946 (calcFunc-arctan (nth 2 x))
947 (calcFunc-arctan (nth 3 x))))
948 ((equal x '(var inf var-inf))
949 (math-quarter-circle t))
950 ((equal x '(neg (var inf var-inf)))
951 (math-neg (math-quarter-circle t)))
952 ((equal x '(var nan var-nan))
954 (t (calc-record-why 'numberp x)
955 (list 'calcFunc-arctan x))))
957 (defun math-arcsin-raw (x) ; [N N]
958 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
959 (if (or (memq (car x) '(cplx polar))
960 (memq (car a) '(cplx polar)))
961 (math-with-extra-prec 2 ; use extra precision for difficult case
962 (math-mul '(cplx 0 -1)
963 (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
964 (math-arctan2-raw x a))))
966 (defun math-arccos-raw (x) ; [N N]
967 (math-sub (math-pi-over-2) (math-arcsin-raw x)))
969 (defun math-arctan-raw (x) ; [N N]
970 (cond ((memq (car x) '(cplx polar))
971 (math-with-extra-prec 2 ; extra-extra
973 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x)))
974 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x))))
976 ((Math-integer-negp (nth 1 x))
977 (math-neg-float (math-arctan-raw (math-neg-float x))))
979 (calc-symbolic-mode (signal 'inexact-result nil))
980 ((math-equal-int x 1) (math-pi-over-4))
981 ((math-equal-int x -1) (math-neg (math-pi-over-4)))
982 ((math-lessp-float '(float 414214 -6) x) ; if x > sqrt(2) - 1, reduce
983 (if (math-lessp-float '(float 1 0) x)
984 (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
985 (math-arctan-raw (math-div-float '(float 1 0) x)))
986 (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
987 (math-arctan-raw (math-div-float
988 (math-sub-float '(float 1 0) x)
989 (math-add-float '(float 1 0)
991 (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x))))))
993 (defun math-arctan-series (sum n x xnegsqr)
994 (math-working "arctan" sum)
995 (let* ((nextx (math-mul-float x xnegsqr))
996 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
997 (if (math-nearly-equal-float sum nextsum)
999 (math-arctan-series nextsum (+ n 2) nextx xnegsqr))))
1001 (defun calcFunc-arctan2 (y x) ; [F R R] [Public]
1004 (math-with-extra-prec 2
1005 (math-from-radians (math-arctan2-raw (math-float y)
1007 (calc-record-why 'anglep x)
1008 (list 'calcFunc-arctan2 y x))
1009 (if (and (or (math-infinitep x) (math-anglep x))
1010 (or (math-infinitep y) (math-anglep y)))
1025 (calcFunc-arctan2 y x)
1026 '(var nan var-nan)))
1027 (calc-record-why 'anglep y)
1028 (list 'calcFunc-arctan2 y x))))
1030 (defun math-arctan2-raw (y x) ; [F R R]
1031 (cond ((math-zerop y)
1032 (if (math-negp x) (math-pi)
1033 (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))
1037 (math-neg (math-pi-over-2))))
1039 (math-arctan-raw (math-div-float y x)))
1041 (math-add-float (math-arctan-raw (math-div-float y x))
1044 (math-sub-float (math-arctan-raw (math-div-float y x))
1047 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1048 (if (and (Math-vectorp x)
1050 (calcFunc-arctan2 (nth 2 x) (nth 1 x))
1051 (math-reject-arg x "*Two-element vector expected")))
1055 ;;; Exponential function.
1057 (defun calcFunc-exp (x) ; [N N] [Public]
1059 ((and (memq x '(1 -1)) calc-symbolic-mode)
1060 (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e))))
1062 (math-with-extra-prec 2 (math-exp-raw (math-float x))))
1063 ((eq (car-safe x) 'sdev)
1064 (let ((ex (calcFunc-exp (nth 1 x))))
1065 (math-make-sdev ex (math-mul (nth 2 x) ex))))
1066 ((eq (car-safe x) 'intv)
1067 (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x))
1068 (calcFunc-exp (nth 3 x))))
1069 ((equal x '(var inf var-inf))
1071 ((equal x '(neg (var inf var-inf)))
1073 ((equal x '(var nan var-nan))
1075 (t (calc-record-why 'numberp x)
1076 (list 'calcFunc-exp x))))
1078 (defun calcFunc-expm1 (x) ; [N N] [Public]
1080 ((math-zerop x) '(float 0 0))
1081 (calc-symbolic-mode (signal 'inexact-result nil))
1083 (math-with-extra-prec 2
1084 (let ((x (math-float x)))
1085 (if (and (eq (car x) 'float)
1086 (math-lessp-float x '(float 1 0))
1087 (math-lessp-float '(float -1 0) x))
1088 (math-exp-minus-1-raw x)
1089 (math-add (math-exp-raw x) -1)))))
1090 ((eq (car-safe x) 'sdev)
1092 (let ((ex (calcFunc-expm1 (nth 1 x))))
1093 (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
1094 (math-make-sdev (calcFunc-expm1 (nth 1 x))
1095 (math-mul (nth 2 x) (calcFunc-exp (nth 1 x))))))
1096 ((eq (car-safe x) 'intv)
1097 (math-make-intv (nth 1 x)
1098 (calcFunc-expm1 (nth 2 x))
1099 (calcFunc-expm1 (nth 3 x))))
1100 ((equal x '(var inf var-inf))
1102 ((equal x '(neg (var inf var-inf)))
1104 ((equal x '(var nan var-nan))
1106 (t (calc-record-why 'numberp x)
1107 (list 'calcFunc-expm1 x))))
1109 (defun calcFunc-exp10 (x) ; [N N] [Public]
1112 (math-pow '(float 1 1) x)))
1114 (defun math-exp-raw (x) ; [N N]
1115 (cond ((math-zerop x) '(float 1 0))
1116 (calc-symbolic-mode (signal 'inexact-result nil))
1118 (let ((expx (math-exp-raw (nth 1 x)))
1119 (sc (math-sin-cos-raw (nth 2 x))))
1121 (math-mul-float expx (cdr sc))
1122 (math-mul-float expx (car sc)))))
1123 ((eq (car x) 'polar)
1124 (let ((xc (math-complex x)))
1126 (math-exp-raw (nth 1 xc))
1127 (math-from-radians (nth 2 xc)))))
1128 ((or (math-lessp-float '(float 5 -1) x)
1129 (math-lessp-float x '(float -5 -1)))
1130 (if (math-lessp-float '(float 921035 1) x)
1132 (if (math-lessp-float x '(float -921035 1))
1134 (let* ((two-x (math-mul-float x '(float 2 0)))
1135 (hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
1136 (hfrac (math-sub-float x (math-mul-float (math-float hint)
1138 (math-mul-float (math-ipow (math-sqrt-e) hint)
1139 (math-add-float '(float 1 0)
1140 (math-exp-minus-1-raw hfrac)))))
1141 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x)))))
1143 (defun math-exp-minus-1-raw (x) ; [F F]
1144 (math-exp-series x 2 3 x x))
1146 (defun math-exp-series (sum nfac n xpow x)
1147 (math-working "exp" sum)
1148 (let* ((nextx (math-mul-float xpow x))
1149 (nextsum (math-add-float sum (math-div-float nextx
1150 (math-float nfac)))))
1151 (if (math-nearly-equal-float sum nextsum)
1153 (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x))))
1159 (defun calcFunc-ln (x) ; [N N] [Public]
1160 (cond ((math-zerop x)
1161 (if calc-infinite-mode
1162 '(neg (var inf var-inf))
1163 (math-reject-arg x "*Logarithm of zero")))
1166 (math-with-extra-prec 2 (math-ln-raw (math-float x))))
1167 ((eq (car-safe x) 'sdev)
1168 (math-make-sdev (calcFunc-ln (nth 1 x))
1169 (math-div (nth 2 x) (nth 1 x))))
1170 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1171 (Math-zerop (nth 2 x))
1172 (not (math-intv-constp x))))
1173 (let ((calc-infinite-mode t))
1174 (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x))
1175 (calcFunc-ln (nth 3 x)))))
1176 ((equal x '(var e var-e))
1178 ((and (eq (car-safe x) '^)
1179 (equal (nth 1 x) '(var e var-e))
1180 (math-known-realp (nth 2 x)))
1183 (if (equal x '(var nan var-nan))
1185 '(var inf var-inf)))
1186 (t (calc-record-why 'numberp x)
1187 (list 'calcFunc-ln x))))
1189 (defun calcFunc-log10 (x) ; [N N] [Public]
1190 (cond ((math-equal-int x 1)
1191 (if (math-floatp x) '(float 0 0) 0))
1192 ((and (Math-integerp x)
1194 (let ((res (math-integer-log x 10)))
1196 (setq x (cdr res)))))
1198 ((and (eq (car-safe x) 'frac)
1200 (let ((res (math-integer-log (nth 2 x) 10)))
1202 (setq x (- (cdr res))))))
1205 (if calc-infinite-mode
1206 '(neg (var inf var-inf))
1207 (math-reject-arg x "*Logarithm of zero")))
1208 (calc-symbolic-mode (signal 'inexact-result nil))
1210 (math-with-extra-prec 2
1211 (let ((xf (math-float x)))
1212 (if (eq (nth 1 xf) 0)
1213 (math-reject-arg x "*Logarithm of zero"))
1214 (if (Math-integer-posp (nth 1 xf))
1215 (if (eq (nth 1 xf) 1) ; log10(1*10^n) = n
1216 (math-float (nth 2 xf))
1217 (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
1219 (math-div-float (math-ln-raw-2
1220 (list 'float (nth 1 xf) (- xdigs)))
1222 (math-float (+ (nth 2 xf) xdigs)))))
1223 (math-div (calcFunc-ln xf) (math-ln-10))))))
1224 ((eq (car-safe x) 'sdev)
1225 (math-make-sdev (calcFunc-log10 (nth 1 x))
1227 (math-mul (nth 1 x) (math-ln-10)))))
1228 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1229 (not (math-intv-constp x))))
1230 (math-make-intv (nth 1 x)
1231 (calcFunc-log10 (nth 2 x))
1232 (calcFunc-log10 (nth 3 x))))
1234 (if (equal x '(var nan var-nan))
1236 '(var inf var-inf)))
1237 (t (calc-record-why 'numberp x)
1238 (list 'calcFunc-log10 x))))
1240 (defun calcFunc-log (x &optional b) ; [N N N] [Public]
1241 (cond ((or (null b) (equal b '(var e var-e)))
1242 (math-normalize (list 'calcFunc-ln x)))
1243 ((or (eq b 10) (equal b '(float 1 1)))
1244 (math-normalize (list 'calcFunc-log10 x)))
1246 (if calc-infinite-mode
1247 (math-div (calcFunc-ln x) (calcFunc-ln b))
1248 (math-reject-arg x "*Logarithm of zero")))
1250 (if calc-infinite-mode
1251 (math-div (calcFunc-ln x) (calcFunc-ln b))
1252 (math-reject-arg b "*Logarithm of zero")))
1253 ((math-equal-int b 1)
1254 (if calc-infinite-mode
1255 (math-div (calcFunc-ln x) 0)
1256 (math-reject-arg b "*Logarithm base one")))
1257 ((math-equal-int x 1)
1258 (if (math-floatp b) '(float 0 0) 0))
1259 ((and (Math-ratp x) (Math-ratp b)
1260 (math-posp x) (math-posp b)
1261 (let* ((sign 1) (inv nil)
1262 (xx (if (Math-lessp 1 x)
1266 (bb (if (Math-lessp 1 b)
1268 (setq sign (- sign))
1270 (res (if (Math-lessp xx bb)
1271 (setq inv (math-integer-log bb xx))
1272 (math-integer-log xx bb))))
1275 (math-div 1 (* sign (cdr res)))
1276 (* sign (cdr res)))))))
1278 (calc-symbolic-mode (signal 'inexact-result nil))
1279 ((and (Math-numberp x) (Math-numberp b))
1280 (math-with-extra-prec 2
1281 (math-div (math-ln-raw (math-float x))
1282 (math-log-base-raw b))))
1283 ((and (eq (car-safe x) 'sdev)
1285 (math-make-sdev (calcFunc-log (nth 1 x) b)
1288 (math-log-base-raw b)))))
1289 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1290 (not (math-intv-constp x)))
1292 (math-make-intv (nth 1 x)
1293 (calcFunc-log (nth 2 x) b)
1294 (calcFunc-log (nth 3 x) b)))
1295 ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv))
1296 (math-div (calcFunc-ln x) (calcFunc-ln b)))
1297 ((or (math-infinitep x)
1299 (math-div (calcFunc-ln x) (calcFunc-ln b)))
1300 (t (if (Math-numberp b)
1301 (calc-record-why 'numberp x)
1302 (calc-record-why 'numberp b))
1303 (list 'calcFunc-log x b))))
1305 (defun calcFunc-alog (x &optional b)
1306 (cond ((or (null b) (equal b '(var e var-e)))
1307 (math-normalize (list 'calcFunc-exp x)))
1308 (t (math-pow b x))))
1310 (defun calcFunc-ilog (x b)
1311 (if (and (math-natnump x) (not (eq x 0))
1312 (math-natnump b) (not (eq b 0)))
1314 (math-reject-arg x "*Logarithm base one")
1315 (if (Math-natnum-lessp x b)
1317 (cdr (math-integer-log x b))))
1318 (math-floor (calcFunc-log x b))))
1320 (defun math-integer-log (x b)
1321 (let ((pows (list b))
1325 (while (not (Math-lessp x pow))
1326 (setq pows (cons pow pows)
1327 pow (math-sqr pow)))
1328 (setq n (lsh 1 (1- (length pows)))
1331 (while (and (setq pows (cdr pows))
1334 next (math-mul pow (car pows)))
1335 (or (Math-lessp x next)
1338 (cons (equal pow x) sum)))
1341 (defvar math-log-base-cache nil)
1342 (defun math-log-base-raw (b) ; [N N]
1343 (if (not (and (equal (car math-log-base-cache) b)
1344 (eq (nth 1 math-log-base-cache) calc-internal-prec)))
1345 (setq math-log-base-cache (list b calc-internal-prec
1346 (math-ln-raw (math-float b)))))
1347 (nth 2 math-log-base-cache))
1349 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1350 (cond ((Math-equal-int x -1)
1351 (if calc-infinite-mode
1352 '(neg (var inf var-inf))
1353 (math-reject-arg x "*Logarithm of zero")))
1355 ((math-zerop x) '(float 0 0))
1356 (calc-symbolic-mode (signal 'inexact-result nil))
1358 (math-with-extra-prec 2
1359 (let ((x (math-float x)))
1360 (if (and (eq (car x) 'float)
1361 (math-lessp-float x '(float 5 -1))
1362 (math-lessp-float '(float -5 -1) x))
1363 (math-ln-plus-1-raw x)
1364 (math-ln-raw (math-add-float x '(float 1 0)))))))
1365 ((eq (car-safe x) 'sdev)
1366 (math-make-sdev (calcFunc-lnp1 (nth 1 x))
1367 (math-div (nth 2 x) (math-add (nth 1 x) 1))))
1368 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1369 (not (math-intv-constp x))))
1370 (math-make-intv (nth 1 x)
1371 (calcFunc-lnp1 (nth 2 x))
1372 (calcFunc-lnp1 (nth 3 x))))
1374 (if (equal x '(var nan var-nan))
1376 '(var inf var-inf)))
1377 (t (calc-record-why 'numberp x)
1378 (list 'calcFunc-lnp1 x))))
1380 (defun math-ln-raw (x) ; [N N] --- must be float format!
1381 (cond ((eq (car-safe x) 'cplx)
1383 (math-mul-float (math-ln-raw
1384 (math-add-float (math-sqr-float (nth 1 x))
1385 (math-sqr-float (nth 2 x))))
1387 (math-arctan2-raw (nth 2 x) (nth 1 x))))
1388 ((eq (car x) 'polar)
1389 (math-polar (list 'cplx
1390 (math-ln-raw (nth 1 x))
1391 (math-to-radians (nth 2 x)))))
1392 ((Math-equal-int x 1)
1394 (calc-symbolic-mode (signal 'inexact-result nil))
1395 ((math-posp (nth 1 x)) ; positive and real
1396 (let ((xdigs (1- (math-numdigs (nth 1 x)))))
1397 (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
1398 (math-mul-float (math-float (+ (nth 2 x) xdigs))
1401 (math-reject-arg x "*Logarithm of zero"))
1402 ((eq calc-complex-mode 'polar) ; negative and real
1404 (list 'cplx ; negative and real
1405 (math-ln-raw (math-neg-float x))
1407 (t (list 'cplx ; negative and real
1408 (math-ln-raw (math-neg-float x))
1411 (defun math-ln-raw-2 (x) ; [F F]
1412 (cond ((math-lessp-float '(float 14 -1) x)
1413 (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
1415 (t ; now .7 < x <= 1.4
1416 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
1417 (math-add-float x '(float 1 0)))))))
1419 (defun math-ln-raw-3 (x) ; [F F]
1420 (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
1423 ;;; Compute ln((1+x)/(1-x))
1424 (defun math-ln-raw-series (sum n x xsqr)
1425 (math-working "log" sum)
1426 (let* ((nextx (math-mul-float x xsqr))
1427 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1428 (if (math-nearly-equal-float sum nextsum)
1430 (math-ln-raw-series nextsum (+ n 2) nextx xsqr))))
1432 (defun math-ln-plus-1-raw (x)
1433 (math-lnp1-series x 2 x (math-neg x)))
1435 (defun math-lnp1-series (sum n xpow x)
1436 (math-working "lnp1" sum)
1437 (let* ((nextx (math-mul-float xpow x))
1438 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1439 (if (math-nearly-equal-float sum nextsum)
1441 (math-lnp1-series nextsum (1+ n) nextx x))))
1443 (math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21)
1444 (math-ln-raw-2 '(float 1 1)))
1446 (math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21)
1447 (math-ln-raw-3 (math-float '(frac 1 3))))
1451 ;;; Hyperbolic functions.
1453 (defun calcFunc-sinh (x) ; [N N] [Public]
1455 (math-expand-formulas
1457 (list '/ (list '- (list 'calcFunc-exp x)
1458 (list 'calcFunc-exp (list 'neg x))) 2)))
1460 (if calc-symbolic-mode (signal 'inexact-result nil))
1461 (math-with-extra-prec 2
1462 (let ((expx (math-exp-raw (math-float x))))
1463 (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
1464 ((eq (car-safe x) 'sdev)
1465 (math-make-sdev (calcFunc-sinh (nth 1 x))
1466 (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x)))))
1468 (math-sort-intv (nth 1 x)
1469 (calcFunc-sinh (nth 2 x))
1470 (calcFunc-sinh (nth 3 x))))
1471 ((or (equal x '(var inf var-inf))
1472 (equal x '(neg (var inf var-inf)))
1473 (equal x '(var nan var-nan)))
1475 (t (calc-record-why 'numberp x)
1476 (list 'calcFunc-sinh x))))
1477 (put 'calcFunc-sinh 'math-expandable t)
1479 (defun calcFunc-cosh (x) ; [N N] [Public]
1481 (math-expand-formulas
1483 (list '/ (list '+ (list 'calcFunc-exp x)
1484 (list 'calcFunc-exp (list 'neg x))) 2)))
1486 (if calc-symbolic-mode (signal 'inexact-result nil))
1487 (math-with-extra-prec 2
1488 (let ((expx (math-exp-raw (math-float x))))
1489 (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
1490 ((eq (car-safe x) 'sdev)
1491 (math-make-sdev (calcFunc-cosh (nth 1 x))
1493 (calcFunc-sinh (nth 1 x)))))
1494 ((and (eq (car x) 'intv) (math-intv-constp x))
1495 (setq x (math-abs x))
1496 (math-sort-intv (nth 1 x)
1497 (calcFunc-cosh (nth 2 x))
1498 (calcFunc-cosh (nth 3 x))))
1499 ((or (equal x '(var inf var-inf))
1500 (equal x '(neg (var inf var-inf)))
1501 (equal x '(var nan var-nan)))
1503 (t (calc-record-why 'numberp x)
1504 (list 'calcFunc-cosh x))))
1505 (put 'calcFunc-cosh 'math-expandable t)
1507 (defun calcFunc-tanh (x) ; [N N] [Public]
1509 (math-expand-formulas
1511 (let ((expx (list 'calcFunc-exp x))
1512 (expmx (list 'calcFunc-exp (list 'neg x))))
1514 (list '/ (list '- expx expmx) (list '+ expx expmx))))))
1516 (if calc-symbolic-mode (signal 'inexact-result nil))
1517 (math-with-extra-prec 2
1518 (let* ((expx (calcFunc-exp (math-float x)))
1519 (expmx (math-div 1 expx)))
1520 (math-div (math-sub expx expmx)
1521 (math-add expx expmx)))))
1522 ((eq (car-safe x) 'sdev)
1523 (math-make-sdev (calcFunc-tanh (nth 1 x))
1525 (math-sqr (calcFunc-cosh (nth 1 x))))))
1527 (math-sort-intv (nth 1 x)
1528 (calcFunc-tanh (nth 2 x))
1529 (calcFunc-tanh (nth 3 x))))
1530 ((equal x '(var inf var-inf))
1532 ((equal x '(neg (var inf var-inf)))
1534 ((equal x '(var nan var-nan))
1536 (t (calc-record-why 'numberp x)
1537 (list 'calcFunc-tanh x))))
1538 (put 'calcFunc-tanh 'math-expandable t)
1540 (defun calcFunc-arcsinh (x) ; [N N] [Public]
1542 (math-expand-formulas
1544 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
1545 (list '+ (list '^ x 2) 1))))))
1547 (if calc-symbolic-mode (signal 'inexact-result nil))
1548 (math-with-extra-prec 2
1549 (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
1551 ((eq (car-safe x) 'sdev)
1552 (math-make-sdev (calcFunc-arcsinh (nth 1 x))
1555 (math-add (math-sqr (nth 1 x)) 1)))))
1557 (math-sort-intv (nth 1 x)
1558 (calcFunc-arcsinh (nth 2 x))
1559 (calcFunc-arcsinh (nth 3 x))))
1560 ((or (equal x '(var inf var-inf))
1561 (equal x '(neg (var inf var-inf)))
1562 (equal x '(var nan var-nan)))
1564 (t (calc-record-why 'numberp x)
1565 (list 'calcFunc-arcsinh x))))
1566 (put 'calcFunc-arcsinh 'math-expandable t)
1568 (defun calcFunc-arccosh (x) ; [N N] [Public]
1570 ((and (eq x -1) calc-symbolic-mode)
1572 ((and (eq x 0) calc-symbolic-mode)
1573 (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2))
1574 (math-expand-formulas
1576 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
1577 (list '- (list '^ x 2) 1))))))
1579 (if calc-symbolic-mode (signal 'inexact-result nil))
1580 (if (Math-equal-int x -1)
1581 (math-imaginary (math-pi))
1582 (math-with-extra-prec 2
1583 (if (or t ; need to do this even in the real case!
1584 (memq (car-safe x) '(cplx polar)))
1585 (let ((xp1 (math-add 1 x))) ; this gets the branch cuts right
1587 (math-add x (math-mul xp1
1594 (math-add x (math-sqrt-raw (math-add (math-sqr x)
1595 '(float -1 0)))))))))
1596 ((eq (car-safe x) 'sdev)
1597 (math-make-sdev (calcFunc-arccosh (nth 1 x))
1600 (math-add (math-sqr (nth 1 x)) -1)))))
1602 (math-sort-intv (nth 1 x)
1603 (calcFunc-arccosh (nth 2 x))
1604 (calcFunc-arccosh (nth 3 x))))
1605 ((or (equal x '(var inf var-inf))
1606 (equal x '(neg (var inf var-inf)))
1607 (equal x '(var nan var-nan)))
1609 (t (calc-record-why 'numberp x)
1610 (list 'calcFunc-arccosh x))))
1611 (put 'calcFunc-arccosh 'math-expandable t)
1613 (defun calcFunc-arctanh (x) ; [N N] [Public]
1615 ((and (Math-equal-int x 1) calc-infinite-mode)
1617 ((and (Math-equal-int x -1) calc-infinite-mode)
1618 '(neg (var inf var-inf)))
1619 (math-expand-formulas
1621 (list 'calcFunc-ln (list '+ 1 x))
1622 (list 'calcFunc-ln (list '- 1 x))) 2))
1624 (if calc-symbolic-mode (signal 'inexact-result nil))
1625 (math-with-extra-prec 2
1626 (if (or (memq (car-safe x) '(cplx polar))
1628 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x))
1629 (math-ln-raw (math-sub '(float 1 0) x)))
1631 (if (and (math-equal-int x 1) calc-infinite-mode)
1633 (if (and (math-equal-int x -1) calc-infinite-mode)
1634 '(neg (var inf var-inf))
1635 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x)
1638 ((eq (car-safe x) 'sdev)
1639 (math-make-sdev (calcFunc-arctanh (nth 1 x))
1641 (math-sub 1 (math-sqr (nth 1 x))))))
1643 (math-sort-intv (nth 1 x)
1644 (calcFunc-arctanh (nth 2 x))
1645 (calcFunc-arctanh (nth 3 x))))
1646 ((equal x '(var nan var-nan))
1648 (t (calc-record-why 'numberp x)
1649 (list 'calcFunc-arctanh x))))
1650 (put 'calcFunc-arctanh 'math-expandable t)
1653 ;;; Convert A from HMS or degrees to radians.
1654 (defun calcFunc-rad (a) ; [R R] [Public]
1655 (cond ((or (Math-numberp a)
1657 (math-with-extra-prec 2
1658 (math-mul a (math-pi-over-180))))
1660 (math-from-hms a 'rad))
1662 (math-make-sdev (calcFunc-rad (nth 1 a))
1663 (calcFunc-rad (nth 2 a))))
1664 (math-expand-formulas
1665 (math-div (math-mul a '(var pi var-pi)) 180))
1666 ((math-infinitep a) a)
1667 (t (list 'calcFunc-rad a))))
1668 (put 'calcFunc-rad 'math-expandable t)
1670 ;;; Convert A from HMS or radians to degrees.
1671 (defun calcFunc-deg (a) ; [R R] [Public]
1672 (cond ((or (Math-numberp a)
1674 (math-with-extra-prec 2
1675 (math-div a (math-pi-over-180))))
1677 (math-from-hms a 'deg))
1679 (math-make-sdev (calcFunc-deg (nth 1 a))
1680 (calcFunc-deg (nth 2 a))))
1681 (math-expand-formulas
1682 (math-div (math-mul 180 a) '(var pi var-pi)))
1683 ((math-infinitep a) a)
1684 (t (list 'calcFunc-deg a))))
1685 (put 'calcFunc-deg 'math-expandable t)
1687 (provide 'calc-math)
1689 ;;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
1690 ;;; calc-math.el ends here