1;;; calc-math.el --- mathematical functions for Calc 2 3;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004, 4;; 2005, 2006, 2007 Free Software Foundation, Inc. 5 6;; Author: David Gillespie <daveg@synaptics.com> 7;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com> 8 9;; This file is part of GNU Emacs. 10 11;; GNU Emacs is free software; you can redistribute it and/or modify 12;; it under the terms of the GNU General Public License as published by 13;; the Free Software Foundation; either version 2, or (at your option) 14;; any later version. 15 16;; GNU Emacs is distributed in the hope that it will be useful, 17;; but WITHOUT ANY WARRANTY; without even the implied warranty of 18;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19;; GNU General Public License for more details. 20 21;; You should have received a copy of the GNU General Public License 22;; along with GNU Emacs; see the file COPYING. If not, write to the 23;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 24;; Boston, MA 02110-1301, USA. 25 26;;; Commentary: 27 28;;; Code: 29 30;; This file is autoloaded from calc-ext.el. 31 32(require 'calc-ext) 33(require 'calc-macs) 34 35(defun calc-sqrt (arg) 36 (interactive "P") 37 (calc-slow-wrapper 38 (if (calc-is-inverse) 39 (calc-unary-op "^2" 'calcFunc-sqr arg) 40 (calc-unary-op "sqrt" 'calcFunc-sqrt arg)))) 41 42(defun calc-isqrt (arg) 43 (interactive "P") 44 (calc-slow-wrapper 45 (if (calc-is-inverse) 46 (calc-unary-op "^2" 'calcFunc-sqr arg) 47 (calc-unary-op "isqt" 'calcFunc-isqrt arg)))) 48 49 50(defun calc-hypot (arg) 51 (interactive "P") 52 (calc-slow-wrapper 53 (calc-binary-op "hypt" 'calcFunc-hypot arg))) 54 55(defun calc-ln (arg) 56 (interactive "P") 57 (calc-invert-func) 58 (calc-exp arg)) 59 60(defun calc-log10 (arg) 61 (interactive "P") 62 (calc-hyperbolic-func) 63 (calc-ln arg)) 64 65(defun calc-log (arg) 66 (interactive "P") 67 (calc-slow-wrapper 68 (if (calc-is-inverse) 69 (calc-binary-op "alog" 'calcFunc-alog arg) 70 (calc-binary-op "log" 'calcFunc-log arg)))) 71 72(defun calc-ilog (arg) 73 (interactive "P") 74 (calc-slow-wrapper 75 (if (calc-is-inverse) 76 (calc-binary-op "alog" 'calcFunc-alog arg) 77 (calc-binary-op "ilog" 'calcFunc-ilog arg)))) 78 79(defun calc-lnp1 (arg) 80 (interactive "P") 81 (calc-invert-func) 82 (calc-expm1 arg)) 83 84(defun calc-exp (arg) 85 (interactive "P") 86 (calc-slow-wrapper 87 (if (calc-is-hyperbolic) 88 (if (calc-is-inverse) 89 (calc-unary-op "lg10" 'calcFunc-log10 arg) 90 (calc-unary-op "10^" 'calcFunc-exp10 arg)) 91 (if (calc-is-inverse) 92 (calc-unary-op "ln" 'calcFunc-ln arg) 93 (calc-unary-op "exp" 'calcFunc-exp arg))))) 94 95(defun calc-expm1 (arg) 96 (interactive "P") 97 (calc-slow-wrapper 98 (if (calc-is-inverse) 99 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg) 100 (calc-unary-op "ex-1" 'calcFunc-expm1 arg)))) 101 102(defun calc-pi () 103 (interactive) 104 (calc-slow-wrapper 105 (if (calc-is-inverse) 106 (if (calc-is-hyperbolic) 107 (if calc-symbolic-mode 108 (calc-pop-push-record 0 "phi" '(var phi var-phi)) 109 (calc-pop-push-record 0 "phi" (math-phi))) 110 (if calc-symbolic-mode 111 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma)) 112 (calc-pop-push-record 0 "gmma" (math-gamma-const)))) 113 (if (calc-is-hyperbolic) 114 (if calc-symbolic-mode 115 (calc-pop-push-record 0 "e" '(var e var-e)) 116 (calc-pop-push-record 0 "e" (math-e))) 117 (if calc-symbolic-mode 118 (calc-pop-push-record 0 "pi" '(var pi var-pi)) 119 (calc-pop-push-record 0 "pi" (math-pi))))))) 120 121(defun calc-sin (arg) 122 (interactive "P") 123 (calc-slow-wrapper 124 (if (calc-is-hyperbolic) 125 (if (calc-is-inverse) 126 (calc-unary-op "asnh" 'calcFunc-arcsinh arg) 127 (calc-unary-op "sinh" 'calcFunc-sinh arg)) 128 (if (calc-is-inverse) 129 (calc-unary-op "asin" 'calcFunc-arcsin arg) 130 (calc-unary-op "sin" 'calcFunc-sin arg))))) 131 132(defun calc-arcsin (arg) 133 (interactive "P") 134 (calc-invert-func) 135 (calc-sin arg)) 136 137(defun calc-sinh (arg) 138 (interactive "P") 139 (calc-hyperbolic-func) 140 (calc-sin arg)) 141 142(defun calc-arcsinh (arg) 143 (interactive "P") 144 (calc-invert-func) 145 (calc-hyperbolic-func) 146 (calc-sin arg)) 147 148(defun calc-sec (arg) 149 (interactive "P") 150 (calc-slow-wrapper 151 (if (calc-is-hyperbolic) 152 (calc-unary-op "sech" 'calcFunc-sech arg) 153 (calc-unary-op "sec" 'calcFunc-sec arg)))) 154 155(defun calc-sech (arg) 156 (interactive "P") 157 (calc-hyperbolic-func) 158 (calc-sec arg)) 159 160(defun calc-cos (arg) 161 (interactive "P") 162 (calc-slow-wrapper 163 (if (calc-is-hyperbolic) 164 (if (calc-is-inverse) 165 (calc-unary-op "acsh" 'calcFunc-arccosh arg) 166 (calc-unary-op "cosh" 'calcFunc-cosh arg)) 167 (if (calc-is-inverse) 168 (calc-unary-op "acos" 'calcFunc-arccos arg) 169 (calc-unary-op "cos" 'calcFunc-cos arg))))) 170 171(defun calc-arccos (arg) 172 (interactive "P") 173 (calc-invert-func) 174 (calc-cos arg)) 175 176(defun calc-cosh (arg) 177 (interactive "P") 178 (calc-hyperbolic-func) 179 (calc-cos arg)) 180 181(defun calc-arccosh (arg) 182 (interactive "P") 183 (calc-invert-func) 184 (calc-hyperbolic-func) 185 (calc-cos arg)) 186 187(defun calc-csc (arg) 188 (interactive "P") 189 (calc-slow-wrapper 190 (if (calc-is-hyperbolic) 191 (calc-unary-op "csch" 'calcFunc-csch arg) 192 (calc-unary-op "csc" 'calcFunc-csc arg)))) 193 194(defun calc-csch (arg) 195 (interactive "P") 196 (calc-hyperbolic-func) 197 (calc-csc arg)) 198 199(defun calc-sincos () 200 (interactive) 201 (calc-slow-wrapper 202 (if (calc-is-inverse) 203 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1))) 204 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1)))))) 205 206(defun calc-tan (arg) 207 (interactive "P") 208 (calc-slow-wrapper 209 (if (calc-is-hyperbolic) 210 (if (calc-is-inverse) 211 (calc-unary-op "atnh" 'calcFunc-arctanh arg) 212 (calc-unary-op "tanh" 'calcFunc-tanh arg)) 213 (if (calc-is-inverse) 214 (calc-unary-op "atan" 'calcFunc-arctan arg) 215 (calc-unary-op "tan" 'calcFunc-tan arg))))) 216 217(defun calc-arctan (arg) 218 (interactive "P") 219 (calc-invert-func) 220 (calc-tan arg)) 221 222(defun calc-tanh (arg) 223 (interactive "P") 224 (calc-hyperbolic-func) 225 (calc-tan arg)) 226 227(defun calc-arctanh (arg) 228 (interactive "P") 229 (calc-invert-func) 230 (calc-hyperbolic-func) 231 (calc-tan arg)) 232 233(defun calc-cot (arg) 234 (interactive "P") 235 (calc-slow-wrapper 236 (if (calc-is-hyperbolic) 237 (calc-unary-op "coth" 'calcFunc-coth arg) 238 (calc-unary-op "cot" 'calcFunc-cot arg)))) 239 240(defun calc-coth (arg) 241 (interactive "P") 242 (calc-hyperbolic-func) 243 (calc-cot arg)) 244 245(defun calc-arctan2 () 246 (interactive) 247 (calc-slow-wrapper 248 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2))))) 249 250(defun calc-conj (arg) 251 (interactive "P") 252 (calc-wrapper 253 (calc-unary-op "conj" 'calcFunc-conj arg))) 254 255(defun calc-imaginary () 256 (interactive) 257 (calc-slow-wrapper 258 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1))))) 259 260(defun calc-to-degrees (arg) 261 (interactive "P") 262 (calc-wrapper 263 (calc-unary-op ">deg" 'calcFunc-deg arg))) 264 265(defun calc-to-radians (arg) 266 (interactive "P") 267 (calc-wrapper 268 (calc-unary-op ">rad" 'calcFunc-rad arg))) 269 270 271(defun calc-degrees-mode (arg) 272 (interactive "p") 273 (cond ((= arg 1) 274 (calc-wrapper 275 (calc-change-mode 'calc-angle-mode 'deg) 276 (message "Angles measured in degrees"))) 277 ((= arg 2) (calc-radians-mode)) 278 ((= arg 3) (calc-hms-mode)) 279 (t (error "Prefix argument out of range")))) 280 281(defun calc-radians-mode () 282 (interactive) 283 (calc-wrapper 284 (calc-change-mode 'calc-angle-mode 'rad) 285 (message "Angles measured in radians"))) 286 287 288;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public] 289;;; This method takes advantage of the fact that Newton's method starting 290;;; with an overestimate always works, even using truncating integer division! 291(defun math-isqrt (a) 292 (cond ((Math-zerop a) a) 293 ((not (math-natnump a)) 294 (math-reject-arg a 'natnump)) 295 ((integerp a) 296 (math-isqrt-small a)) 297 (t 298 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a)))))))) 299 300(defun calcFunc-isqrt (a) 301 (if (math-realp a) 302 (math-isqrt (math-floor a)) 303 (math-floor (math-sqrt a)))) 304 305 306;;; This returns (flag . result) where the flag is t if A is a perfect square. 307(defun math-isqrt-bignum (a) ; [P.l L] 308 (let ((len (length a))) 309 (if (= (% len 2) 0) 310 (let* ((top (nthcdr (- len 2) a))) 311 (math-isqrt-bignum-iter 312 a 313 (math-scale-bignum-3 314 (math-bignum-big 315 (1+ (math-isqrt-small 316 (+ (* (nth 1 top) 1000) (car top))))) 317 (1- (/ len 2))))) 318 (let* ((top (nth (1- len) a))) 319 (math-isqrt-bignum-iter 320 a 321 (math-scale-bignum-3 322 (list (1+ (math-isqrt-small top))) 323 (/ len 2))))))) 324 325(defun math-isqrt-bignum-iter (a guess) ; [l L l] 326 (math-working "isqrt" (cons 'bigpos guess)) 327 (let* ((q (math-div-bignum a guess)) 328 (s (math-add-bignum (car q) guess)) 329 (g2 (math-div2-bignum s)) 330 (comp (math-compare-bignum g2 guess))) 331 (if (< comp 0) 332 (math-isqrt-bignum-iter a g2) 333 (cons (and (= comp 0) 334 (math-zerop-bignum (cdr q)) 335 (= (% (car s) 2) 0)) 336 guess)))) 337 338(defun math-zerop-bignum (a) 339 (and (eq (car a) 0) 340 (progn 341 (while (eq (car (setq a (cdr a))) 0)) 342 (null a)))) 343 344(defun math-scale-bignum-3 (a n) ; [L L S] 345 (while (> n 0) 346 (setq a (cons 0 a) 347 n (1- n))) 348 a) 349 350(defun math-isqrt-small (a) ; A > 0. [S S] 351 (let ((g (cond ((>= a 10000) 1000) 352 ((>= a 100) 100) 353 (t 10))) 354 g2) 355 (while (< (setq g2 (/ (+ g (/ a g)) 2)) g) 356 (setq g g2)) 357 g)) 358 359 360 361 362;;; Compute the square root of a number. 363;;; [T N] if possible, else [F N] if possible, else [C N]. [Public] 364(defun math-sqrt (a) 365 (or 366 (and (Math-zerop a) a) 367 (and (math-known-nonposp a) 368 (math-imaginary (math-sqrt (math-neg a)))) 369 (and (integerp a) 370 (let ((sqrt (math-isqrt-small a))) 371 (if (= (* sqrt sqrt) a) 372 sqrt 373 (if calc-symbolic-mode 374 (list 'calcFunc-sqrt a) 375 (math-sqrt-float (math-float a) (math-float sqrt)))))) 376 (and (eq (car-safe a) 'bigpos) 377 (let* ((res (math-isqrt-bignum (cdr a))) 378 (sqrt (math-normalize (cons 'bigpos (cdr res))))) 379 (if (car res) 380 sqrt 381 (if calc-symbolic-mode 382 (list 'calcFunc-sqrt a) 383 (math-sqrt-float (math-float a) (math-float sqrt)))))) 384 (and (eq (car-safe a) 'frac) 385 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a))))) 386 (num-sqrt (math-normalize (cons 'bigpos (cdr num-res)))) 387 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a))))) 388 (den-sqrt (math-normalize (cons 'bigpos (cdr den-res))))) 389 (if (and (car num-res) (car den-res)) 390 (list 'frac num-sqrt den-sqrt) 391 (if calc-symbolic-mode 392 (if (or (car num-res) (car den-res)) 393 (math-div (if (car num-res) 394 num-sqrt (list 'calcFunc-sqrt (nth 1 a))) 395 (if (car den-res) 396 den-sqrt (list 'calcFunc-sqrt (nth 2 a)))) 397 (list 'calcFunc-sqrt a)) 398 (math-sqrt-float (math-float a) 399 (math-div (math-float num-sqrt) den-sqrt)))))) 400 (and (eq (car-safe a) 'float) 401 (if calc-symbolic-mode 402 (if (= (% (nth 2 a) 2) 0) 403 (let ((res (math-isqrt-bignum 404 (cdr (Math-bignum-test (nth 1 a)))))) 405 (if (car res) 406 (math-make-float (math-normalize 407 (cons 'bigpos (cdr res))) 408 (/ (nth 2 a) 2)) 409 (signal 'inexact-result nil))) 410 (signal 'inexact-result nil)) 411 (math-sqrt-float a))) 412 (and (eq (car-safe a) 'cplx) 413 (math-with-extra-prec 2 414 (let* ((d (math-abs a)) 415 (imag (math-sqrt (math-mul (math-sub d (nth 1 a)) 416 '(float 5 -1))))) 417 (list 'cplx 418 (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1))) 419 (if (math-negp (nth 2 a)) (math-neg imag) imag))))) 420 (and (eq (car-safe a) 'polar) 421 (list 'polar 422 (math-sqrt (nth 1 a)) 423 (math-mul (nth 2 a) '(float 5 -1)))) 424 (and (eq (car-safe a) 'sdev) 425 (let ((sqrt (math-sqrt (nth 1 a)))) 426 (math-make-sdev sqrt 427 (math-div (nth 2 a) (math-mul sqrt 2))))) 428 (and (eq (car-safe a) 'intv) 429 (not (math-negp (nth 2 a))) 430 (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a)))) 431 (and (eq (car-safe a) '*) 432 (or (math-known-nonnegp (nth 1 a)) 433 (math-known-nonnegp (nth 2 a))) 434 (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a)))) 435 (and (eq (car-safe a) '/) 436 (or (and (math-known-nonnegp (nth 2 a)) 437 (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a)))) 438 (and (math-known-nonnegp (nth 1 a)) 439 (not (math-equal-int (nth 1 a) 1)) 440 (math-mul (math-sqrt (nth 1 a)) 441 (math-sqrt (math-div 1 (nth 2 a))))))) 442 (and (eq (car-safe a) '^) 443 (math-known-evenp (nth 2 a)) 444 (math-known-realp (nth 1 a)) 445 (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2)))) 446 (let ((inf (math-infinitep a))) 447 (and inf 448 (math-mul (math-sqrt (math-infinite-dir a inf)) inf))) 449 (progn 450 (calc-record-why 'numberp a) 451 (list 'calcFunc-sqrt a)))) 452(defalias 'calcFunc-sqrt 'math-sqrt) 453 454(defun math-infinite-dir (a &optional inf) 455 (or inf (setq inf (math-infinitep a))) 456 (math-normalize (math-expr-subst a inf 1))) 457 458(defun math-sqrt-float (a &optional guess) ; [F F F] 459 (if calc-symbolic-mode 460 (signal 'inexact-result nil) 461 (math-with-extra-prec 1 (math-sqrt-raw a guess)))) 462 463(defun math-sqrt-raw (a &optional guess) ; [F F F] 464 (if (not (Math-posp a)) 465 (math-sqrt a) 466 (if (null guess) 467 (let ((ldiff (- (math-numdigs (nth 1 a)) 6))) 468 (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff))) 469 (setq guess (math-make-float (math-isqrt-small 470 (math-scale-int (nth 1 a) (- ldiff))) 471 (/ (+ (nth 2 a) ldiff) 2))))) 472 (math-sqrt-float-iter a guess))) 473 474(defun math-sqrt-float-iter (a guess) ; [F F F] 475 (math-working "sqrt" guess) 476 (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess)) 477 '(float 5 -1)))) 478 (if (math-nearly-equal-float g2 guess) 479 g2 480 (math-sqrt-float-iter a g2)))) 481 482;;; True if A and B differ only in the last digit of precision. [P F F] 483(defun math-nearly-equal-float (a b) 484 (let ((ediff (- (nth 2 a) (nth 2 b)))) 485 (cond ((= ediff 0) ;; Expanded out for speed 486 (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b))) 487 (or (eq ediff 0) 488 (and (not (consp ediff)) 489 (< ediff 10) 490 (> ediff -10) 491 (= (math-numdigs (nth 1 a)) calc-internal-prec)))) 492 ((= ediff 1) 493 (setq ediff (math-add (Math-integer-neg (nth 1 b)) 494 (math-scale-int (nth 1 a) 1))) 495 (and (not (consp ediff)) 496 (< ediff 10) 497 (> ediff -10) 498 (= (math-numdigs (nth 1 b)) calc-internal-prec))) 499 ((= ediff -1) 500 (setq ediff (math-add (Math-integer-neg (nth 1 a)) 501 (math-scale-int (nth 1 b) 1))) 502 (and (not (consp ediff)) 503 (< ediff 10) 504 (> ediff -10) 505 (= (math-numdigs (nth 1 a)) calc-internal-prec)))))) 506 507(defun math-nearly-equal (a b) ; [P N N] [Public] 508 (setq a (math-float a)) 509 (setq b (math-float b)) 510 (if (eq (car a) 'polar) (setq a (math-complex a))) 511 (if (eq (car b) 'polar) (setq b (math-complex b))) 512 (if (eq (car a) 'cplx) 513 (if (eq (car b) 'cplx) 514 (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b)) 515 (and (math-nearly-zerop-float (nth 1 a) (nth 2 a)) 516 (math-nearly-zerop-float (nth 1 b) (nth 2 b)))) 517 (or (math-nearly-equal-float (nth 2 a) (nth 2 b)) 518 (and (math-nearly-zerop-float (nth 2 a) (nth 1 a)) 519 (math-nearly-zerop-float (nth 2 b) (nth 1 b))))) 520 (and (math-nearly-equal-float (nth 1 a) b) 521 (math-nearly-zerop-float (nth 2 a) b))) 522 (if (eq (car b) 'cplx) 523 (and (math-nearly-equal-float a (nth 1 b)) 524 (math-nearly-zerop-float a (nth 2 b))) 525 (math-nearly-equal-float a b)))) 526 527;;; True if A is nearly zero compared to B. [P F F] 528(defun math-nearly-zerop-float (a b) 529 (or (eq (nth 1 a) 0) 530 (<= (+ (math-numdigs (nth 1 a)) (nth 2 a)) 531 (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec))))) 532 533(defun math-nearly-zerop (a b) ; [P N R] [Public] 534 (setq a (math-float a)) 535 (setq b (math-float b)) 536 (if (eq (car a) 'cplx) 537 (and (math-nearly-zerop-float (nth 1 a) b) 538 (math-nearly-zerop-float (nth 2 a) b)) 539 (if (eq (car a) 'polar) 540 (math-nearly-zerop-float (nth 1 a) b) 541 (math-nearly-zerop-float a b)))) 542 543;;; This implementation could be improved, accuracy-wise. 544(defun math-hypot (a b) 545 (cond ((Math-zerop a) (math-abs b)) 546 ((Math-zerop b) (math-abs a)) 547 ((not (Math-scalarp a)) 548 (if (math-infinitep a) 549 (if (math-infinitep b) 550 (if (equal a b) 551 a 552 '(var nan var-nan)) 553 a) 554 (calc-record-why 'scalarp a) 555 (list 'calcFunc-hypot a b))) 556 ((not (Math-scalarp b)) 557 (if (math-infinitep b) 558 b 559 (calc-record-why 'scalarp b) 560 (list 'calcFunc-hypot a b))) 561 ((and (Math-numberp a) (Math-numberp b)) 562 (math-with-extra-prec 1 563 (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b))))) 564 ((eq (car-safe a) 'hms) 565 (if (eq (car-safe b) 'hms) ; this helps sdev's of hms forms 566 (math-to-hms (math-hypot (math-from-hms a 'deg) 567 (math-from-hms b 'deg))) 568 (math-to-hms (math-hypot (math-from-hms a 'deg) b)))) 569 ((eq (car-safe b) 'hms) 570 (math-to-hms (math-hypot a (math-from-hms b 'deg)))) 571 (t nil))) 572(defalias 'calcFunc-hypot 'math-hypot) 573 574(defun calcFunc-sqr (x) 575 (math-pow x 2)) 576 577 578 579(defun math-nth-root (a n) 580 (cond ((= n 2) (math-sqrt a)) 581 ((Math-zerop a) a) 582 ((Math-negp a) nil) 583 ((Math-integerp a) 584 (let ((root (math-nth-root-integer a n))) 585 (if (car root) 586 (cdr root) 587 (and (not calc-symbolic-mode) 588 (math-nth-root-float (math-float a) n 589 (math-float (cdr root))))))) 590 ((eq (car-safe a) 'frac) 591 (let* ((num-root (math-nth-root-integer (nth 1 a) n)) 592 (den-root (math-nth-root-integer (nth 2 a) n))) 593 (if (and (car num-root) (car den-root)) 594 (list 'frac (cdr num-root) (cdr den-root)) 595 (and (not calc-symbolic-mode) 596 (math-nth-root-float 597 (math-float a) n 598 (math-div-float (math-float (cdr num-root)) 599 (math-float (cdr den-root)))))))) 600 ((eq (car-safe a) 'float) 601 (and (not calc-symbolic-mode) 602 (math-nth-root-float a n))) 603 ((eq (car-safe a) 'polar) 604 (let ((root (math-nth-root (nth 1 a) n))) 605 (and root (list 'polar root (math-div (nth 2 a) n))))) 606 (t nil))) 607 608;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local 609;; to math-nth-root-float, but are used by math-nth-root-float-iter, 610;; which is called by math-nth-root-float. 611(defvar math-nrf-n) 612(defvar math-nrf-nf) 613(defvar math-nrf-nfm1) 614 615(defun math-nth-root-float (a math-nrf-n &optional guess) 616 (math-inexact-result) 617 (math-with-extra-prec 1 618 (let ((math-nrf-nf (math-float math-nrf-n)) 619 (math-nrf-nfm1 (math-float (1- math-nrf-n)))) 620 (math-nth-root-float-iter a (or guess 621 (math-make-float 622 1 (/ (+ (math-numdigs (nth 1 a)) 623 (nth 2 a) 624 (/ math-nrf-n 2)) 625 math-nrf-n))))))) 626 627(defun math-nth-root-float-iter (a guess) 628 (math-working "root" guess) 629 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess) 630 (math-div-float 631 a (math-ipow guess (1- math-nrf-n)))) 632 math-nrf-nf))) 633 (if (math-nearly-equal-float g2 guess) 634 g2 635 (math-nth-root-float-iter a g2)))) 636 637;; The variable math-nri-n is local to math-nth-root-integer, but 638;; is used by math-nth-root-int-iter, which is called by 639;; math-nth-root-int. 640(defvar math-nri-n) 641 642(defun math-nth-root-integer (a math-nri-n &optional guess) ; [I I S] 643 (math-nth-root-int-iter a (or guess 644 (math-scale-int 1 (/ (+ (math-numdigs a) 645 (1- math-nri-n)) 646 math-nri-n))))) 647 648(defun math-nth-root-int-iter (a guess) 649 (math-working "root" guess) 650 (let* ((q (math-idivmod a (math-ipow guess (1- math-nri-n)))) 651 (s (math-add (car q) (math-mul (1- math-nri-n) guess))) 652 (g2 (math-idivmod s math-nri-n))) 653 (if (Math-natnum-lessp (car g2) guess) 654 (math-nth-root-int-iter a (car g2)) 655 (cons (and (equal (car g2) guess) 656 (eq (cdr q) 0) 657 (eq (cdr g2) 0)) 658 guess)))) 659 660(defun calcFunc-nroot (x n) 661 (calcFunc-pow x (if (integerp n) 662 (math-make-frac 1 n) 663 (math-div 1 n)))) 664 665 666 667 668;;;; Transcendental functions. 669 670;;; All of these functions are defined on the complex plane. 671;;; (Branch cuts, etc. follow Steele's Common Lisp book.) 672 673;;; Most functions increase calc-internal-prec by 2 digits, then round 674;;; down afterward. "-raw" functions use the current precision, require 675;;; their arguments to be in float (or complex float) format, and always 676;;; work in radians (where applicable). 677 678(defun math-to-radians (a) ; [N N] 679 (cond ((eq (car-safe a) 'hms) 680 (math-from-hms a 'rad)) 681 ((memq calc-angle-mode '(deg hms)) 682 (math-mul a (math-pi-over-180))) 683 (t a))) 684 685(defun math-from-radians (a) ; [N N] 686 (cond ((eq calc-angle-mode 'deg) 687 (if (math-constp a) 688 (math-div a (math-pi-over-180)) 689 (list 'calcFunc-deg a))) 690 ((eq calc-angle-mode 'hms) 691 (math-to-hms a 'rad)) 692 (t a))) 693 694(defun math-to-radians-2 (a) ; [N N] 695 (cond ((eq (car-safe a) 'hms) 696 (math-from-hms a 'rad)) 697 ((memq calc-angle-mode '(deg hms)) 698 (if calc-symbolic-mode 699 (math-div (math-mul a '(var pi var-pi)) 180) 700 (math-mul a (math-pi-over-180)))) 701 (t a))) 702 703(defun math-from-radians-2 (a) ; [N N] 704 (cond ((memq calc-angle-mode '(deg hms)) 705 (if calc-symbolic-mode 706 (math-div (math-mul 180 a) '(var pi var-pi)) 707 (math-div a (math-pi-over-180)))) 708 (t a))) 709 710 711 712;;; Sine, cosine, and tangent. 713 714(defun calcFunc-sin (x) ; [N N] [Public] 715 (cond ((and (integerp x) 716 (if (eq calc-angle-mode 'deg) 717 (= (% x 90) 0) 718 (= x 0))) 719 (aref [0 1 0 -1] (math-mod (/ x 90) 4))) 720 ((Math-scalarp x) 721 (math-with-extra-prec 2 722 (math-sin-raw (math-to-radians (math-float x))))) 723 ((eq (car x) 'sdev) 724 (if (math-constp x) 725 (math-with-extra-prec 2 726 (let* ((xx (math-to-radians (math-float (nth 1 x)))) 727 (xs (math-to-radians (math-float (nth 2 x)))) 728 (sc (math-sin-cos-raw xx))) 729 (math-make-sdev (car sc) (math-mul xs (cdr sc))))) 730 (math-make-sdev (calcFunc-sin (nth 1 x)) 731 (math-mul (nth 2 x) (calcFunc-cos (nth 1 x)))))) 732 ((and (eq (car x) 'intv) (math-intv-constp x)) 733 (calcFunc-cos (math-sub x (math-quarter-circle nil)))) 734 ((equal x '(var nan var-nan)) 735 x) 736 (t (calc-record-why 'scalarp x) 737 (list 'calcFunc-sin x)))) 738 739(defun calcFunc-cos (x) ; [N N] [Public] 740 (cond ((and (integerp x) 741 (if (eq calc-angle-mode 'deg) 742 (= (% x 90) 0) 743 (= x 0))) 744 (aref [1 0 -1 0] (math-mod (/ x 90) 4))) 745 ((Math-scalarp x) 746 (math-with-extra-prec 2 747 (math-cos-raw (math-to-radians (math-float x))))) 748 ((eq (car x) 'sdev) 749 (if (math-constp x) 750 (math-with-extra-prec 2 751 (let* ((xx (math-to-radians (math-float (nth 1 x)))) 752 (xs (math-to-radians (math-float (nth 2 x)))) 753 (sc (math-sin-cos-raw xx))) 754 (math-make-sdev (cdr sc) (math-mul xs (car sc))))) 755 (math-make-sdev (calcFunc-cos (nth 1 x)) 756 (math-mul (nth 2 x) (calcFunc-sin (nth 1 x)))))) 757 ((and (eq (car x) 'intv) (math-intv-constp x)) 758 (math-with-extra-prec 2 759 (let* ((xx (math-to-radians (math-float x))) 760 (na (math-floor (math-div (nth 2 xx) (math-pi)))) 761 (nb (math-floor (math-div (nth 3 xx) (math-pi)))) 762 (span (math-sub nb na))) 763 (if (memq span '(0 1)) 764 (let ((int (math-sort-intv (nth 1 x) 765 (math-cos-raw (nth 2 xx)) 766 (math-cos-raw (nth 3 xx))))) 767 (if (eq span 1) 768 (if (math-evenp na) 769 (math-make-intv (logior (nth 1 x) 2) 770 -1 771 (nth 3 int)) 772 (math-make-intv (logior (nth 1 x) 1) 773 (nth 2 int) 774 1)) 775 int)) 776 (list 'intv 3 -1 1))))) 777 ((equal x '(var nan var-nan)) 778 x) 779 (t (calc-record-why 'scalarp x) 780 (list 'calcFunc-cos x)))) 781 782(defun calcFunc-sincos (x) ; [V N] [Public] 783 (if (Math-scalarp x) 784 (math-with-extra-prec 2 785 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x))))) 786 (list 'vec (cdr sc) (car sc)))) ; the vector [cos, sin] 787 (list 'vec (calcFunc-sin x) (calcFunc-cos x)))) 788 789(defun calcFunc-tan (x) ; [N N] [Public] 790 (cond ((and (integerp x) 791 (if (eq calc-angle-mode 'deg) 792 (= (% x 180) 0) 793 (= x 0))) 794 0) 795 ((Math-scalarp x) 796 (math-with-extra-prec 2 797 (math-tan-raw (math-to-radians (math-float x))))) 798 ((eq (car x) 'sdev) 799 (if (math-constp x) 800 (math-with-extra-prec 2 801 (let* ((xx (math-to-radians (math-float (nth 1 x)))) 802 (xs (math-to-radians (math-float (nth 2 x)))) 803 (sc (math-sin-cos-raw xx))) 804 (if (and (math-zerop (cdr sc)) (not calc-infinite-mode)) 805 (progn 806 (calc-record-why "*Division by zero") 807 (list 'calcFunc-tan x)) 808 (math-make-sdev (math-div-float (car sc) (cdr sc)) 809 (math-div-float xs (math-sqr (cdr sc))))))) 810 (math-make-sdev (calcFunc-tan (nth 1 x)) 811 (math-div (nth 2 x) 812 (math-sqr (calcFunc-cos (nth 1 x))))))) 813 ((and (eq (car x) 'intv) (math-intv-constp x)) 814 (or (math-with-extra-prec 2 815 (let* ((xx (math-to-radians (math-float x))) 816 (na (math-floor (math-div (math-sub (nth 2 xx) 817 (math-pi-over-2)) 818 (math-pi)))) 819 (nb (math-floor (math-div (math-sub (nth 3 xx) 820 (math-pi-over-2)) 821 (math-pi))))) 822 (and (equal na nb) 823 (math-sort-intv (nth 1 x) 824 (math-tan-raw (nth 2 xx)) 825 (math-tan-raw (nth 3 xx)))))) 826 '(intv 3 (neg (var inf var-inf)) (var inf var-inf)))) 827 ((equal x '(var nan var-nan)) 828 x) 829 (t (calc-record-why 'scalarp x) 830 (list 'calcFunc-tan x)))) 831 832(defun calcFunc-sec (x) 833 (cond ((and (integerp x) 834 (eq calc-angle-mode 'deg) 835 (= (% x 180) 0)) 836 (if (= (% x 360) 0) 837 1 838 -1)) 839 ((and (integerp x) 840 (eq calc-angle-mode 'rad) 841 (= x 0)) 842 1) 843 ((Math-scalarp x) 844 (math-with-extra-prec 2 845 (math-sec-raw (math-to-radians (math-float x))))) 846 ((eq (car x) 'sdev) 847 (if (math-constp x) 848 (math-with-extra-prec 2 849 (let* ((xx (math-to-radians (math-float (nth 1 x)))) 850 (xs (math-to-radians (math-float (nth 2 x)))) 851 (sc (math-sin-cos-raw xx))) 852 (if (and (math-zerop (cdr sc)) 853 (not calc-infinite-mode)) 854 (progn 855 (calc-record-why "*Division by zero") 856 (list 'calcFunc-sec x)) 857 (math-make-sdev (math-div-float '(float 1 0) (cdr sc)) 858 (math-div-float 859 (math-mul xs (car sc)) 860 (math-sqr (cdr sc))))))) 861 (math-make-sdev (calcFunc-sec (nth 1 x)) 862 (math-div 863 (math-mul (nth 2 x) 864 (calcFunc-sin (nth 1 x))) 865 (math-sqr (calcFunc-cos (nth 1 x))))))) 866 ((and (eq (car x) 'intv) 867 (math-intv-constp x)) 868 (math-with-extra-prec 2 869 (let* ((xx (math-to-radians (math-float x))) 870 (na (math-floor (math-div (math-sub (nth 2 xx) 871 (math-pi-over-2)) 872 (math-pi)))) 873 (nb (math-floor (math-div (math-sub (nth 3 xx) 874 (math-pi-over-2)) 875 (math-pi)))) 876 (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2)))) 877 (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2)))) 878 (span (math-sub nbb naa))) 879 (if (not (equal na nb)) 880 '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) 881 (let ((int (math-sort-intv (nth 1 x) 882 (math-sec-raw (nth 2 xx)) 883 (math-sec-raw (nth 3 xx))))) 884 (if (eq span 1) 885 (if (math-evenp (math-div (math-add naa 1) 2)) 886 (math-make-intv (logior (nth 1 int) 2) 887 1 888 (nth 3 int)) 889 (math-make-intv (logior (nth 1 int) 1) 890 (nth 2 int) 891 -1)) 892 int)))))) 893 ((equal x '(var nan var-nan)) 894 x) 895 (t (calc-record-why 'scalarp x) 896 (list 'calcFunc-sec x)))) 897 898(defun calcFunc-csc (x) 899 (cond ((and (integerp x) 900 (eq calc-angle-mode 'deg) 901 (= (% (- x 90) 180) 0)) 902 (if (= (% (- x 90) 360) 0) 903 1 904 -1)) 905 ((Math-scalarp x) 906 (math-with-extra-prec 2 907 (math-csc-raw (math-to-radians (math-float x))))) 908 ((eq (car x) 'sdev) 909 (if (math-constp x) 910 (math-with-extra-prec 2 911 (let* ((xx (math-to-radians (math-float (nth 1 x)))) 912 (xs (math-to-radians (math-float (nth 2 x)))) 913 (sc (math-sin-cos-raw xx))) 914 (if (and (math-zerop (car sc)) 915 (not calc-infinite-mode)) 916 (progn 917 (calc-record-why "*Division by zero") 918 (list 'calcFunc-csc x)) 919 (math-make-sdev (math-div-float '(float 1 0) (car sc)) 920 (math-div-float 921 (math-mul xs (cdr sc)) 922 (math-sqr (car sc))))))) 923 (math-make-sdev (calcFunc-csc (nth 1 x)) 924 (math-div 925 (math-mul (nth 2 x) 926 (calcFunc-cos (nth 1 x))) 927 (math-sqr (calcFunc-sin (nth 1 x))))))) 928 ((and (eq (car x) 'intv) 929 (math-intv-constp x)) 930 (math-with-extra-prec 2 931 (let* ((xx (math-to-radians (math-float x))) 932 (na (math-floor (math-div (nth 2 xx) (math-pi)))) 933 (nb (math-floor (math-div (nth 3 xx) (math-pi)))) 934 (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2)))) 935 (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2)))) 936 (span (math-sub nbb naa))) 937 (if (not (equal na nb)) 938 '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) 939 (let ((int (math-sort-intv (nth 1 x) 940 (math-csc-raw (nth 2 xx)) 941 (math-csc-raw (nth 3 xx))))) 942 (if (eq span 1) 943 (if (math-evenp (math-div naa 2)) 944 (math-make-intv (logior (nth 1 int) 2) 945 1 946 (nth 3 int)) 947 (math-make-intv (logior (nth 1 int) 1) 948 (nth 2 int) 949 -1)) 950 int)))))) 951 ((equal x '(var nan var-nan)) 952 x) 953 (t (calc-record-why 'scalarp x) 954 (list 'calcFunc-csc x)))) 955 956(defun calcFunc-cot (x) ; [N N] [Public] 957 (cond ((and (integerp x) 958 (if (eq calc-angle-mode 'deg) 959 (= (% (- x 90) 180) 0) 960 (= x 0))) 961 0) 962 ((Math-scalarp x) 963 (math-with-extra-prec 2 964 (math-cot-raw (math-to-radians (math-float x))))) 965 ((eq (car x) 'sdev) 966 (if (math-constp x) 967 (math-with-extra-prec 2 968 (let* ((xx (math-to-radians (math-float (nth 1 x)))) 969 (xs (math-to-radians (math-float (nth 2 x)))) 970 (sc (math-sin-cos-raw xx))) 971 (if (and (math-zerop (car sc)) (not calc-infinite-mode)) 972 (progn 973 (calc-record-why "*Division by zero") 974 (list 'calcFunc-cot x)) 975 (math-make-sdev (math-div-float (cdr sc) (car sc)) 976 (math-div-float xs (math-sqr (car sc))))))) 977 (math-make-sdev (calcFunc-cot (nth 1 x)) 978 (math-div (nth 2 x) 979 (math-sqr (calcFunc-sin (nth 1 x))))))) 980 ((and (eq (car x) 'intv) (math-intv-constp x)) 981 (or (math-with-extra-prec 2 982 (let* ((xx (math-to-radians (math-float x))) 983 (na (math-floor (math-div (nth 2 xx) (math-pi)))) 984 (nb (math-floor (math-div (nth 3 xx) (math-pi))))) 985 (and (equal na nb) 986 (math-sort-intv (nth 1 x) 987 (math-cot-raw (nth 2 xx)) 988 (math-cot-raw (nth 3 xx)))))) 989 '(intv 3 (neg (var inf var-inf)) (var inf var-inf)))) 990 ((equal x '(var nan var-nan)) 991 x) 992 (t (calc-record-why 'scalarp x) 993 (list 'calcFunc-cot x)))) 994 995(defun math-sin-raw (x) ; [N N] 996 (cond ((eq (car x) 'cplx) 997 (let* ((expx (math-exp-raw (nth 2 x))) 998 (expmx (math-div-float '(float 1 0) expx)) 999 (sc (math-sin-cos-raw (nth 1 x)))) 1000 (list 'cplx 1001 (math-mul-float (car sc) 1002 (math-mul-float (math-add-float expx expmx) 1003 '(float 5 -1))) 1004 (math-mul-float (cdr sc) 1005 (math-mul-float (math-sub-float expx expmx) 1006 '(float 5 -1)))))) 1007 ((eq (car x) 'polar) 1008 (math-polar (math-sin-raw (math-complex x)))) 1009 ((Math-integer-negp (nth 1 x)) 1010 (math-neg-float (math-sin-raw (math-neg-float x)))) 1011 ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff 1012 (math-sin-raw (math-mod x (math-two-pi)))) 1013 (t (math-sin-raw-2 x x)))) 1014 1015(defun math-cos-raw (x) ; [N N] 1016 (if (eq (car-safe x) 'polar) 1017 (math-polar (math-cos-raw (math-complex x))) 1018 (math-sin-raw (math-sub (math-pi-over-2) x)))) 1019 1020(defun math-sec-raw (x) ; [N N] 1021 (cond ((eq (car x) 'cplx) 1022 (let* ((x (math-mul x '(float 1 0))) 1023 (expx (math-exp-raw (nth 2 x))) 1024 (expmx (math-div-float '(float 1 0) expx)) 1025 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1))) 1026 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1))) 1027 (sc (math-sin-cos-raw (nth 1 x))) 1028 (d (math-add-float 1029 (math-mul-float (math-sqr (car sc)) 1030 (math-sqr sh)) 1031 (math-mul-float (math-sqr (cdr sc)) 1032 (math-sqr ch))))) 1033 (and (not (eq (nth 1 d) 0)) 1034 (list 'cplx 1035 (math-div-float (math-mul-float (cdr sc) ch) d) 1036 (math-div-float (math-mul-float (car sc) sh) d))))) 1037 ((eq (car x) 'polar) 1038 (math-polar (math-sec-raw (math-complex x)))) 1039 (t 1040 (let ((cs (math-cos-raw x))) 1041 (if (eq cs 0) 1042 (math-div 1 0) 1043 (math-div-float '(float 1 0) cs)))))) 1044 1045(defun math-csc-raw (x) ; [N N] 1046 (cond ((eq (car x) 'cplx) 1047 (let* ((x (math-mul x '(float 1 0))) 1048 (expx (math-exp-raw (nth 2 x))) 1049 (expmx (math-div-float '(float 1 0) expx)) 1050 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1))) 1051 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1))) 1052 (sc (math-sin-cos-raw (nth 1 x))) 1053 (d (math-add-float 1054 (math-mul-float (math-sqr (car sc)) 1055 (math-sqr ch)) 1056 (math-mul-float (math-sqr (cdr sc)) 1057 (math-sqr sh))))) 1058 (and (not (eq (nth 1 d) 0)) 1059 (list 'cplx 1060 (math-div-float (math-mul-float (car sc) ch) d) 1061 (math-div-float (math-mul-float (cdr sc) sh) d))))) 1062 ((eq (car x) 'polar) 1063 (math-polar (math-csc-raw (math-complex x)))) 1064 (t 1065 (let ((sn (math-sin-raw x))) 1066 (if (eq sn 0) 1067 (math-div 1 0) 1068 (math-div-float '(float 1 0) sn)))))) 1069 1070(defun math-cot-raw (x) ; [N N] 1071 (cond ((eq (car x) 'cplx) 1072 (let* ((x (math-mul x '(float 1 0))) 1073 (expx (math-exp-raw (nth 2 x))) 1074 (expmx (math-div-float '(float 1 0) expx)) 1075 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1))) 1076 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1))) 1077 (sc (math-sin-cos-raw (nth 1 x))) 1078 (d (math-add-float 1079 (math-sqr (car sc)) 1080 (math-sqr sh)))) 1081 (and (not (eq (nth 1 d) 0)) 1082 (list 'cplx 1083 (math-div-float 1084 (math-mul-float (car sc) (cdr sc)) 1085 d) 1086 (math-neg 1087 (math-div-float 1088 (math-mul-float sh ch) 1089 d)))))) 1090 ((eq (car x) 'polar) 1091 (math-polar (math-cot-raw (math-complex x)))) 1092 (t 1093 (let ((sc (math-sin-cos-raw x))) 1094 (if (eq (nth 1 (car sc)) 0) 1095 (math-div (cdr sc) 0) 1096 (math-div-float (cdr sc) (car sc))))))) 1097 1098 1099;;; This could use a smarter method: Reduce x as in math-sin-raw, then 1100;;; compute either sin(x) or cos(x), whichever is smaller, and compute 1101;;; the other using the identity sin(x)^2 + cos(x)^2 = 1. 1102(defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x)) 1103 (cons (math-sin-raw x) (math-cos-raw x))) 1104 1105(defun math-tan-raw (x) ; [N N] 1106 (cond ((eq (car x) 'cplx) 1107 (let* ((x (math-mul x '(float 2 0))) 1108 (expx (math-exp-raw (nth 2 x))) 1109 (expmx (math-div-float '(float 1 0) expx)) 1110 (sc (math-sin-cos-raw (nth 1 x))) 1111 (d (math-add-float (cdr sc) 1112 (math-mul-float (math-add-float expx expmx) 1113 '(float 5 -1))))) 1114 (and (not (eq (nth 1 d) 0)) 1115 (list 'cplx 1116 (math-div-float (car sc) d) 1117 (math-div-float (math-mul-float (math-sub-float expx 1118 expmx) 1119 '(float 5 -1)) d))))) 1120 ((eq (car x) 'polar) 1121 (math-polar (math-tan-raw (math-complex x)))) 1122 (t 1123 (let ((sc (math-sin-cos-raw x))) 1124 (if (eq (nth 1 (cdr sc)) 0) 1125 (math-div (car sc) 0) 1126 (math-div-float (car sc) (cdr sc))))))) 1127 1128(defun math-sin-raw-2 (x orgx) ; This avoids poss of inf recursion. [F F] 1129 (let ((xmpo2 (math-sub-float (math-pi-over-2) x))) 1130 (cond ((Math-integer-negp (nth 1 xmpo2)) 1131 (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi)) 1132 orgx))) 1133 ((math-lessp-float (math-pi-over-4) x) 1134 (math-cos-raw-2 xmpo2 orgx)) 1135 ((math-lessp-float x (math-neg (math-pi-over-4))) 1136 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx))) 1137 ((math-nearly-zerop-float x orgx) '(float 0 0)) 1138 (calc-symbolic-mode (signal 'inexact-result nil)) 1139 (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x))))))) 1140 1141(defun math-cos-raw-2 (x orgx) ; [F F] 1142 (cond ((math-nearly-zerop-float x orgx) '(float 1 0)) 1143 (calc-symbolic-mode (signal 'inexact-result nil)) 1144 (t (let ((xnegsqr (math-neg-float (math-sqr-float x)))) 1145 (math-sin-series 1146 (math-add-float '(float 1 0) 1147 (math-mul-float xnegsqr '(float 5 -1))) 1148 24 5 xnegsqr xnegsqr))))) 1149 1150(defun math-sin-series (sum nfac n x xnegsqr) 1151 (math-working "sin" sum) 1152 (let* ((nextx (math-mul-float x xnegsqr)) 1153 (nextsum (math-add-float sum (math-div-float nextx 1154 (math-float nfac))))) 1155 (if (math-nearly-equal-float sum nextsum) 1156 sum 1157 (math-sin-series nextsum (math-mul nfac (* n (1+ n))) 1158 (+ n 2) nextx xnegsqr)))) 1159 1160 1161;;; Inverse sine, cosine, tangent. 1162 1163(defun calcFunc-arcsin (x) ; [N N] [Public] 1164 (cond ((eq x 0) 0) 1165 ((and (eq x 1) (eq calc-angle-mode 'deg)) 90) 1166 ((and (eq x -1) (eq calc-angle-mode 'deg)) -90) 1167 (calc-symbolic-mode (signal 'inexact-result nil)) 1168 ((Math-numberp x) 1169 (math-with-extra-prec 2 1170 (math-from-radians (math-arcsin-raw (math-float x))))) 1171 ((eq (car x) 'sdev) 1172 (math-make-sdev (calcFunc-arcsin (nth 1 x)) 1173 (math-from-radians 1174 (math-div (nth 2 x) 1175 (math-sqrt 1176 (math-sub 1 (math-sqr (nth 1 x)))))))) 1177 ((eq (car x) 'intv) 1178 (math-sort-intv (nth 1 x) 1179 (calcFunc-arcsin (nth 2 x)) 1180 (calcFunc-arcsin (nth 3 x)))) 1181 ((equal x '(var nan var-nan)) 1182 x) 1183 (t (calc-record-why 'numberp x) 1184 (list 'calcFunc-arcsin x)))) 1185 1186(defun calcFunc-arccos (x) ; [N N] [Public] 1187 (cond ((eq x 1) 0) 1188 ((and (eq x 0) (eq calc-angle-mode 'deg)) 90) 1189 ((and (eq x -1) (eq calc-angle-mode 'deg)) 180) 1190 (calc-symbolic-mode (signal 'inexact-result nil)) 1191 ((Math-numberp x) 1192 (math-with-extra-prec 2 1193 (math-from-radians (math-arccos-raw (math-float x))))) 1194 ((eq (car x) 'sdev) 1195 (math-make-sdev (calcFunc-arccos (nth 1 x)) 1196 (math-from-radians 1197 (math-div (nth 2 x) 1198 (math-sqrt 1199 (math-sub 1 (math-sqr (nth 1 x)))))))) 1200 ((eq (car x) 'intv) 1201 (math-sort-intv (nth 1 x) 1202 (calcFunc-arccos (nth 2 x)) 1203 (calcFunc-arccos (nth 3 x)))) 1204 ((equal x '(var nan var-nan)) 1205 x) 1206 (t (calc-record-why 'numberp x) 1207 (list 'calcFunc-arccos x)))) 1208 1209(defun calcFunc-arctan (x) ; [N N] [Public] 1210 (cond ((eq x 0) 0) 1211 ((and (eq x 1) (eq calc-angle-mode 'deg)) 45) 1212 ((and (eq x -1) (eq calc-angle-mode 'deg)) -45) 1213 ((Math-numberp x) 1214 (math-with-extra-prec 2 1215 (math-from-radians (math-arctan-raw (math-float x))))) 1216 ((eq (car x) 'sdev) 1217 (math-make-sdev (calcFunc-arctan (nth 1 x)) 1218 (math-from-radians 1219 (math-div (nth 2 x) 1220 (math-add 1 (math-sqr (nth 1 x))))))) 1221 ((eq (car x) 'intv) 1222 (math-sort-intv (nth 1 x) 1223 (calcFunc-arctan (nth 2 x)) 1224 (calcFunc-arctan (nth 3 x)))) 1225 ((equal x '(var inf var-inf)) 1226 (math-quarter-circle t)) 1227 ((equal x '(neg (var inf var-inf))) 1228 (math-neg (math-quarter-circle t))) 1229 ((equal x '(var nan var-nan)) 1230 x) 1231 (t (calc-record-why 'numberp x) 1232 (list 'calcFunc-arctan x)))) 1233 1234(defun math-arcsin-raw (x) ; [N N] 1235 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x))))) 1236 (if (or (memq (car x) '(cplx polar)) 1237 (memq (car a) '(cplx polar))) 1238 (math-with-extra-prec 2 ; use extra precision for difficult case 1239 (math-mul '(cplx 0 -1) 1240 (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a)))) 1241 (math-arctan2-raw x a)))) 1242 1243(defun math-arccos-raw (x) ; [N N] 1244 (math-sub (math-pi-over-2) (math-arcsin-raw x))) 1245 1246(defun math-arctan-raw (x) ; [N N] 1247 (cond ((memq (car x) '(cplx polar)) 1248 (math-with-extra-prec 2 ; extra-extra 1249 (math-div (math-sub 1250 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x))) 1251 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x)))) 1252 '(cplx 0 2)))) 1253 ((Math-integer-negp (nth 1 x)) 1254 (math-neg-float (math-arctan-raw (math-neg-float x)))) 1255 ((math-zerop x) x) 1256 (calc-symbolic-mode (signal 'inexact-result nil)) 1257 ((math-equal-int x 1) (math-pi-over-4)) 1258 ((math-equal-int x -1) (math-neg (math-pi-over-4))) 1259 ((math-lessp-float '(float 414214 -6) x) ; if x > sqrt(2) - 1, reduce 1260 (if (math-lessp-float '(float 1 0) x) 1261 (math-sub-float (math-mul-float (math-pi) '(float 5 -1)) 1262 (math-arctan-raw (math-div-float '(float 1 0) x))) 1263 (math-sub-float (math-mul-float (math-pi) '(float 25 -2)) 1264 (math-arctan-raw (math-div-float 1265 (math-sub-float '(float 1 0) x) 1266 (math-add-float '(float 1 0) 1267 x)))))) 1268 (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x)))))) 1269 1270(defun math-arctan-series (sum n x xnegsqr) 1271 (math-working "arctan" sum) 1272 (let* ((nextx (math-mul-float x xnegsqr)) 1273 (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) 1274 (if (math-nearly-equal-float sum nextsum) 1275 sum 1276 (math-arctan-series nextsum (+ n 2) nextx xnegsqr)))) 1277 1278(defun calcFunc-arctan2 (y x) ; [F R R] [Public] 1279 (if (Math-anglep y) 1280 (if (Math-anglep x) 1281 (math-with-extra-prec 2 1282 (math-from-radians (math-arctan2-raw (math-float y) 1283 (math-float x)))) 1284 (calc-record-why 'anglep x) 1285 (list 'calcFunc-arctan2 y x)) 1286 (if (and (or (math-infinitep x) (math-anglep x)) 1287 (or (math-infinitep y) (math-anglep y))) 1288 (progn 1289 (if (math-posp x) 1290 (setq x 1) 1291 (if (math-negp x) 1292 (setq x -1) 1293 (or (math-zerop x) 1294 (setq x nil)))) 1295 (if (math-posp y) 1296 (setq y 1) 1297 (if (math-negp y) 1298 (setq y -1) 1299 (or (math-zerop y) 1300 (setq y nil)))) 1301 (if (and y x) 1302 (calcFunc-arctan2 y x) 1303 '(var nan var-nan))) 1304 (calc-record-why 'anglep y) 1305 (list 'calcFunc-arctan2 y x)))) 1306 1307(defun math-arctan2-raw (y x) ; [F R R] 1308 (cond ((math-zerop y) 1309 (if (math-negp x) (math-pi) 1310 (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0))) 1311 ((math-zerop x) 1312 (if (math-posp y) 1313 (math-pi-over-2) 1314 (math-neg (math-pi-over-2)))) 1315 ((math-posp x) 1316 (math-arctan-raw (math-div-float y x))) 1317 ((math-posp y) 1318 (math-add-float (math-arctan-raw (math-div-float y x)) 1319 (math-pi))) 1320 (t 1321 (math-sub-float (math-arctan-raw (math-div-float y x)) 1322 (math-pi))))) 1323 1324(defun calcFunc-arcsincos (x) ; [V N] [Public] 1325 (if (and (Math-vectorp x) 1326 (= (length x) 3)) 1327 (calcFunc-arctan2 (nth 2 x) (nth 1 x)) 1328 (math-reject-arg x "*Two-element vector expected"))) 1329 1330 1331 1332;;; Exponential function. 1333 1334(defun calcFunc-exp (x) ; [N N] [Public] 1335 (cond ((eq x 0) 1) 1336 ((and (memq x '(1 -1)) calc-symbolic-mode) 1337 (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e)))) 1338 ((Math-numberp x) 1339 (math-with-extra-prec 2 (math-exp-raw (math-float x)))) 1340 ((eq (car-safe x) 'sdev) 1341 (let ((ex (calcFunc-exp (nth 1 x)))) 1342 (math-make-sdev ex (math-mul (nth 2 x) ex)))) 1343 ((eq (car-safe x) 'intv) 1344 (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x)) 1345 (calcFunc-exp (nth 3 x)))) 1346 ((equal x '(var inf var-inf)) 1347 x) 1348 ((equal x '(neg (var inf var-inf))) 1349 0) 1350 ((equal x '(var nan var-nan)) 1351 x) 1352 (t (calc-record-why 'numberp x) 1353 (list 'calcFunc-exp x)))) 1354 1355(defun calcFunc-expm1 (x) ; [N N] [Public] 1356 (cond ((eq x 0) 0) 1357 ((math-zerop x) '(float 0 0)) 1358 (calc-symbolic-mode (signal 'inexact-result nil)) 1359 ((Math-numberp x) 1360 (math-with-extra-prec 2 1361 (let ((x (math-float x))) 1362 (if (and (eq (car x) 'float) 1363 (math-lessp-float x '(float 1 0)) 1364 (math-lessp-float '(float -1 0) x)) 1365 (math-exp-minus-1-raw x) 1366 (math-add (math-exp-raw x) -1))))) 1367 ((eq (car-safe x) 'sdev) 1368 (if (math-constp x) 1369 (let ((ex (calcFunc-expm1 (nth 1 x)))) 1370 (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1)))) 1371 (math-make-sdev (calcFunc-expm1 (nth 1 x)) 1372 (math-mul (nth 2 x) (calcFunc-exp (nth 1 x)))))) 1373 ((eq (car-safe x) 'intv) 1374 (math-make-intv (nth 1 x) 1375 (calcFunc-expm1 (nth 2 x)) 1376 (calcFunc-expm1 (nth 3 x)))) 1377 ((equal x '(var inf var-inf)) 1378 x) 1379 ((equal x '(neg (var inf var-inf))) 1380 -1) 1381 ((equal x '(var nan var-nan)) 1382 x) 1383 (t (calc-record-why 'numberp x) 1384 (list 'calcFunc-expm1 x)))) 1385 1386(defun calcFunc-exp10 (x) ; [N N] [Public] 1387 (if (eq x 0) 1388 1 1389 (math-pow '(float 1 1) x))) 1390 1391(defun math-exp-raw (x) ; [N N] 1392 (cond ((math-zerop x) '(float 1 0)) 1393 (calc-symbolic-mode (signal 'inexact-result nil)) 1394 ((eq (car x) 'cplx) 1395 (let ((expx (math-exp-raw (nth 1 x))) 1396 (sc (math-sin-cos-raw (nth 2 x)))) 1397 (list 'cplx 1398 (math-mul-float expx (cdr sc)) 1399 (math-mul-float expx (car sc))))) 1400 ((eq (car x) 'polar) 1401 (let ((xc (math-complex x))) 1402 (list 'polar 1403 (math-exp-raw (nth 1 xc)) 1404 (math-from-radians (nth 2 xc))))) 1405 ((or (math-lessp-float '(float 5 -1) x) 1406 (math-lessp-float x '(float -5 -1))) 1407 (if (math-lessp-float '(float 921035 1) x) 1408 (math-overflow) 1409 (if (math-lessp-float x '(float -921035 1)) 1410 (math-underflow))) 1411 (let* ((two-x (math-mul-float x '(float 2 0))) 1412 (hint (math-scale-int (nth 1 two-x) (nth 2 two-x))) 1413 (hfrac (math-sub-float x (math-mul-float (math-float hint) 1414 '(float 5 -1))))) 1415 (math-mul-float (math-ipow (math-sqrt-e) hint) 1416 (math-add-float '(float 1 0) 1417 (math-exp-minus-1-raw hfrac))))) 1418 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x))))) 1419 1420(defun math-exp-minus-1-raw (x) ; [F F] 1421 (math-exp-series x 2 3 x x)) 1422 1423(defun math-exp-series (sum nfac n xpow x) 1424 (math-working "exp" sum) 1425 (let* ((nextx (math-mul-float xpow x)) 1426 (nextsum (math-add-float sum (math-div-float nextx 1427 (math-float nfac))))) 1428 (if (math-nearly-equal-float sum nextsum) 1429 sum 1430 (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x)))) 1431 1432 1433 1434;;; Logarithms. 1435 1436(defun calcFunc-ln (x) ; [N N] [Public] 1437 (cond ((math-zerop x) 1438 (if calc-infinite-mode 1439 '(neg (var inf var-inf)) 1440 (math-reject-arg x "*Logarithm of zero"))) 1441 ((eq x 1) 0) 1442 ((Math-numberp x) 1443 (math-with-extra-prec 2 (math-ln-raw (math-float x)))) 1444 ((eq (car-safe x) 'sdev) 1445 (math-make-sdev (calcFunc-ln (nth 1 x)) 1446 (math-div (nth 2 x) (nth 1 x)))) 1447 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) 1448 (Math-zerop (nth 2 x)) 1449 (not (math-intv-constp x)))) 1450 (let ((calc-infinite-mode t)) 1451 (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x)) 1452 (calcFunc-ln (nth 3 x))))) 1453 ((equal x '(var e var-e)) 1454 1) 1455 ((and (eq (car-safe x) '^) 1456 (equal (nth 1 x) '(var e var-e)) 1457 (math-known-realp (nth 2 x))) 1458 (nth 2 x)) 1459 ((math-infinitep x) 1460 (if (equal x '(var nan var-nan)) 1461 x 1462 '(var inf var-inf))) 1463 (t (calc-record-why 'numberp x) 1464 (list 'calcFunc-ln x)))) 1465 1466(defun calcFunc-log10 (x) ; [N N] [Public] 1467 (cond ((math-equal-int x 1) 1468 (if (math-floatp x) '(float 0 0) 0)) 1469 ((and (Math-integerp x) 1470 (math-posp x) 1471 (let ((res (math-integer-log x 10))) 1472 (and (car res) 1473 (setq x (cdr res))))) 1474 x) 1475 ((and (eq (car-safe x) 'frac) 1476 (eq (nth 1 x) 1) 1477 (let ((res (math-integer-log (nth 2 x) 10))) 1478 (and (car res) 1479 (setq x (- (cdr res)))))) 1480 x) 1481 ((math-zerop x) 1482 (if calc-infinite-mode 1483 '(neg (var inf var-inf)) 1484 (math-reject-arg x "*Logarithm of zero"))) 1485 (calc-symbolic-mode (signal 'inexact-result nil)) 1486 ((Math-numberp x) 1487 (math-with-extra-prec 2 1488 (let ((xf (math-float x))) 1489 (if (eq (nth 1 xf) 0) 1490 (math-reject-arg x "*Logarithm of zero")) 1491 (if (Math-integer-posp (nth 1 xf)) 1492 (if (eq (nth 1 xf) 1) ; log10(1*10^n) = n 1493 (math-float (nth 2 xf)) 1494 (let ((xdigs (1- (math-numdigs (nth 1 xf))))) 1495 (math-add-float 1496 (math-div-float (math-ln-raw-2 1497 (list 'float (nth 1 xf) (- xdigs))) 1498 (math-ln-10)) 1499 (math-float (+ (nth 2 xf) xdigs))))) 1500 (math-div (calcFunc-ln xf) (math-ln-10)))))) 1501 ((eq (car-safe x) 'sdev) 1502 (math-make-sdev (calcFunc-log10 (nth 1 x)) 1503 (math-div (nth 2 x) 1504 (math-mul (nth 1 x) (math-ln-10))))) 1505 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) 1506 (not (math-intv-constp x)))) 1507 (math-make-intv (nth 1 x) 1508 (calcFunc-log10 (nth 2 x)) 1509 (calcFunc-log10 (nth 3 x)))) 1510 ((math-infinitep x) 1511 (if (equal x '(var nan var-nan)) 1512 x 1513 '(var inf var-inf))) 1514 (t (calc-record-why 'numberp x) 1515 (list 'calcFunc-log10 x)))) 1516 1517(defun calcFunc-log (x &optional b) ; [N N N] [Public] 1518 (cond ((or (null b) (equal b '(var e var-e))) 1519 (math-normalize (list 'calcFunc-ln x))) 1520 ((or (eq b 10) (equal b '(float 1 1))) 1521 (math-normalize (list 'calcFunc-log10 x))) 1522 ((math-zerop x) 1523 (if calc-infinite-mode 1524 (math-div (calcFunc-ln x) (calcFunc-ln b)) 1525 (math-reject-arg x "*Logarithm of zero"))) 1526 ((math-zerop b) 1527 (if calc-infinite-mode 1528 (math-div (calcFunc-ln x) (calcFunc-ln b)) 1529 (math-reject-arg b "*Logarithm of zero"))) 1530 ((math-equal-int b 1) 1531 (if calc-infinite-mode 1532 (math-div (calcFunc-ln x) 0) 1533 (math-reject-arg b "*Logarithm base one"))) 1534 ((math-equal-int x 1) 1535 (if (math-floatp b) '(float 0 0) 0)) 1536 ((and (Math-ratp x) (Math-ratp b) 1537 (math-posp x) (math-posp b) 1538 (let* ((sign 1) (inv nil) 1539 (xx (if (Math-lessp 1 x) 1540 x 1541 (setq sign -1) 1542 (math-div 1 x))) 1543 (bb (if (Math-lessp 1 b) 1544 b 1545 (setq sign (- sign)) 1546 (math-div 1 b))) 1547 (res (if (Math-lessp xx bb) 1548 (setq inv (math-integer-log bb xx)) 1549 (math-integer-log xx bb)))) 1550 (and (car res) 1551 (setq x (if inv 1552 (math-div 1 (* sign (cdr res))) 1553 (* sign (cdr res))))))) 1554 x) 1555 (calc-symbolic-mode (signal 'inexact-result nil)) 1556 ((and (Math-numberp x) (Math-numberp b)) 1557 (math-with-extra-prec 2 1558 (math-div (math-ln-raw (math-float x)) 1559 (math-log-base-raw b)))) 1560 ((and (eq (car-safe x) 'sdev) 1561 (Math-numberp b)) 1562 (math-make-sdev (calcFunc-log (nth 1 x) b) 1563 (math-div (nth 2 x) 1564 (math-mul (nth 1 x) 1565 (math-log-base-raw b))))) 1566 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) 1567 (not (math-intv-constp x))) 1568 (math-realp b)) 1569 (math-make-intv (nth 1 x) 1570 (calcFunc-log (nth 2 x) b) 1571 (calcFunc-log (nth 3 x) b))) 1572 ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv)) 1573 (math-div (calcFunc-ln x) (calcFunc-ln b))) 1574 ((or (math-infinitep x) 1575 (math-infinitep b)) 1576 (math-div (calcFunc-ln x) (calcFunc-ln b))) 1577 (t (if (Math-numberp b) 1578 (calc-record-why 'numberp x) 1579 (calc-record-why 'numberp b)) 1580 (list 'calcFunc-log x b)))) 1581 1582(defun calcFunc-alog (x &optional b) 1583 (cond ((or (null b) (equal b '(var e var-e))) 1584 (math-normalize (list 'calcFunc-exp x))) 1585 (t (math-pow b x)))) 1586 1587(defun calcFunc-ilog (x b) 1588 (if (and (math-natnump x) (not (eq x 0)) 1589 (math-natnump b) (not (eq b 0))) 1590 (if (eq b 1) 1591 (math-reject-arg x "*Logarithm base one") 1592 (if (Math-natnum-lessp x b) 1593 0 1594 (cdr (math-integer-log x b)))) 1595 (math-floor (calcFunc-log x b)))) 1596 1597(defun math-integer-log (x b) 1598 (let ((pows (list b)) 1599 (pow (math-sqr b)) 1600 next 1601 sum n) 1602 (while (not (Math-lessp x pow)) 1603 (setq pows (cons pow pows) 1604 pow (math-sqr pow))) 1605 (setq n (lsh 1 (1- (length pows))) 1606 sum n 1607 pow (car pows)) 1608 (while (and (setq pows (cdr pows)) 1609 (Math-lessp pow x)) 1610 (setq n (/ n 2) 1611 next (math-mul pow (car pows))) 1612 (or (Math-lessp x next) 1613 (setq pow next 1614 sum (+ sum n)))) 1615 (cons (equal pow x) sum))) 1616 1617 1618(defvar math-log-base-cache nil) 1619(defun math-log-base-raw (b) ; [N N] 1620 (if (not (and (equal (car math-log-base-cache) b) 1621 (eq (nth 1 math-log-base-cache) calc-internal-prec))) 1622 (setq math-log-base-cache (list b calc-internal-prec 1623 (math-ln-raw (math-float b))))) 1624 (nth 2 math-log-base-cache)) 1625 1626(defun calcFunc-lnp1 (x) ; [N N] [Public] 1627 (cond ((Math-equal-int x -1) 1628 (if calc-infinite-mode 1629 '(neg (var inf var-inf)) 1630 (math-reject-arg x "*Logarithm of zero"))) 1631 ((eq x 0) 0) 1632 ((math-zerop x) '(float 0 0)) 1633 (calc-symbolic-mode (signal 'inexact-result nil)) 1634 ((Math-numberp x) 1635 (math-with-extra-prec 2 1636 (let ((x (math-float x))) 1637 (if (and (eq (car x) 'float) 1638 (math-lessp-float x '(float 5 -1)) 1639 (math-lessp-float '(float -5 -1) x)) 1640 (math-ln-plus-1-raw x) 1641 (math-ln-raw (math-add-float x '(float 1 0))))))) 1642 ((eq (car-safe x) 'sdev) 1643 (math-make-sdev (calcFunc-lnp1 (nth 1 x)) 1644 (math-div (nth 2 x) (math-add (nth 1 x) 1)))) 1645 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) 1646 (not (math-intv-constp x)))) 1647 (math-make-intv (nth 1 x) 1648 (calcFunc-lnp1 (nth 2 x)) 1649 (calcFunc-lnp1 (nth 3 x)))) 1650 ((math-infinitep x) 1651 (if (equal x '(var nan var-nan)) 1652 x 1653 '(var inf var-inf))) 1654 (t (calc-record-why 'numberp x) 1655 (list 'calcFunc-lnp1 x)))) 1656 1657(defun math-ln-raw (x) ; [N N] --- must be float format! 1658 (cond ((eq (car-safe x) 'cplx) 1659 (list 'cplx 1660 (math-mul-float (math-ln-raw 1661 (math-add-float (math-sqr-float (nth 1 x)) 1662 (math-sqr-float (nth 2 x)))) 1663 '(float 5 -1)) 1664 (math-arctan2-raw (nth 2 x) (nth 1 x)))) 1665 ((eq (car x) 'polar) 1666 (math-polar (list 'cplx 1667 (math-ln-raw (nth 1 x)) 1668 (math-to-radians (nth 2 x))))) 1669 ((Math-equal-int x 1) 1670 '(float 0 0)) 1671 (calc-symbolic-mode (signal 'inexact-result nil)) 1672 ((math-posp (nth 1 x)) ; positive and real 1673 (let ((xdigs (1- (math-numdigs (nth 1 x))))) 1674 (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs))) 1675 (math-mul-float (math-float (+ (nth 2 x) xdigs)) 1676 (math-ln-10))))) 1677 ((math-zerop x) 1678 (math-reject-arg x "*Logarithm of zero")) 1679 ((eq calc-complex-mode 'polar) ; negative and real 1680 (math-polar 1681 (list 'cplx ; negative and real 1682 (math-ln-raw (math-neg-float x)) 1683 (math-pi)))) 1684 (t (list 'cplx ; negative and real 1685 (math-ln-raw (math-neg-float x)) 1686 (math-pi))))) 1687 1688(defun math-ln-raw-2 (x) ; [F F] 1689 (cond ((math-lessp-float '(float 14 -1) x) 1690 (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1))) 1691 (math-ln-2))) 1692 (t ; now .7 < x <= 1.4 1693 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0)) 1694 (math-add-float x '(float 1 0))))))) 1695 1696(defun math-ln-raw-3 (x) ; [F F] 1697 (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x)) 1698 '(float 2 0))) 1699 1700;;; Compute ln((1+x)/(1-x)) 1701(defun math-ln-raw-series (sum n x xsqr) 1702 (math-working "log" sum) 1703 (let* ((nextx (math-mul-float x xsqr)) 1704 (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) 1705 (if (math-nearly-equal-float sum nextsum) 1706 sum 1707 (math-ln-raw-series nextsum (+ n 2) nextx xsqr)))) 1708 1709(defun math-ln-plus-1-raw (x) 1710 (math-lnp1-series x 2 x (math-neg x))) 1711 1712(defun math-lnp1-series (sum n xpow x) 1713 (math-working "lnp1" sum) 1714 (let* ((nextx (math-mul-float xpow x)) 1715 (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) 1716 (if (math-nearly-equal-float sum nextsum) 1717 sum 1718 (math-lnp1-series nextsum (1+ n) nextx x)))) 1719 1720(math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21) 1721 (math-ln-raw-2 '(float 1 1))) 1722 1723(math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21) 1724 (math-ln-raw-3 (math-float '(frac 1 3)))) 1725 1726 1727 1728;;; Hyperbolic functions. 1729 1730(defun calcFunc-sinh (x) ; [N N] [Public] 1731 (cond ((eq x 0) 0) 1732 (math-expand-formulas 1733 (math-normalize 1734 (list '/ (list '- (list 'calcFunc-exp x) 1735 (list 'calcFunc-exp (list 'neg x))) 2))) 1736 ((Math-numberp x) 1737 (if calc-symbolic-mode (signal 'inexact-result nil)) 1738 (math-with-extra-prec 2 1739 (let ((expx (math-exp-raw (math-float x)))) 1740 (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1))))) 1741 ((eq (car-safe x) 'sdev) 1742 (math-make-sdev (calcFunc-sinh (nth 1 x)) 1743 (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x))))) 1744 ((eq (car x) 'intv) 1745 (math-sort-intv (nth 1 x) 1746 (calcFunc-sinh (nth 2 x)) 1747 (calcFunc-sinh (nth 3 x)))) 1748 ((or (equal x '(var inf var-inf)) 1749 (equal x '(neg (var inf var-inf))) 1750 (equal x '(var nan var-nan))) 1751 x) 1752 (t (calc-record-why 'numberp x) 1753 (list 'calcFunc-sinh x)))) 1754(put 'calcFunc-sinh 'math-expandable t) 1755 1756(defun calcFunc-cosh (x) ; [N N] [Public] 1757 (cond ((eq x 0) 1) 1758 (math-expand-formulas 1759 (math-normalize 1760 (list '/ (list '+ (list 'calcFunc-exp x) 1761 (list 'calcFunc-exp (list 'neg x))) 2))) 1762 ((Math-numberp x) 1763 (if calc-symbolic-mode (signal 'inexact-result nil)) 1764 (math-with-extra-prec 2 1765 (let ((expx (math-exp-raw (math-float x)))) 1766 (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1))))) 1767 ((eq (car-safe x) 'sdev) 1768 (math-make-sdev (calcFunc-cosh (nth 1 x)) 1769 (math-mul (nth 2 x) 1770 (calcFunc-sinh (nth 1 x))))) 1771 ((and (eq (car x) 'intv) (math-intv-constp x)) 1772 (setq x (math-abs x)) 1773 (math-sort-intv (nth 1 x) 1774 (calcFunc-cosh (nth 2 x)) 1775 (calcFunc-cosh (nth 3 x)))) 1776 ((or (equal x '(var inf var-inf)) 1777 (equal x '(neg (var inf var-inf))) 1778 (equal x '(var nan var-nan))) 1779 (math-abs x)) 1780 (t (calc-record-why 'numberp x) 1781 (list 'calcFunc-cosh x)))) 1782(put 'calcFunc-cosh 'math-expandable t) 1783 1784(defun calcFunc-tanh (x) ; [N N] [Public] 1785 (cond ((eq x 0) 0) 1786 (math-expand-formulas 1787 (math-normalize 1788 (let ((expx (list 'calcFunc-exp x)) 1789 (expmx (list 'calcFunc-exp (list 'neg x)))) 1790 (math-normalize 1791 (list '/ (list '- expx expmx) (list '+ expx expmx)))))) 1792 ((Math-numberp x) 1793 (if calc-symbolic-mode (signal 'inexact-result nil)) 1794 (math-with-extra-prec 2 1795 (let* ((expx (calcFunc-exp (math-float x))) 1796 (expmx (math-div 1 expx))) 1797 (math-div (math-sub expx expmx) 1798 (math-add expx expmx))))) 1799 ((eq (car-safe x) 'sdev) 1800 (math-make-sdev (calcFunc-tanh (nth 1 x)) 1801 (math-div (nth 2 x) 1802 (math-sqr (calcFunc-cosh (nth 1 x)))))) 1803 ((eq (car x) 'intv) 1804 (math-sort-intv (nth 1 x) 1805 (calcFunc-tanh (nth 2 x)) 1806 (calcFunc-tanh (nth 3 x)))) 1807 ((equal x '(var inf var-inf)) 1808 1) 1809 ((equal x '(neg (var inf var-inf))) 1810 -1) 1811 ((equal x '(var nan var-nan)) 1812 x) 1813 (t (calc-record-why 'numberp x) 1814 (list 'calcFunc-tanh x)))) 1815(put 'calcFunc-tanh 'math-expandable t) 1816 1817(defun calcFunc-sech (x) ; [N N] [Public] 1818 (cond ((eq x 0) 1) 1819 (math-expand-formulas 1820 (math-normalize 1821 (list '/ 2 (list '+ (list 'calcFunc-exp x) 1822 (list 'calcFunc-exp (list 'neg x)))))) 1823 ((Math-numberp x) 1824 (if calc-symbolic-mode (signal 'inexact-result nil)) 1825 (math-with-extra-prec 2 1826 (let ((expx (math-exp-raw (math-float x)))) 1827 (math-div '(float 2 0) (math-add expx (math-div 1 expx)))))) 1828 ((eq (car-safe x) 'sdev) 1829 (math-make-sdev (calcFunc-sech (nth 1 x)) 1830 (math-mul (nth 2 x) 1831 (math-mul (calcFunc-sech (nth 1 x)) 1832 (calcFunc-tanh (nth 1 x)))))) 1833 ((and (eq (car x) 'intv) (math-intv-constp x)) 1834 (setq x (math-abs x)) 1835 (math-sort-intv (nth 1 x) 1836 (calcFunc-sech (nth 2 x)) 1837 (calcFunc-sech (nth 3 x)))) 1838 ((or (equal x '(var inf var-inf)) 1839 (equal x '(neg (var inf var-inf)))) 1840 0) 1841 ((equal x '(var nan var-nan)) 1842 x) 1843 (t (calc-record-why 'numberp x) 1844 (list 'calcFunc-sech x)))) 1845(put 'calcFunc-sech 'math-expandable t) 1846 1847(defun calcFunc-csch (x) ; [N N] [Public] 1848 (cond ((eq x 0) (math-div 1 0)) 1849 (math-expand-formulas 1850 (math-normalize 1851 (list '/ 2 (list '- (list 'calcFunc-exp x) 1852 (list 'calcFunc-exp (list 'neg x)))))) 1853 ((Math-numberp x) 1854 (if calc-symbolic-mode (signal 'inexact-result nil)) 1855 (math-with-extra-prec 2 1856 (let ((expx (math-exp-raw (math-float x)))) 1857 (math-div '(float 2 0) (math-add expx (math-div -1 expx)))))) 1858 ((eq (car-safe x) 'sdev) 1859 (math-make-sdev (calcFunc-csch (nth 1 x)) 1860 (math-mul (nth 2 x) 1861 (math-mul (calcFunc-csch (nth 1 x)) 1862 (calcFunc-coth (nth 1 x)))))) 1863 ((eq (car x) 'intv) 1864 (if (and (Math-negp (nth 2 x)) 1865 (Math-posp (nth 3 x))) 1866 '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) 1867 (math-sort-intv (nth 1 x) 1868 (calcFunc-csch (nth 2 x)) 1869 (calcFunc-csch (nth 3 x))))) 1870 ((or (equal x '(var inf var-inf)) 1871 (equal x '(neg (var inf var-inf)))) 1872 0) 1873 ((equal x '(var nan var-nan)) 1874 x) 1875 (t (calc-record-why 'numberp x) 1876 (list 'calcFunc-csch x)))) 1877(put 'calcFunc-csch 'math-expandable t) 1878 1879(defun calcFunc-coth (x) ; [N N] [Public] 1880 (cond ((eq x 0) (math-div 1 0)) 1881 (math-expand-formulas 1882 (math-normalize 1883 (let ((expx (list 'calcFunc-exp x)) 1884 (expmx (list 'calcFunc-exp (list 'neg x)))) 1885 (math-normalize 1886 (list '/ (list '+ expx expmx) (list '- expx expmx)))))) 1887 ((Math-numberp x) 1888 (if calc-symbolic-mode (signal 'inexact-result nil)) 1889 (math-with-extra-prec 2 1890 (let* ((expx (calcFunc-exp (math-float x))) 1891 (expmx (math-div 1 expx))) 1892 (math-div (math-add expx expmx) 1893 (math-sub expx expmx))))) 1894 ((eq (car-safe x) 'sdev) 1895 (math-make-sdev (calcFunc-coth (nth 1 x)) 1896 (math-div (nth 2 x) 1897 (math-sqr (calcFunc-sinh (nth 1 x)))))) 1898 ((eq (car x) 'intv) 1899 (if (and (Math-negp (nth 2 x)) 1900 (Math-posp (nth 3 x))) 1901 '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) 1902 (math-sort-intv (nth 1 x) 1903 (calcFunc-coth (nth 2 x)) 1904 (calcFunc-coth (nth 3 x))))) 1905 ((equal x '(var inf var-inf)) 1906 1) 1907 ((equal x '(neg (var inf var-inf))) 1908 -1) 1909 ((equal x '(var nan var-nan)) 1910 x) 1911 (t (calc-record-why 'numberp x) 1912 (list 'calcFunc-coth x)))) 1913(put 'calcFunc-coth 'math-expandable t) 1914 1915(defun calcFunc-arcsinh (x) ; [N N] [Public] 1916 (cond ((eq x 0) 0) 1917 (math-expand-formulas 1918 (math-normalize 1919 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt 1920 (list '+ (list '^ x 2) 1)))))) 1921 ((Math-numberp x) 1922 (if calc-symbolic-mode (signal 'inexact-result nil)) 1923 (math-with-extra-prec 2 1924 (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x) 1925 '(float 1 0))))))) 1926 ((eq (car-safe x) 'sdev) 1927 (math-make-sdev (calcFunc-arcsinh (nth 1 x)) 1928 (math-div (nth 2 x) 1929 (math-sqrt 1930 (math-add (math-sqr (nth 1 x)) 1))))) 1931 ((eq (car x) 'intv) 1932 (math-sort-intv (nth 1 x) 1933 (calcFunc-arcsinh (nth 2 x)) 1934 (calcFunc-arcsinh (nth 3 x)))) 1935 ((or (equal x '(var inf var-inf)) 1936 (equal x '(neg (var inf var-inf))) 1937 (equal x '(var nan var-nan))) 1938 x) 1939 (t (calc-record-why 'numberp x) 1940 (list 'calcFunc-arcsinh x)))) 1941(put 'calcFunc-arcsinh 'math-expandable t) 1942 1943(defun calcFunc-arccosh (x) ; [N N] [Public] 1944 (cond ((eq x 1) 0) 1945 ((and (eq x -1) calc-symbolic-mode) 1946 '(var pi var-pi)) 1947 ((and (eq x 0) calc-symbolic-mode) 1948 (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2)) 1949 (math-expand-formulas 1950 (math-normalize 1951 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt 1952 (list '- (list '^ x 2) 1)))))) 1953 ((Math-numberp x) 1954 (if calc-symbolic-mode (signal 'inexact-result nil)) 1955 (if (Math-equal-int x -1) 1956 (math-imaginary (math-pi)) 1957 (math-with-extra-prec 2 1958 (if (or t ; need to do this even in the real case! 1959 (memq (car-safe x) '(cplx polar))) 1960 (let ((xp1 (math-add 1 x))) ; this gets the branch cuts right 1961 (math-ln-raw 1962 (math-add x (math-mul xp1 1963 (math-sqrt-raw 1964 (math-div (math-sub 1965 x 1966 '(float 1 0)) 1967 xp1)))))) 1968 (math-ln-raw 1969 (math-add x (math-sqrt-raw (math-add (math-sqr x) 1970 '(float -1 0))))))))) 1971 ((eq (car-safe x) 'sdev) 1972 (math-make-sdev (calcFunc-arccosh (nth 1 x)) 1973 (math-div (nth 2 x) 1974 (math-sqrt 1975 (math-add (math-sqr (nth 1 x)) -1))))) 1976 ((eq (car x) 'intv) 1977 (math-sort-intv (nth 1 x) 1978 (calcFunc-arccosh (nth 2 x)) 1979 (calcFunc-arccosh (nth 3 x)))) 1980 ((or (equal x '(var inf var-inf)) 1981 (equal x '(neg (var inf var-inf))) 1982 (equal x '(var nan var-nan))) 1983 x) 1984 (t (calc-record-why 'numberp x) 1985 (list 'calcFunc-arccosh x)))) 1986(put 'calcFunc-arccosh 'math-expandable t) 1987 1988(defun calcFunc-arctanh (x) ; [N N] [Public] 1989 (cond ((eq x 0) 0) 1990 ((and (Math-equal-int x 1) calc-infinite-mode) 1991 '(var inf var-inf)) 1992 ((and (Math-equal-int x -1) calc-infinite-mode) 1993 '(neg (var inf var-inf))) 1994 (math-expand-formulas 1995 (list '/ (list '- 1996 (list 'calcFunc-ln (list '+ 1 x)) 1997 (list 'calcFunc-ln (list '- 1 x))) 2)) 1998 ((Math-numberp x) 1999 (if calc-symbolic-mode (signal 'inexact-result nil)) 2000 (math-with-extra-prec 2 2001 (if (or (memq (car-safe x) '(cplx polar)) 2002 (Math-lessp 1 x)) 2003 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x)) 2004 (math-ln-raw (math-sub '(float 1 0) x))) 2005 '(float 5 -1)) 2006 (if (and (math-equal-int x 1) calc-infinite-mode) 2007 '(var inf var-inf) 2008 (if (and (math-equal-int x -1) calc-infinite-mode) 2009 '(neg (var inf var-inf)) 2010 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x) 2011 (math-sub 1 x))) 2012 '(float 5 -1))))))) 2013 ((eq (car-safe x) 'sdev) 2014 (math-make-sdev (calcFunc-arctanh (nth 1 x)) 2015 (math-div (nth 2 x) 2016 (math-sub 1 (math-sqr (nth 1 x)))))) 2017 ((eq (car x) 'intv) 2018 (math-sort-intv (nth 1 x) 2019 (calcFunc-arctanh (nth 2 x)) 2020 (calcFunc-arctanh (nth 3 x)))) 2021 ((equal x '(var nan var-nan)) 2022 x) 2023 (t (calc-record-why 'numberp x) 2024 (list 'calcFunc-arctanh x)))) 2025(put 'calcFunc-arctanh 'math-expandable t) 2026 2027 2028;;; Convert A from HMS or degrees to radians. 2029(defun calcFunc-rad (a) ; [R R] [Public] 2030 (cond ((or (Math-numberp a) 2031 (eq (car a) 'intv)) 2032 (math-with-extra-prec 2 2033 (math-mul a (math-pi-over-180)))) 2034 ((eq (car a) 'hms) 2035 (math-from-hms a 'rad)) 2036 ((eq (car a) 'sdev) 2037 (math-make-sdev (calcFunc-rad (nth 1 a)) 2038 (calcFunc-rad (nth 2 a)))) 2039 (math-expand-formulas 2040 (math-div (math-mul a '(var pi var-pi)) 180)) 2041 ((math-infinitep a) a) 2042 (t (list 'calcFunc-rad a)))) 2043(put 'calcFunc-rad 'math-expandable t) 2044 2045;;; Convert A from HMS or radians to degrees. 2046(defun calcFunc-deg (a) ; [R R] [Public] 2047 (cond ((or (Math-numberp a) 2048 (eq (car a) 'intv)) 2049 (math-with-extra-prec 2 2050 (math-div a (math-pi-over-180)))) 2051 ((eq (car a) 'hms) 2052 (math-from-hms a 'deg)) 2053 ((eq (car a) 'sdev) 2054 (math-make-sdev (calcFunc-deg (nth 1 a)) 2055 (calcFunc-deg (nth 2 a)))) 2056 (math-expand-formulas 2057 (math-div (math-mul 180 a) '(var pi var-pi))) 2058 ((math-infinitep a) a) 2059 (t (list 'calcFunc-deg a)))) 2060(put 'calcFunc-deg 'math-expandable t) 2061 2062(provide 'calc-math) 2063 2064;;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c 2065;;; calc-math.el ends here 2066