1;;; calc-funcs.el --- well-known 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-inc-gamma (arg) 36 (interactive "P") 37 (calc-slow-wrapper 38 (if (calc-is-inverse) 39 (if (calc-is-hyperbolic) 40 (calc-binary-op "gamG" 'calcFunc-gammaG arg) 41 (calc-binary-op "gamQ" 'calcFunc-gammaQ arg)) 42 (if (calc-is-hyperbolic) 43 (calc-binary-op "gamg" 'calcFunc-gammag arg) 44 (calc-binary-op "gamP" 'calcFunc-gammaP arg))))) 45 46(defun calc-erf (arg) 47 (interactive "P") 48 (calc-slow-wrapper 49 (if (calc-is-inverse) 50 (calc-unary-op "erfc" 'calcFunc-erfc arg) 51 (calc-unary-op "erf" 'calcFunc-erf arg)))) 52 53(defun calc-erfc (arg) 54 (interactive "P") 55 (calc-invert-func) 56 (calc-erf arg)) 57 58(defun calc-beta (arg) 59 (interactive "P") 60 (calc-slow-wrapper 61 (calc-binary-op "beta" 'calcFunc-beta arg))) 62 63(defun calc-inc-beta () 64 (interactive) 65 (calc-slow-wrapper 66 (if (calc-is-hyperbolic) 67 (calc-enter-result 3 "betB" (cons 'calcFunc-betaB (calc-top-list-n 3))) 68 (calc-enter-result 3 "betI" (cons 'calcFunc-betaI (calc-top-list-n 3)))))) 69 70(defun calc-bessel-J (arg) 71 (interactive "P") 72 (calc-slow-wrapper 73 (calc-binary-op "besJ" 'calcFunc-besJ arg))) 74 75(defun calc-bessel-Y (arg) 76 (interactive "P") 77 (calc-slow-wrapper 78 (calc-binary-op "besY" 'calcFunc-besY arg))) 79 80(defun calc-bernoulli-number (arg) 81 (interactive "P") 82 (calc-slow-wrapper 83 (if (calc-is-hyperbolic) 84 (calc-binary-op "bern" 'calcFunc-bern arg) 85 (calc-unary-op "bern" 'calcFunc-bern arg)))) 86 87(defun calc-euler-number (arg) 88 (interactive "P") 89 (calc-slow-wrapper 90 (if (calc-is-hyperbolic) 91 (calc-binary-op "eulr" 'calcFunc-euler arg) 92 (calc-unary-op "eulr" 'calcFunc-euler arg)))) 93 94(defun calc-stirling-number (arg) 95 (interactive "P") 96 (calc-slow-wrapper 97 (if (calc-is-hyperbolic) 98 (calc-binary-op "str2" 'calcFunc-stir2 arg) 99 (calc-binary-op "str1" 'calcFunc-stir1 arg)))) 100 101(defun calc-utpb () 102 (interactive) 103 (calc-prob-dist "b" 3)) 104 105(defun calc-utpc () 106 (interactive) 107 (calc-prob-dist "c" 2)) 108 109(defun calc-utpf () 110 (interactive) 111 (calc-prob-dist "f" 3)) 112 113(defun calc-utpn () 114 (interactive) 115 (calc-prob-dist "n" 3)) 116 117(defun calc-utpp () 118 (interactive) 119 (calc-prob-dist "p" 2)) 120 121(defun calc-utpt () 122 (interactive) 123 (calc-prob-dist "t" 2)) 124 125(defun calc-prob-dist (letter nargs) 126 (calc-slow-wrapper 127 (if (calc-is-inverse) 128 (calc-enter-result nargs (concat "ltp" letter) 129 (append (list (intern (concat "calcFunc-ltp" letter)) 130 (calc-top-n 1)) 131 (calc-top-list-n (1- nargs) 2))) 132 (calc-enter-result nargs (concat "utp" letter) 133 (append (list (intern (concat "calcFunc-utp" letter)) 134 (calc-top-n 1)) 135 (calc-top-list-n (1- nargs) 2)))))) 136 137 138 139 140;;; Sources: Numerical Recipes, Press et al; 141;;; Handbook of Mathematical Functions, Abramowitz & Stegun. 142 143 144;;; Gamma function. 145 146(defun calcFunc-gamma (x) 147 (or (math-numberp x) (math-reject-arg x 'numberp)) 148 (calcFunc-fact (math-add x -1))) 149 150(defun math-gammap1-raw (x &optional fprec nfprec) ; compute gamma(1 + x) 151 (or fprec 152 (setq fprec (math-float calc-internal-prec) 153 nfprec (math-float (- calc-internal-prec)))) 154 (cond ((math-lessp-float (calcFunc-re x) fprec) 155 (if (math-lessp-float (calcFunc-re x) nfprec) 156 (math-neg (math-div 157 (math-pi) 158 (math-mul (math-gammap1-raw 159 (math-add (math-neg x) 160 '(float -1 0)) 161 fprec nfprec) 162 (math-sin-raw 163 (math-mul (math-pi) x))))) 164 (let ((xplus1 (math-add x '(float 1 0)))) 165 (math-div (math-gammap1-raw xplus1 fprec nfprec) xplus1)))) 166 ((and (math-realp x) 167 (math-lessp-float '(float 736276 0) x)) 168 (math-overflow)) 169 (t ; re(x) now >= 10.0 170 (let ((xinv (math-div 1 x)) 171 (lnx (math-ln-raw x))) 172 (math-mul (math-sqrt-two-pi) 173 (math-exp-raw 174 (math-gamma-series 175 (math-sub (math-mul (math-add x '(float 5 -1)) 176 lnx) 177 x) 178 xinv 179 (math-sqr xinv) 180 '(float 0 0) 181 2))))))) 182 183(defun math-gamma-series (sum x xinvsqr oterm n) 184 (math-working "gamma" sum) 185 (let* ((bn (math-bernoulli-number n)) 186 (term (math-mul (math-div-float (math-float (nth 1 bn)) 187 (math-float (* (nth 2 bn) 188 (* n (1- n))))) 189 x)) 190 (next (math-add sum term))) 191 (if (math-nearly-equal sum next) 192 next 193 (if (> n (* 2 calc-internal-prec)) 194 (progn 195 ;; Need this because series eventually diverges for large enough n. 196 (calc-record-why 197 "*Gamma computation stopped early, not all digits may be valid") 198 next) 199 (math-gamma-series next (math-mul x xinvsqr) xinvsqr term (+ n 2)))))) 200 201 202;;; Incomplete gamma function. 203 204(defvar math-current-gamma-value nil) 205(defun calcFunc-gammaP (a x) 206 (if (equal x '(var inf var-inf)) 207 '(float 1 0) 208 (math-inexact-result) 209 (or (Math-numberp a) (math-reject-arg a 'numberp)) 210 (or (math-numberp x) (math-reject-arg x 'numberp)) 211 (if (and (math-num-integerp a) 212 (integerp (setq a (math-trunc a))) 213 (> a 0) (< a 20)) 214 (math-sub 1 (calcFunc-gammaQ a x)) 215 (let ((math-current-gamma-value (calcFunc-gamma a))) 216 (math-div (calcFunc-gammag a x) math-current-gamma-value))))) 217 218(defun calcFunc-gammaQ (a x) 219 (if (equal x '(var inf var-inf)) 220 '(float 0 0) 221 (math-inexact-result) 222 (or (Math-numberp a) (math-reject-arg a 'numberp)) 223 (or (math-numberp x) (math-reject-arg x 'numberp)) 224 (if (and (math-num-integerp a) 225 (integerp (setq a (math-trunc a))) 226 (> a 0) (< a 20)) 227 (let ((n 0) 228 (sum '(float 1 0)) 229 (term '(float 1 0))) 230 (math-with-extra-prec 1 231 (while (< (setq n (1+ n)) a) 232 (setq term (math-div (math-mul term x) n) 233 sum (math-add sum term)) 234 (math-working "gamma" sum)) 235 (math-mul sum (calcFunc-exp (math-neg x))))) 236 (let ((math-current-gamma-value (calcFunc-gamma a))) 237 (math-div (calcFunc-gammaG a x) math-current-gamma-value))))) 238 239(defun calcFunc-gammag (a x) 240 (if (equal x '(var inf var-inf)) 241 (calcFunc-gamma a) 242 (math-inexact-result) 243 (or (Math-numberp a) (math-reject-arg a 'numberp)) 244 (or (Math-numberp x) (math-reject-arg x 'numberp)) 245 (math-with-extra-prec 2 246 (setq a (math-float a)) 247 (setq x (math-float x)) 248 (if (or (math-negp (calcFunc-re a)) 249 (math-lessp-float (calcFunc-re x) 250 (math-add-float (calcFunc-re a) 251 '(float 1 0)))) 252 (math-inc-gamma-series a x) 253 (math-sub (or math-current-gamma-value (calcFunc-gamma a)) 254 (math-inc-gamma-cfrac a x)))))) 255 256(defun calcFunc-gammaG (a x) 257 (if (equal x '(var inf var-inf)) 258 '(float 0 0) 259 (math-inexact-result) 260 (or (Math-numberp a) (math-reject-arg a 'numberp)) 261 (or (Math-numberp x) (math-reject-arg x 'numberp)) 262 (math-with-extra-prec 2 263 (setq a (math-float a)) 264 (setq x (math-float x)) 265 (if (or (math-negp (calcFunc-re a)) 266 (math-lessp-float (calcFunc-re x) 267 (math-add-float (math-abs-approx a) 268 '(float 1 0)))) 269 (math-sub (or math-current-gamma-value (calcFunc-gamma a)) 270 (math-inc-gamma-series a x)) 271 (math-inc-gamma-cfrac a x))))) 272 273(defun math-inc-gamma-series (a x) 274 (if (Math-zerop x) 275 '(float 0 0) 276 (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x)) 277 (math-with-extra-prec 2 278 (let ((start (math-div '(float 1 0) a))) 279 (math-inc-gamma-series-step start start a x)))))) 280 281(defun math-inc-gamma-series-step (sum term a x) 282 (math-working "gamma" sum) 283 (setq a (math-add a '(float 1 0)) 284 term (math-div (math-mul term x) a)) 285 (let ((next (math-add sum term))) 286 (if (math-nearly-equal sum next) 287 next 288 (math-inc-gamma-series-step next term a x)))) 289 290(defun math-inc-gamma-cfrac (a x) 291 (if (Math-zerop x) 292 (or math-current-gamma-value (calcFunc-gamma a)) 293 (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x)) 294 (math-inc-gamma-cfrac-step '(float 1 0) x 295 '(float 0 0) '(float 1 0) 296 '(float 1 0) '(float 1 0) '(float 0 0) 297 a x)))) 298 299(defun math-inc-gamma-cfrac-step (a0 a1 b0 b1 n fac g a x) 300 (let ((ana (math-sub n a)) 301 (anf (math-mul n fac))) 302 (setq n (math-add n '(float 1 0)) 303 a0 (math-mul (math-add a1 (math-mul a0 ana)) fac) 304 b0 (math-mul (math-add b1 (math-mul b0 ana)) fac) 305 a1 (math-add (math-mul x a0) (math-mul anf a1)) 306 b1 (math-add (math-mul x b0) (math-mul anf b1))) 307 (if (math-zerop a1) 308 (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac g a x) 309 (setq fac (math-div '(float 1 0) a1)) 310 (let ((next (math-mul b1 fac))) 311 (math-working "gamma" next) 312 (if (math-nearly-equal next g) 313 next 314 (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac next a x)))))) 315 316 317;;; Error function. 318 319(defun calcFunc-erf (x) 320 (if (equal x '(var inf var-inf)) 321 '(float 1 0) 322 (if (equal x '(neg (var inf var-inf))) 323 '(float -1 0) 324 (if (Math-zerop x) 325 x 326 (let ((math-current-gamma-value (math-sqrt-pi))) 327 (math-to-same-complex-quad 328 (math-div (calcFunc-gammag '(float 5 -1) 329 (math-sqr (math-to-complex-quad-one x))) 330 math-current-gamma-value) 331 x)))))) 332 333(defun calcFunc-erfc (x) 334 (if (equal x '(var inf var-inf)) 335 '(float 0 0) 336 (if (math-posp x) 337 (let ((math-current-gamma-value (math-sqrt-pi))) 338 (math-div (calcFunc-gammaG '(float 5 -1) (math-sqr x)) 339 math-current-gamma-value)) 340 (math-sub 1 (calcFunc-erf x))))) 341 342(defun math-to-complex-quad-one (x) 343 (if (eq (car-safe x) 'polar) (setq x (math-complex x))) 344 (if (eq (car-safe x) 'cplx) 345 (list 'cplx (math-abs (nth 1 x)) (math-abs (nth 2 x))) 346 x)) 347 348(defun math-to-same-complex-quad (x y) 349 (if (eq (car-safe y) 'cplx) 350 (if (eq (car-safe x) 'cplx) 351 (list 'cplx 352 (if (math-negp (nth 1 y)) (math-neg (nth 1 x)) (nth 1 x)) 353 (if (math-negp (nth 2 y)) (math-neg (nth 2 x)) (nth 2 x))) 354 (if (math-negp (nth 1 y)) (math-neg x) x)) 355 (if (math-negp y) 356 (if (eq (car-safe x) 'cplx) 357 (list 'cplx (math-neg (nth 1 x)) (nth 2 x)) 358 (math-neg x)) 359 x))) 360 361 362;;; Beta function. 363 364(defun calcFunc-beta (a b) 365 (if (math-num-integerp a) 366 (let ((am (math-add a -1))) 367 (or (math-numberp b) (math-reject-arg b 'numberp)) 368 (math-div 1 (math-mul b (calcFunc-choose (math-add b am) am)))) 369 (if (math-num-integerp b) 370 (calcFunc-beta b a) 371 (math-div (math-mul (calcFunc-gamma a) (calcFunc-gamma b)) 372 (calcFunc-gamma (math-add a b)))))) 373 374 375;;; Incomplete beta function. 376 377(defvar math-current-beta-value nil) 378(defun calcFunc-betaI (x a b) 379 (cond ((math-zerop x) 380 '(float 0 0)) 381 ((math-equal-int x 1) 382 '(float 1 0)) 383 ((or (math-zerop a) 384 (and (math-num-integerp a) 385 (math-negp a))) 386 (if (or (math-zerop b) 387 (and (math-num-integerp b) 388 (math-negp b))) 389 (math-reject-arg b 'range) 390 '(float 1 0))) 391 ((or (math-zerop b) 392 (and (math-num-integerp b) 393 (math-negp b))) 394 '(float 0 0)) 395 ((not (math-numberp a)) (math-reject-arg a 'numberp)) 396 ((not (math-numberp b)) (math-reject-arg b 'numberp)) 397 ((math-inexact-result)) 398 (t (let ((math-current-beta-value (calcFunc-beta a b))) 399 (math-div (calcFunc-betaB x a b) math-current-beta-value))))) 400 401(defun calcFunc-betaB (x a b) 402 (cond 403 ((math-zerop x) 404 '(float 0 0)) 405 ((math-equal-int x 1) 406 (calcFunc-beta a b)) 407 ((not (math-numberp x)) (math-reject-arg x 'numberp)) 408 ((not (math-numberp a)) (math-reject-arg a 'numberp)) 409 ((not (math-numberp b)) (math-reject-arg b 'numberp)) 410 ((math-zerop a) (math-reject-arg a 'nonzerop)) 411 ((math-zerop b) (math-reject-arg b 'nonzerop)) 412 ((and (math-num-integerp b) 413 (if (math-negp b) 414 (math-reject-arg b 'range) 415 (Math-natnum-lessp (setq b (math-trunc b)) 20))) 416 (and calc-symbolic-mode (or (math-floatp a) (math-floatp b)) 417 (math-inexact-result)) 418 (math-mul 419 (math-with-extra-prec 2 420 (let* ((i 0) 421 (term 1) 422 (sum (math-div term a))) 423 (while (< (setq i (1+ i)) b) 424 (setq term (math-mul (math-div (math-mul term (- i b)) i) x) 425 sum (math-add sum (math-div term (math-add a i)))) 426 (math-working "beta" sum)) 427 sum)) 428 (math-pow x a))) 429 ((and (math-num-integerp a) 430 (if (math-negp a) 431 (math-reject-arg a 'range) 432 (Math-natnum-lessp (setq a (math-trunc a)) 20))) 433 (math-sub (or math-current-beta-value (calcFunc-beta a b)) 434 (calcFunc-betaB (math-sub 1 x) b a))) 435 (t 436 (math-inexact-result) 437 (math-with-extra-prec 2 438 (setq x (math-float x)) 439 (setq a (math-float a)) 440 (setq b (math-float b)) 441 (let ((bt (math-exp-raw (math-add (math-mul a (math-ln-raw x)) 442 (math-mul b (math-ln-raw 443 (math-sub '(float 1 0) 444 x))))))) 445 (if (Math-lessp x (math-div (math-add a '(float 1 0)) 446 (math-add (math-add a b) '(float 2 0)))) 447 (math-div (math-mul bt (math-beta-cfrac a b x)) a) 448 (math-sub (or math-current-beta-value (calcFunc-beta a b)) 449 (math-div (math-mul bt 450 (math-beta-cfrac b a (math-sub 1 x))) 451 b)))))))) 452 453(defun math-beta-cfrac (a b x) 454 (let ((qab (math-add a b)) 455 (qap (math-add a '(float 1 0))) 456 (qam (math-add a '(float -1 0)))) 457 (math-beta-cfrac-step '(float 1 0) 458 (math-sub '(float 1 0) 459 (math-div (math-mul qab x) qap)) 460 '(float 1 0) '(float 1 0) 461 '(float 1 0) 462 qab qap qam a b x))) 463 464(defun math-beta-cfrac-step (az bz am bm m qab qap qam a b x) 465 (let* ((two-m (math-mul m '(float 2 0))) 466 (d (math-div (math-mul (math-mul (math-sub b m) m) x) 467 (math-mul (math-add qam two-m) (math-add a two-m)))) 468 (ap (math-add az (math-mul d am))) 469 (bp (math-add bz (math-mul d bm))) 470 (d2 (math-neg 471 (math-div (math-mul (math-mul (math-add a m) (math-add qab m)) x) 472 (math-mul (math-add qap two-m) (math-add a two-m))))) 473 (app (math-add ap (math-mul d2 az))) 474 (bpp (math-add bp (math-mul d2 bz))) 475 (next (math-div app bpp))) 476 (math-working "beta" next) 477 (if (math-nearly-equal next az) 478 next 479 (math-beta-cfrac-step next '(float 1 0) 480 (math-div ap bpp) (math-div bp bpp) 481 (math-add m '(float 1 0)) 482 qab qap qam a b x)))) 483 484 485;;; Bessel functions. 486 487;;; Should generalize this to handle arbitrary precision! 488 489(defun calcFunc-besJ (v x) 490 (or (math-numberp v) (math-reject-arg v 'numberp)) 491 (or (math-numberp x) (math-reject-arg x 'numberp)) 492 (let ((calc-internal-prec (min 8 calc-internal-prec))) 493 (math-with-extra-prec 3 494 (setq x (math-float (math-normalize x))) 495 (setq v (math-float (math-normalize v))) 496 (cond ((math-zerop x) 497 (if (math-zerop v) 498 '(float 1 0) 499 '(float 0 0))) 500 ((math-inexact-result)) 501 ((not (math-num-integerp v)) 502 (let ((start (math-div 1 (calcFunc-fact v)))) 503 (math-mul (math-besJ-series start start 504 0 505 (math-mul '(float -25 -2) 506 (math-sqr x)) 507 v) 508 (math-pow (math-div x 2) v)))) 509 ((math-negp (setq v (math-trunc v))) 510 (if (math-oddp v) 511 (math-neg (calcFunc-besJ (math-neg v) x)) 512 (calcFunc-besJ (math-neg v) x))) 513 ((eq v 0) 514 (math-besJ0 x)) 515 ((eq v 1) 516 (math-besJ1 x)) 517 ((Math-lessp v (math-abs-approx x)) 518 (let ((j 0) 519 (bjm (math-besJ0 x)) 520 (bj (math-besJ1 x)) 521 (two-over-x (math-div 2 x)) 522 bjp) 523 (while (< (setq j (1+ j)) v) 524 (setq bjp (math-sub (math-mul (math-mul j two-over-x) bj) 525 bjm) 526 bjm bj 527 bj bjp)) 528 bj)) 529 (t 530 (if (Math-lessp 100 v) (math-reject-arg v 'range)) 531 (let* ((j (logior (+ v (math-isqrt-small (* 40 v))) 1)) 532 (two-over-x (math-div 2 x)) 533 (jsum nil) 534 (bjp '(float 0 0)) 535 (sum '(float 0 0)) 536 (bj '(float 1 0)) 537 bjm ans) 538 (while (> (setq j (1- j)) 0) 539 (setq bjm (math-sub (math-mul (math-mul j two-over-x) bj) 540 bjp) 541 bjp bj 542 bj bjm) 543 (if (> (nth 2 (math-abs-approx bj)) 10) 544 (setq bj (math-mul bj '(float 1 -10)) 545 bjp (math-mul bjp '(float 1 -10)) 546 ans (and ans (math-mul ans '(float 1 -10))) 547 sum (math-mul sum '(float 1 -10)))) 548 (or (setq jsum (not jsum)) 549 (setq sum (math-add sum bj))) 550 (if (= j v) 551 (setq ans bjp))) 552 (math-div ans (math-sub (math-mul 2 sum) bj)))))))) 553 554(defun math-besJ-series (sum term k zz vk) 555 (math-working "besJ" sum) 556 (setq k (1+ k) 557 vk (math-add 1 vk) 558 term (math-div (math-mul term zz) (math-mul k vk))) 559 (let ((next (math-add sum term))) 560 (if (math-nearly-equal next sum) 561 next 562 (math-besJ-series next term k zz vk)))) 563 564(defun math-besJ0 (x &optional yflag) 565 (cond ((and (not yflag) (math-negp (calcFunc-re x))) 566 (math-besJ0 (math-neg x))) 567 ((Math-lessp '(float 8 0) (math-abs-approx x)) 568 (let* ((z (math-div '(float 8 0) x)) 569 (y (math-sqr z)) 570 (xx (math-add x '(float (bigneg 164 398 785) -9))) 571 (a1 (math-poly-eval y 572 '((float (bigpos 211 887 093 2) -16) 573 (float (bigneg 639 370 073 2) -15) 574 (float (bigpos 407 510 734 2) -14) 575 (float (bigneg 627 628 098 1) -12) 576 (float 1 0)))) 577 (a2 (math-poly-eval y 578 '((float (bigneg 152 935 934) -16) 579 (float (bigpos 161 095 621 7) -16) 580 (float (bigneg 651 147 911 6) -15) 581 (float (bigpos 765 488 430 1) -13) 582 (float (bigneg 995 499 562 1) -11)))) 583 (sc (math-sin-cos-raw xx))) 584 (if yflag 585 (setq sc (cons (math-neg (cdr sc)) (car sc)))) 586 (math-mul (math-sqrt 587 (math-div '(float (bigpos 722 619 636) -9) x)) 588 (math-sub (math-mul (cdr sc) a1) 589 (math-mul (car sc) (math-mul z a2)))))) 590 (t 591 (let ((y (math-sqr x))) 592 (math-div (math-poly-eval y 593 '((float (bigneg 456 052 849 1) -7) 594 (float (bigpos 017 233 739 7) -5) 595 (float (bigneg 418 442 121 1) -2) 596 (float (bigpos 407 196 516 6) -1) 597 (float (bigneg 354 590 362 13) 0) 598 (float (bigpos 574 490 568 57) 0))) 599 (math-poly-eval y 600 '((float 1 0) 601 (float (bigpos 712 532 678 2) -7) 602 (float (bigpos 853 264 927 5) -5) 603 (float (bigpos 718 680 494 9) -3) 604 (float (bigpos 985 532 029 1) 0) 605 (float (bigpos 411 490 568 57) 0)))))))) 606 607(defun math-besJ1 (x &optional yflag) 608 (cond ((and (math-negp (calcFunc-re x)) (not yflag)) 609 (math-neg (math-besJ1 (math-neg x)))) 610 ((Math-lessp '(float 8 0) (math-abs-approx x)) 611 (let* ((z (math-div '(float 8 0) x)) 612 (y (math-sqr z)) 613 (xx (math-add x '(float (bigneg 491 194 356 2) -9))) 614 (a1 (math-poly-eval y 615 '((float (bigneg 019 337 240) -15) 616 (float (bigpos 174 520 457 2) -15) 617 (float (bigneg 496 396 516 3) -14) 618 (float 183105 -8) 619 (float 1 0)))) 620 (a2 (math-poly-eval y 621 '((float (bigpos 412 787 105) -15) 622 (float (bigneg 987 228 88) -14) 623 (float (bigpos 096 199 449 8) -15) 624 (float (bigneg 873 690 002 2) -13) 625 (float (bigpos 995 499 687 4) -11)))) 626 (sc (math-sin-cos-raw xx))) 627 (if yflag 628 (setq sc (cons (math-neg (cdr sc)) (car sc))) 629 (if (math-negp x) 630 (setq sc (cons (math-neg (car sc)) (math-neg (cdr sc)))))) 631 (math-mul (math-sqrt (math-div '(float (bigpos 722 619 636) -9) x)) 632 (math-sub (math-mul (cdr sc) a1) 633 (math-mul (car sc) (math-mul z a2)))))) 634 (t 635 (let ((y (math-sqr x))) 636 (math-mul 637 x 638 (math-div (math-poly-eval y 639 '((float (bigneg 606 036 016 3) -8) 640 (float (bigpos 826 044 157) -4) 641 (float (bigneg 439 611 972 2) -3) 642 (float (bigpos 531 968 423 2) -1) 643 (float (bigneg 235 059 895 7) 0) 644 (float (bigpos 232 614 362 72) 0))) 645 (math-poly-eval y 646 '((float 1 0) 647 (float (bigpos 397 991 769 3) -7) 648 (float (bigpos 394 743 944 9) -5) 649 (float (bigpos 474 330 858 1) -2) 650 (float (bigpos 178 535 300 2) 0) 651 (float (bigpos 442 228 725 144) 652 0))))))))) 653 654(defun calcFunc-besY (v x) 655 (math-inexact-result) 656 (or (math-numberp v) (math-reject-arg v 'numberp)) 657 (or (math-numberp x) (math-reject-arg x 'numberp)) 658 (let ((calc-internal-prec (min 8 calc-internal-prec))) 659 (math-with-extra-prec 3 660 (setq x (math-float (math-normalize x))) 661 (setq v (math-float (math-normalize v))) 662 (cond ((not (math-num-integerp v)) 663 (let ((sc (math-sin-cos-raw (math-mul v (math-pi))))) 664 (math-div (math-sub (math-mul (calcFunc-besJ v x) (cdr sc)) 665 (calcFunc-besJ (math-neg v) x)) 666 (car sc)))) 667 ((math-negp (setq v (math-trunc v))) 668 (if (math-oddp v) 669 (math-neg (calcFunc-besY (math-neg v) x)) 670 (calcFunc-besY (math-neg v) x))) 671 ((eq v 0) 672 (math-besY0 x)) 673 ((eq v 1) 674 (math-besY1 x)) 675 (t 676 (let ((j 0) 677 (bym (math-besY0 x)) 678 (by (math-besY1 x)) 679 (two-over-x (math-div 2 x)) 680 byp) 681 (while (< (setq j (1+ j)) v) 682 (setq byp (math-sub (math-mul (math-mul j two-over-x) by) 683 bym) 684 bym by 685 by byp)) 686 by)))))) 687 688(defun math-besY0 (x) 689 (cond ((Math-lessp (math-abs-approx x) '(float 8 0)) 690 (let ((y (math-sqr x))) 691 (math-add 692 (math-div (math-poly-eval y 693 '((float (bigpos 733 622 284 2) -7) 694 (float (bigneg 757 792 632 8) -5) 695 (float (bigpos 129 988 087 1) -2) 696 (float (bigneg 036 598 123 5) -1) 697 (float (bigpos 065 834 062 7) 0) 698 (float (bigneg 389 821 957 2) 0))) 699 (math-poly-eval y 700 '((float 1 0) 701 (float (bigpos 244 030 261 2) -7) 702 (float (bigpos 647 472 474) -4) 703 (float (bigpos 438 466 189 7) -3) 704 (float (bigpos 648 499 452 7) -1) 705 (float (bigpos 269 544 076 40) 0)))) 706 (math-mul '(float (bigpos 772 619 636) -9) 707 (math-mul (math-besJ0 x) (math-ln-raw x)))))) 708 ((math-negp (calcFunc-re x)) 709 (math-add (math-besJ0 (math-neg x) t) 710 (math-mul '(cplx 0 2) 711 (math-besJ0 (math-neg x))))) 712 (t 713 (math-besJ0 x t)))) 714 715(defun math-besY1 (x) 716 (cond ((Math-lessp (math-abs-approx x) '(float 8 0)) 717 (let ((y (math-sqr x))) 718 (math-add 719 (math-mul 720 x 721 (math-div (math-poly-eval y 722 '((float (bigpos 935 937 511 8) -6) 723 (float (bigneg 726 922 237 4) -3) 724 (float (bigpos 551 264 349 7) -1) 725 (float (bigneg 139 438 153 5) 1) 726 (float (bigpos 439 527 127) 4) 727 (float (bigneg 943 604 900 4) 3))) 728 (math-poly-eval y 729 '((float 1 0) 730 (float (bigpos 885 632 549 3) -7) 731 (float (bigpos 605 042 102) -3) 732 (float (bigpos 002 904 245 2) -2) 733 (float (bigpos 367 650 733 3) 0) 734 (float (bigpos 664 419 244 4) 2) 735 (float (bigpos 057 958 249) 5))))) 736 (math-mul '(float (bigpos 772 619 636) -9) 737 (math-sub (math-mul (math-besJ1 x) (math-ln-raw x)) 738 (math-div 1 x)))))) 739 ((math-negp (calcFunc-re x)) 740 (math-neg 741 (math-add (math-besJ1 (math-neg x) t) 742 (math-mul '(cplx 0 2) 743 (math-besJ1 (math-neg x)))))) 744 (t 745 (math-besJ1 x t)))) 746 747(defun math-poly-eval (x coefs) 748 (let ((accum (car coefs))) 749 (while (setq coefs (cdr coefs)) 750 (setq accum (math-add (car coefs) (math-mul accum x)))) 751 accum)) 752 753 754;;;; Bernoulli and Euler polynomials and numbers. 755 756(defun calcFunc-bern (n &optional x) 757 (if (and x (not (math-zerop x))) 758 (if (and calc-symbolic-mode (math-floatp x)) 759 (math-inexact-result) 760 (math-build-polynomial-expr (math-bernoulli-coefs n) x)) 761 (or (math-num-natnump n) (math-reject-arg n 'natnump)) 762 (if (consp n) 763 (progn 764 (math-inexact-result) 765 (math-float (math-bernoulli-number (math-trunc n)))) 766 (math-bernoulli-number n)))) 767 768(defun calcFunc-euler (n &optional x) 769 (or (math-num-natnump n) (math-reject-arg n 'natnump)) 770 (if x 771 (let* ((n1 (math-add n 1)) 772 (coefs (math-bernoulli-coefs n1)) 773 (fac (math-div (math-pow 2 n1) n1)) 774 (k -1) 775 (x1 (math-div (math-add x 1) 2)) 776 (x2 (math-div x 2))) 777 (if (math-numberp x) 778 (if (and calc-symbolic-mode (math-floatp x)) 779 (math-inexact-result) 780 (math-mul fac 781 (math-sub (math-build-polynomial-expr coefs x1) 782 (math-build-polynomial-expr coefs x2)))) 783 (calcFunc-collect 784 (math-reduce-vec 785 'math-add 786 (cons 'vec 787 (mapcar (function 788 (lambda (c) 789 (setq k (1+ k)) 790 (math-mul (math-mul fac c) 791 (math-sub (math-pow x1 k) 792 (math-pow x2 k))))) 793 coefs))) 794 x))) 795 (math-mul (math-pow 2 n) 796 (if (consp n) 797 (progn 798 (math-inexact-result) 799 (calcFunc-euler n '(float 5 -1))) 800 (calcFunc-euler n '(frac 1 2)))))) 801 802(defvar math-bernoulli-b-cache '((frac -174611 803 (bigpos 0 200 291 698 662 857 802)) 804 (frac 43867 (bigpos 0 944 170 217 94 109 5)) 805 (frac -3617 (bigpos 0 880 842 622 670 10)) 806 (frac 1 (bigpos 600 249 724 74)) 807 (frac -691 (bigpos 0 368 674 307 1)) 808 (frac 1 (bigpos 160 900 47)) 809 (frac -1 (bigpos 600 209 1)) 810 (frac 1 30240) (frac -1 720) 811 (frac 1 12) 1 )) 812 813(defvar math-bernoulli-B-cache '((frac -174611 330) (frac 43867 798) 814 (frac -3617 510) (frac 7 6) (frac -691 2730) 815 (frac 5 66) (frac -1 30) (frac 1 42) 816 (frac -1 30) (frac 1 6) 1 )) 817 818(defvar math-bernoulli-cache-size 11) 819(defun math-bernoulli-coefs (n) 820 (let* ((coefs (list (calcFunc-bern n))) 821 (nn (math-trunc n)) 822 (k nn) 823 (term nn) 824 coef 825 (calc-prefer-frac (or (integerp n) calc-prefer-frac))) 826 (while (>= (setq k (1- k)) 0) 827 (setq term (math-div term (- nn k)) 828 coef (math-mul term (math-bernoulli-number k)) 829 coefs (cons (if (consp n) (math-float coef) coef) coefs) 830 term (math-mul term k))) 831 (nreverse coefs))) 832 833(defun math-bernoulli-number (n) 834 (if (= (% n 2) 1) 835 (if (= n 1) 836 '(frac -1 2) 837 0) 838 (setq n (/ n 2)) 839 (while (>= n math-bernoulli-cache-size) 840 (let* ((sum 0) 841 (nk 1) ; nk = n-k+1 842 (fact 1) ; fact = (n-k+1)! 843 ofact 844 (p math-bernoulli-b-cache) 845 (calc-prefer-frac t)) 846 (math-working "bernoulli B" (* 2 math-bernoulli-cache-size)) 847 (while p 848 (setq nk (+ nk 2) 849 ofact fact 850 fact (math-mul fact (* nk (1- nk))) 851 sum (math-add sum (math-div (car p) fact)) 852 p (cdr p))) 853 (setq ofact (math-mul ofact (1- nk)) 854 sum (math-sub (math-div '(frac 1 2) ofact) sum) 855 math-bernoulli-b-cache (cons sum math-bernoulli-b-cache) 856 math-bernoulli-B-cache (cons (math-mul sum ofact) 857 math-bernoulli-B-cache) 858 math-bernoulli-cache-size (1+ math-bernoulli-cache-size)))) 859 (nth (- math-bernoulli-cache-size n 1) math-bernoulli-B-cache))) 860 861;;; Bn = n! bn 862;;; bn = - sum_k=0^n-1 bk / (n-k+1)! 863 864;;; A faster method would be to use "tangent numbers", c.f., Concrete 865;;; Mathematics pg. 273. 866 867 868;;; Probability distributions. 869 870;;; Binomial. 871(defun calcFunc-utpb (x n p) 872 (if math-expand-formulas 873 (math-normalize (list 'calcFunc-betaI p x (list '+ (list '- n x) 1))) 874 (calcFunc-betaI p x (math-add (math-sub n x) 1)))) 875(put 'calcFunc-utpb 'math-expandable t) 876 877(defun calcFunc-ltpb (x n p) 878 (math-sub 1 (calcFunc-utpb x n p))) 879(put 'calcFunc-ltpb 'math-expandable t) 880 881;;; Chi-square. 882(defun calcFunc-utpc (chisq v) 883 (if math-expand-formulas 884 (math-normalize (list 'calcFunc-gammaQ (list '/ v 2) (list '/ chisq 2))) 885 (calcFunc-gammaQ (math-div v 2) (math-div chisq 2)))) 886(put 'calcFunc-utpc 'math-expandable t) 887 888(defun calcFunc-ltpc (chisq v) 889 (if math-expand-formulas 890 (math-normalize (list 'calcFunc-gammaP (list '/ v 2) (list '/ chisq 2))) 891 (calcFunc-gammaP (math-div v 2) (math-div chisq 2)))) 892(put 'calcFunc-ltpc 'math-expandable t) 893 894;;; F-distribution. 895(defun calcFunc-utpf (f v1 v2) 896 (if math-expand-formulas 897 (math-normalize (list 'calcFunc-betaI 898 (list '/ v2 (list '+ v2 (list '* v1 f))) 899 (list '/ v2 2) 900 (list '/ v1 2))) 901 (calcFunc-betaI (math-div v2 (math-add v2 (math-mul v1 f))) 902 (math-div v2 2) 903 (math-div v1 2)))) 904(put 'calcFunc-utpf 'math-expandable t) 905 906(defun calcFunc-ltpf (f v1 v2) 907 (math-sub 1 (calcFunc-utpf f v1 v2))) 908(put 'calcFunc-ltpf 'math-expandable t) 909 910;;; Normal. 911(defun calcFunc-utpn (x mean sdev) 912 (if math-expand-formulas 913 (math-normalize 914 (list '/ 915 (list '+ 1 916 (list 'calcFunc-erf 917 (list '/ (list '- mean x) 918 (list '* sdev (list 'calcFunc-sqrt 2))))) 919 2)) 920 (math-mul (math-add '(float 1 0) 921 (calcFunc-erf 922 (math-div (math-sub mean x) 923 (math-mul sdev (math-sqrt-2))))) 924 '(float 5 -1)))) 925(put 'calcFunc-utpn 'math-expandable t) 926 927(defun calcFunc-ltpn (x mean sdev) 928 (if math-expand-formulas 929 (math-normalize 930 (list '/ 931 (list '+ 1 932 (list 'calcFunc-erf 933 (list '/ (list '- x mean) 934 (list '* sdev (list 'calcFunc-sqrt 2))))) 935 2)) 936 (math-mul (math-add '(float 1 0) 937 (calcFunc-erf 938 (math-div (math-sub x mean) 939 (math-mul sdev (math-sqrt-2))))) 940 '(float 5 -1)))) 941(put 'calcFunc-ltpn 'math-expandable t) 942 943;;; Poisson. 944(defun calcFunc-utpp (n x) 945 (if math-expand-formulas 946 (math-normalize (list 'calcFunc-gammaP x n)) 947 (calcFunc-gammaP x n))) 948(put 'calcFunc-utpp 'math-expandable t) 949 950(defun calcFunc-ltpp (n x) 951 (if math-expand-formulas 952 (math-normalize (list 'calcFunc-gammaQ x n)) 953 (calcFunc-gammaQ x n))) 954(put 'calcFunc-ltpp 'math-expandable t) 955 956;;; Student's t. (As defined in Abramowitz & Stegun and Numerical Recipes.) 957(defun calcFunc-utpt (tt v) 958 (if math-expand-formulas 959 (math-normalize (list 'calcFunc-betaI 960 (list '/ v (list '+ v (list '^ tt 2))) 961 (list '/ v 2) 962 '(float 5 -1))) 963 (calcFunc-betaI (math-div v (math-add v (math-sqr tt))) 964 (math-div v 2) 965 '(float 5 -1)))) 966(put 'calcFunc-utpt 'math-expandable t) 967 968(defun calcFunc-ltpt (tt v) 969 (math-sub 1 (calcFunc-utpt tt v))) 970(put 'calcFunc-ltpt 'math-expandable t) 971 972(provide 'calc-funcs) 973 974;;; arch-tag: 421ddb7a-550f-4dda-a31c-06638ebfc43a 975;;; calc-funcs.el ends here 976