Mercurial > emacs
comparison lisp/calc/calc-funcs.el @ 41047:73f364fd8aaa
Style cleanup; don't put closing parens on their
own line, add "foo.el ends here" to each file, and update
copyright date.
author | Colin Walters <walters@gnu.org> |
---|---|
date | Wed, 14 Nov 2001 09:09:09 +0000 |
parents | 2fb9d407ae73 |
children | 0759b2de09c1 |
comparison
equal
deleted
inserted
replaced
41046:14b73d89514a | 41047:73f364fd8aaa |
---|---|
1 ;; Calculator for GNU Emacs, part II [calc-funcs.el] | 1 ;; Calculator for GNU Emacs, part II [calc-funcs.el] |
2 ;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc. | 2 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc. |
3 ;; Written by Dave Gillespie, daveg@synaptics.com. | 3 ;; Written by Dave Gillespie, daveg@synaptics.com. |
4 | 4 |
5 ;; This file is part of GNU Emacs. | 5 ;; This file is part of GNU Emacs. |
6 | 6 |
7 ;; GNU Emacs is distributed in the hope that it will be useful, | 7 ;; GNU Emacs is distributed in the hope that it will be useful, |
36 (if (calc-is-hyperbolic) | 36 (if (calc-is-hyperbolic) |
37 (calc-binary-op "gamG" 'calcFunc-gammaG arg) | 37 (calc-binary-op "gamG" 'calcFunc-gammaG arg) |
38 (calc-binary-op "gamQ" 'calcFunc-gammaQ arg)) | 38 (calc-binary-op "gamQ" 'calcFunc-gammaQ arg)) |
39 (if (calc-is-hyperbolic) | 39 (if (calc-is-hyperbolic) |
40 (calc-binary-op "gamg" 'calcFunc-gammag arg) | 40 (calc-binary-op "gamg" 'calcFunc-gammag arg) |
41 (calc-binary-op "gamP" 'calcFunc-gammaP arg)))) | 41 (calc-binary-op "gamP" 'calcFunc-gammaP arg))))) |
42 ) | |
43 | 42 |
44 (defun calc-erf (arg) | 43 (defun calc-erf (arg) |
45 (interactive "P") | 44 (interactive "P") |
46 (calc-slow-wrapper | 45 (calc-slow-wrapper |
47 (if (calc-is-inverse) | 46 (if (calc-is-inverse) |
48 (calc-unary-op "erfc" 'calcFunc-erfc arg) | 47 (calc-unary-op "erfc" 'calcFunc-erfc arg) |
49 (calc-unary-op "erf" 'calcFunc-erf arg))) | 48 (calc-unary-op "erf" 'calcFunc-erf arg)))) |
50 ) | |
51 | 49 |
52 (defun calc-erfc (arg) | 50 (defun calc-erfc (arg) |
53 (interactive "P") | 51 (interactive "P") |
54 (calc-invert-func) | 52 (calc-invert-func) |
55 (calc-erf arg) | 53 (calc-erf arg)) |
56 ) | |
57 | 54 |
58 (defun calc-beta (arg) | 55 (defun calc-beta (arg) |
59 (interactive "P") | 56 (interactive "P") |
60 (calc-slow-wrapper | 57 (calc-slow-wrapper |
61 (calc-binary-op "beta" 'calcFunc-beta arg)) | 58 (calc-binary-op "beta" 'calcFunc-beta arg))) |
62 ) | |
63 | 59 |
64 (defun calc-inc-beta () | 60 (defun calc-inc-beta () |
65 (interactive) | 61 (interactive) |
66 (calc-slow-wrapper | 62 (calc-slow-wrapper |
67 (if (calc-is-hyperbolic) | 63 (if (calc-is-hyperbolic) |
68 (calc-enter-result 3 "betB" (cons 'calcFunc-betaB (calc-top-list-n 3))) | 64 (calc-enter-result 3 "betB" (cons 'calcFunc-betaB (calc-top-list-n 3))) |
69 (calc-enter-result 3 "betI" (cons 'calcFunc-betaI (calc-top-list-n 3))))) | 65 (calc-enter-result 3 "betI" (cons 'calcFunc-betaI (calc-top-list-n 3)))))) |
70 ) | |
71 | 66 |
72 (defun calc-bessel-J (arg) | 67 (defun calc-bessel-J (arg) |
73 (interactive "P") | 68 (interactive "P") |
74 (calc-slow-wrapper | 69 (calc-slow-wrapper |
75 (calc-binary-op "besJ" 'calcFunc-besJ arg)) | 70 (calc-binary-op "besJ" 'calcFunc-besJ arg))) |
76 ) | |
77 | 71 |
78 (defun calc-bessel-Y (arg) | 72 (defun calc-bessel-Y (arg) |
79 (interactive "P") | 73 (interactive "P") |
80 (calc-slow-wrapper | 74 (calc-slow-wrapper |
81 (calc-binary-op "besY" 'calcFunc-besY arg)) | 75 (calc-binary-op "besY" 'calcFunc-besY arg))) |
82 ) | |
83 | 76 |
84 (defun calc-bernoulli-number (arg) | 77 (defun calc-bernoulli-number (arg) |
85 (interactive "P") | 78 (interactive "P") |
86 (calc-slow-wrapper | 79 (calc-slow-wrapper |
87 (if (calc-is-hyperbolic) | 80 (if (calc-is-hyperbolic) |
88 (calc-binary-op "bern" 'calcFunc-bern arg) | 81 (calc-binary-op "bern" 'calcFunc-bern arg) |
89 (calc-unary-op "bern" 'calcFunc-bern arg))) | 82 (calc-unary-op "bern" 'calcFunc-bern arg)))) |
90 ) | |
91 | 83 |
92 (defun calc-euler-number (arg) | 84 (defun calc-euler-number (arg) |
93 (interactive "P") | 85 (interactive "P") |
94 (calc-slow-wrapper | 86 (calc-slow-wrapper |
95 (if (calc-is-hyperbolic) | 87 (if (calc-is-hyperbolic) |
96 (calc-binary-op "eulr" 'calcFunc-euler arg) | 88 (calc-binary-op "eulr" 'calcFunc-euler arg) |
97 (calc-unary-op "eulr" 'calcFunc-euler arg))) | 89 (calc-unary-op "eulr" 'calcFunc-euler arg)))) |
98 ) | |
99 | 90 |
100 (defun calc-stirling-number (arg) | 91 (defun calc-stirling-number (arg) |
101 (interactive "P") | 92 (interactive "P") |
102 (calc-slow-wrapper | 93 (calc-slow-wrapper |
103 (if (calc-is-hyperbolic) | 94 (if (calc-is-hyperbolic) |
104 (calc-binary-op "str2" 'calcFunc-stir2 arg) | 95 (calc-binary-op "str2" 'calcFunc-stir2 arg) |
105 (calc-binary-op "str1" 'calcFunc-stir1 arg))) | 96 (calc-binary-op "str1" 'calcFunc-stir1 arg)))) |
106 ) | |
107 | 97 |
108 (defun calc-utpb () | 98 (defun calc-utpb () |
109 (interactive) | 99 (interactive) |
110 (calc-prob-dist "b" 3) | 100 (calc-prob-dist "b" 3)) |
111 ) | |
112 | 101 |
113 (defun calc-utpc () | 102 (defun calc-utpc () |
114 (interactive) | 103 (interactive) |
115 (calc-prob-dist "c" 2) | 104 (calc-prob-dist "c" 2)) |
116 ) | |
117 | 105 |
118 (defun calc-utpf () | 106 (defun calc-utpf () |
119 (interactive) | 107 (interactive) |
120 (calc-prob-dist "f" 3) | 108 (calc-prob-dist "f" 3)) |
121 ) | |
122 | 109 |
123 (defun calc-utpn () | 110 (defun calc-utpn () |
124 (interactive) | 111 (interactive) |
125 (calc-prob-dist "n" 3) | 112 (calc-prob-dist "n" 3)) |
126 ) | |
127 | 113 |
128 (defun calc-utpp () | 114 (defun calc-utpp () |
129 (interactive) | 115 (interactive) |
130 (calc-prob-dist "p" 2) | 116 (calc-prob-dist "p" 2)) |
131 ) | |
132 | 117 |
133 (defun calc-utpt () | 118 (defun calc-utpt () |
134 (interactive) | 119 (interactive) |
135 (calc-prob-dist "t" 2) | 120 (calc-prob-dist "t" 2)) |
136 ) | |
137 | 121 |
138 (defun calc-prob-dist (letter nargs) | 122 (defun calc-prob-dist (letter nargs) |
139 (calc-slow-wrapper | 123 (calc-slow-wrapper |
140 (if (calc-is-inverse) | 124 (if (calc-is-inverse) |
141 (calc-enter-result nargs (concat "ltp" letter) | 125 (calc-enter-result nargs (concat "ltp" letter) |
143 (calc-top-n 1)) | 127 (calc-top-n 1)) |
144 (calc-top-list-n (1- nargs) 2))) | 128 (calc-top-list-n (1- nargs) 2))) |
145 (calc-enter-result nargs (concat "utp" letter) | 129 (calc-enter-result nargs (concat "utp" letter) |
146 (append (list (intern (concat "calcFunc-utp" letter)) | 130 (append (list (intern (concat "calcFunc-utp" letter)) |
147 (calc-top-n 1)) | 131 (calc-top-n 1)) |
148 (calc-top-list-n (1- nargs) 2))))) | 132 (calc-top-list-n (1- nargs) 2)))))) |
149 ) | |
150 | 133 |
151 | 134 |
152 | 135 |
153 | 136 |
154 ;;; Sources: Numerical Recipes, Press et al; | 137 ;;; Sources: Numerical Recipes, Press et al; |
157 | 140 |
158 ;;; Gamma function. | 141 ;;; Gamma function. |
159 | 142 |
160 (defun calcFunc-gamma (x) | 143 (defun calcFunc-gamma (x) |
161 (or (math-numberp x) (math-reject-arg x 'numberp)) | 144 (or (math-numberp x) (math-reject-arg x 'numberp)) |
162 (calcFunc-fact (math-add x -1)) | 145 (calcFunc-fact (math-add x -1))) |
163 ) | |
164 | 146 |
165 (defun math-gammap1-raw (x &optional fprec nfprec) ; compute gamma(1 + x) | 147 (defun math-gammap1-raw (x &optional fprec nfprec) ; compute gamma(1 + x) |
166 (or fprec | 148 (or fprec |
167 (setq fprec (math-float calc-internal-prec) | 149 (setq fprec (math-float calc-internal-prec) |
168 nfprec (math-float (- calc-internal-prec)))) | 150 nfprec (math-float (- calc-internal-prec)))) |
191 lnx) | 173 lnx) |
192 x) | 174 x) |
193 xinv | 175 xinv |
194 (math-sqr xinv) | 176 (math-sqr xinv) |
195 '(float 0 0) | 177 '(float 0 0) |
196 2)))))) | 178 2))))))) |
197 ) | |
198 | 179 |
199 (defun math-gamma-series (sum x xinvsqr oterm n) | 180 (defun math-gamma-series (sum x xinvsqr oterm n) |
200 (math-working "gamma" sum) | 181 (math-working "gamma" sum) |
201 (let* ((bn (math-bernoulli-number n)) | 182 (let* ((bn (math-bernoulli-number n)) |
202 (term (math-mul (math-div-float (math-float (nth 1 bn)) | 183 (term (math-mul (math-div-float (math-float (nth 1 bn)) |
210 (progn | 191 (progn |
211 ;; Need this because series eventually diverges for large enough n. | 192 ;; Need this because series eventually diverges for large enough n. |
212 (calc-record-why | 193 (calc-record-why |
213 "*Gamma computation stopped early, not all digits may be valid") | 194 "*Gamma computation stopped early, not all digits may be valid") |
214 next) | 195 next) |
215 (math-gamma-series next (math-mul x xinvsqr) xinvsqr term (+ n 2))))) | 196 (math-gamma-series next (math-mul x xinvsqr) xinvsqr term (+ n 2)))))) |
216 ) | |
217 | 197 |
218 | 198 |
219 ;;; Incomplete gamma function. | 199 ;;; Incomplete gamma function. |
220 | 200 |
221 (defun calcFunc-gammaP (a x) | 201 (defun calcFunc-gammaP (a x) |
227 (if (and (math-num-integerp a) | 207 (if (and (math-num-integerp a) |
228 (integerp (setq a (math-trunc a))) | 208 (integerp (setq a (math-trunc a))) |
229 (> a 0) (< a 20)) | 209 (> a 0) (< a 20)) |
230 (math-sub 1 (calcFunc-gammaQ a x)) | 210 (math-sub 1 (calcFunc-gammaQ a x)) |
231 (let ((math-current-gamma-value (calcFunc-gamma a))) | 211 (let ((math-current-gamma-value (calcFunc-gamma a))) |
232 (math-div (calcFunc-gammag a x) math-current-gamma-value)))) | 212 (math-div (calcFunc-gammag a x) math-current-gamma-value))))) |
233 ) | |
234 | 213 |
235 (defun calcFunc-gammaQ (a x) | 214 (defun calcFunc-gammaQ (a x) |
236 (if (equal x '(var inf var-inf)) | 215 (if (equal x '(var inf var-inf)) |
237 '(float 0 0) | 216 '(float 0 0) |
238 (math-inexact-result) | 217 (math-inexact-result) |
249 (setq term (math-div (math-mul term x) n) | 228 (setq term (math-div (math-mul term x) n) |
250 sum (math-add sum term)) | 229 sum (math-add sum term)) |
251 (math-working "gamma" sum)) | 230 (math-working "gamma" sum)) |
252 (math-mul sum (calcFunc-exp (math-neg x))))) | 231 (math-mul sum (calcFunc-exp (math-neg x))))) |
253 (let ((math-current-gamma-value (calcFunc-gamma a))) | 232 (let ((math-current-gamma-value (calcFunc-gamma a))) |
254 (math-div (calcFunc-gammaG a x) math-current-gamma-value)))) | 233 (math-div (calcFunc-gammaG a x) math-current-gamma-value))))) |
255 ) | |
256 | 234 |
257 (defun calcFunc-gammag (a x) | 235 (defun calcFunc-gammag (a x) |
258 (if (equal x '(var inf var-inf)) | 236 (if (equal x '(var inf var-inf)) |
259 (calcFunc-gamma a) | 237 (calcFunc-gamma a) |
260 (math-inexact-result) | 238 (math-inexact-result) |
267 (math-lessp-float (calcFunc-re x) | 245 (math-lessp-float (calcFunc-re x) |
268 (math-add-float (calcFunc-re a) | 246 (math-add-float (calcFunc-re a) |
269 '(float 1 0)))) | 247 '(float 1 0)))) |
270 (math-inc-gamma-series a x) | 248 (math-inc-gamma-series a x) |
271 (math-sub (or math-current-gamma-value (calcFunc-gamma a)) | 249 (math-sub (or math-current-gamma-value (calcFunc-gamma a)) |
272 (math-inc-gamma-cfrac a x))))) | 250 (math-inc-gamma-cfrac a x)))))) |
273 ) | |
274 (setq math-current-gamma-value nil) | 251 (setq math-current-gamma-value nil) |
275 | 252 |
276 (defun calcFunc-gammaG (a x) | 253 (defun calcFunc-gammaG (a x) |
277 (if (equal x '(var inf var-inf)) | 254 (if (equal x '(var inf var-inf)) |
278 '(float 0 0) | 255 '(float 0 0) |
286 (math-lessp-float (calcFunc-re x) | 263 (math-lessp-float (calcFunc-re x) |
287 (math-add-float (math-abs-approx a) | 264 (math-add-float (math-abs-approx a) |
288 '(float 1 0)))) | 265 '(float 1 0)))) |
289 (math-sub (or math-current-gamma-value (calcFunc-gamma a)) | 266 (math-sub (or math-current-gamma-value (calcFunc-gamma a)) |
290 (math-inc-gamma-series a x)) | 267 (math-inc-gamma-series a x)) |
291 (math-inc-gamma-cfrac a x)))) | 268 (math-inc-gamma-cfrac a x))))) |
292 ) | |
293 | 269 |
294 (defun math-inc-gamma-series (a x) | 270 (defun math-inc-gamma-series (a x) |
295 (if (Math-zerop x) | 271 (if (Math-zerop x) |
296 '(float 0 0) | 272 '(float 0 0) |
297 (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x)) | 273 (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x)) |
298 (math-with-extra-prec 2 | 274 (math-with-extra-prec 2 |
299 (let ((start (math-div '(float 1 0) a))) | 275 (let ((start (math-div '(float 1 0) a))) |
300 (math-inc-gamma-series-step start start a x))))) | 276 (math-inc-gamma-series-step start start a x)))))) |
301 ) | |
302 | 277 |
303 (defun math-inc-gamma-series-step (sum term a x) | 278 (defun math-inc-gamma-series-step (sum term a x) |
304 (math-working "gamma" sum) | 279 (math-working "gamma" sum) |
305 (setq a (math-add a '(float 1 0)) | 280 (setq a (math-add a '(float 1 0)) |
306 term (math-div (math-mul term x) a)) | 281 term (math-div (math-mul term x) a)) |
307 (let ((next (math-add sum term))) | 282 (let ((next (math-add sum term))) |
308 (if (math-nearly-equal sum next) | 283 (if (math-nearly-equal sum next) |
309 next | 284 next |
310 (math-inc-gamma-series-step next term a x))) | 285 (math-inc-gamma-series-step next term a x)))) |
311 ) | |
312 | 286 |
313 (defun math-inc-gamma-cfrac (a x) | 287 (defun math-inc-gamma-cfrac (a x) |
314 (if (Math-zerop x) | 288 (if (Math-zerop x) |
315 (or math-current-gamma-value (calcFunc-gamma a)) | 289 (or math-current-gamma-value (calcFunc-gamma a)) |
316 (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x)) | 290 (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x)) |
317 (math-inc-gamma-cfrac-step '(float 1 0) x | 291 (math-inc-gamma-cfrac-step '(float 1 0) x |
318 '(float 0 0) '(float 1 0) | 292 '(float 0 0) '(float 1 0) |
319 '(float 1 0) '(float 1 0) '(float 0 0) | 293 '(float 1 0) '(float 1 0) '(float 0 0) |
320 a x))) | 294 a x)))) |
321 ) | |
322 | 295 |
323 (defun math-inc-gamma-cfrac-step (a0 a1 b0 b1 n fac g a x) | 296 (defun math-inc-gamma-cfrac-step (a0 a1 b0 b1 n fac g a x) |
324 (let ((ana (math-sub n a)) | 297 (let ((ana (math-sub n a)) |
325 (anf (math-mul n fac))) | 298 (anf (math-mul n fac))) |
326 (setq n (math-add n '(float 1 0)) | 299 (setq n (math-add n '(float 1 0)) |
333 (setq fac (math-div '(float 1 0) a1)) | 306 (setq fac (math-div '(float 1 0) a1)) |
334 (let ((next (math-mul b1 fac))) | 307 (let ((next (math-mul b1 fac))) |
335 (math-working "gamma" next) | 308 (math-working "gamma" next) |
336 (if (math-nearly-equal next g) | 309 (if (math-nearly-equal next g) |
337 next | 310 next |
338 (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac next a x))))) | 311 (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac next a x)))))) |
339 ) | |
340 | 312 |
341 | 313 |
342 ;;; Error function. | 314 ;;; Error function. |
343 | 315 |
344 (defun calcFunc-erf (x) | 316 (defun calcFunc-erf (x) |
351 (let ((math-current-gamma-value (math-sqrt-pi))) | 323 (let ((math-current-gamma-value (math-sqrt-pi))) |
352 (math-to-same-complex-quad | 324 (math-to-same-complex-quad |
353 (math-div (calcFunc-gammag '(float 5 -1) | 325 (math-div (calcFunc-gammag '(float 5 -1) |
354 (math-sqr (math-to-complex-quad-one x))) | 326 (math-sqr (math-to-complex-quad-one x))) |
355 math-current-gamma-value) | 327 math-current-gamma-value) |
356 x))))) | 328 x)))))) |
357 ) | |
358 | 329 |
359 (defun calcFunc-erfc (x) | 330 (defun calcFunc-erfc (x) |
360 (if (equal x '(var inf var-inf)) | 331 (if (equal x '(var inf var-inf)) |
361 '(float 0 0) | 332 '(float 0 0) |
362 (if (math-posp x) | 333 (if (math-posp x) |
363 (let ((math-current-gamma-value (math-sqrt-pi))) | 334 (let ((math-current-gamma-value (math-sqrt-pi))) |
364 (math-div (calcFunc-gammaG '(float 5 -1) (math-sqr x)) | 335 (math-div (calcFunc-gammaG '(float 5 -1) (math-sqr x)) |
365 math-current-gamma-value)) | 336 math-current-gamma-value)) |
366 (math-sub 1 (calcFunc-erf x)))) | 337 (math-sub 1 (calcFunc-erf x))))) |
367 ) | |
368 | 338 |
369 (defun math-to-complex-quad-one (x) | 339 (defun math-to-complex-quad-one (x) |
370 (if (eq (car-safe x) 'polar) (setq x (math-complex x))) | 340 (if (eq (car-safe x) 'polar) (setq x (math-complex x))) |
371 (if (eq (car-safe x) 'cplx) | 341 (if (eq (car-safe x) 'cplx) |
372 (list 'cplx (math-abs (nth 1 x)) (math-abs (nth 2 x))) | 342 (list 'cplx (math-abs (nth 1 x)) (math-abs (nth 2 x))) |
373 x) | 343 x)) |
374 ) | |
375 | 344 |
376 (defun math-to-same-complex-quad (x y) | 345 (defun math-to-same-complex-quad (x y) |
377 (if (eq (car-safe y) 'cplx) | 346 (if (eq (car-safe y) 'cplx) |
378 (if (eq (car-safe x) 'cplx) | 347 (if (eq (car-safe x) 'cplx) |
379 (list 'cplx | 348 (list 'cplx |
382 (if (math-negp (nth 1 y)) (math-neg x) x)) | 351 (if (math-negp (nth 1 y)) (math-neg x) x)) |
383 (if (math-negp y) | 352 (if (math-negp y) |
384 (if (eq (car-safe x) 'cplx) | 353 (if (eq (car-safe x) 'cplx) |
385 (list 'cplx (math-neg (nth 1 x)) (nth 2 x)) | 354 (list 'cplx (math-neg (nth 1 x)) (nth 2 x)) |
386 (math-neg x)) | 355 (math-neg x)) |
387 x)) | 356 x))) |
388 ) | |
389 | 357 |
390 | 358 |
391 ;;; Beta function. | 359 ;;; Beta function. |
392 | 360 |
393 (defun calcFunc-beta (a b) | 361 (defun calcFunc-beta (a b) |
396 (or (math-numberp b) (math-reject-arg b 'numberp)) | 364 (or (math-numberp b) (math-reject-arg b 'numberp)) |
397 (math-div 1 (math-mul b (calcFunc-choose (math-add b am) am)))) | 365 (math-div 1 (math-mul b (calcFunc-choose (math-add b am) am)))) |
398 (if (math-num-integerp b) | 366 (if (math-num-integerp b) |
399 (calcFunc-beta b a) | 367 (calcFunc-beta b a) |
400 (math-div (math-mul (calcFunc-gamma a) (calcFunc-gamma b)) | 368 (math-div (math-mul (calcFunc-gamma a) (calcFunc-gamma b)) |
401 (calcFunc-gamma (math-add a b))))) | 369 (calcFunc-gamma (math-add a b)))))) |
402 ) | |
403 | 370 |
404 | 371 |
405 ;;; Incomplete beta function. | 372 ;;; Incomplete beta function. |
406 | 373 |
407 (defun calcFunc-betaI (x a b) | 374 (defun calcFunc-betaI (x a b) |
423 '(float 0 0)) | 390 '(float 0 0)) |
424 ((not (math-numberp a)) (math-reject-arg a 'numberp)) | 391 ((not (math-numberp a)) (math-reject-arg a 'numberp)) |
425 ((not (math-numberp b)) (math-reject-arg b 'numberp)) | 392 ((not (math-numberp b)) (math-reject-arg b 'numberp)) |
426 ((math-inexact-result)) | 393 ((math-inexact-result)) |
427 (t (let ((math-current-beta-value (calcFunc-beta a b))) | 394 (t (let ((math-current-beta-value (calcFunc-beta a b))) |
428 (math-div (calcFunc-betaB x a b) math-current-beta-value)))) | 395 (math-div (calcFunc-betaB x a b) math-current-beta-value))))) |
429 ) | |
430 | 396 |
431 (defun calcFunc-betaB (x a b) | 397 (defun calcFunc-betaB (x a b) |
432 (cond | 398 (cond |
433 ((math-zerop x) | 399 ((math-zerop x) |
434 '(float 0 0)) | 400 '(float 0 0)) |
476 (math-add (math-add a b) '(float 2 0)))) | 442 (math-add (math-add a b) '(float 2 0)))) |
477 (math-div (math-mul bt (math-beta-cfrac a b x)) a) | 443 (math-div (math-mul bt (math-beta-cfrac a b x)) a) |
478 (math-sub (or math-current-beta-value (calcFunc-beta a b)) | 444 (math-sub (or math-current-beta-value (calcFunc-beta a b)) |
479 (math-div (math-mul bt | 445 (math-div (math-mul bt |
480 (math-beta-cfrac b a (math-sub 1 x))) | 446 (math-beta-cfrac b a (math-sub 1 x))) |
481 b))))))) | 447 b)))))))) |
482 ) | |
483 (setq math-current-beta-value nil) | 448 (setq math-current-beta-value nil) |
484 | 449 |
485 (defun math-beta-cfrac (a b x) | 450 (defun math-beta-cfrac (a b x) |
486 (let ((qab (math-add a b)) | 451 (let ((qab (math-add a b)) |
487 (qap (math-add a '(float 1 0))) | 452 (qap (math-add a '(float 1 0))) |
489 (math-beta-cfrac-step '(float 1 0) | 454 (math-beta-cfrac-step '(float 1 0) |
490 (math-sub '(float 1 0) | 455 (math-sub '(float 1 0) |
491 (math-div (math-mul qab x) qap)) | 456 (math-div (math-mul qab x) qap)) |
492 '(float 1 0) '(float 1 0) | 457 '(float 1 0) '(float 1 0) |
493 '(float 1 0) | 458 '(float 1 0) |
494 qab qap qam a b x)) | 459 qab qap qam a b x))) |
495 ) | |
496 | 460 |
497 (defun math-beta-cfrac-step (az bz am bm m qab qap qam a b x) | 461 (defun math-beta-cfrac-step (az bz am bm m qab qap qam a b x) |
498 (let* ((two-m (math-mul m '(float 2 0))) | 462 (let* ((two-m (math-mul m '(float 2 0))) |
499 (d (math-div (math-mul (math-mul (math-sub b m) m) x) | 463 (d (math-div (math-mul (math-mul (math-sub b m) m) x) |
500 (math-mul (math-add qam two-m) (math-add a two-m)))) | 464 (math-mul (math-add qam two-m) (math-add a two-m)))) |
510 (if (math-nearly-equal next az) | 474 (if (math-nearly-equal next az) |
511 next | 475 next |
512 (math-beta-cfrac-step next '(float 1 0) | 476 (math-beta-cfrac-step next '(float 1 0) |
513 (math-div ap bpp) (math-div bp bpp) | 477 (math-div ap bpp) (math-div bp bpp) |
514 (math-add m '(float 1 0)) | 478 (math-add m '(float 1 0)) |
515 qab qap qam a b x))) | 479 qab qap qam a b x)))) |
516 ) | |
517 | 480 |
518 | 481 |
519 ;;; Bessel functions. | 482 ;;; Bessel functions. |
520 | 483 |
521 ;;; Should generalize this to handle arbitrary precision! | 484 ;;; Should generalize this to handle arbitrary precision! |
581 sum (math-mul sum '(float 1 -10)))) | 544 sum (math-mul sum '(float 1 -10)))) |
582 (or (setq jsum (not jsum)) | 545 (or (setq jsum (not jsum)) |
583 (setq sum (math-add sum bj))) | 546 (setq sum (math-add sum bj))) |
584 (if (= j v) | 547 (if (= j v) |
585 (setq ans bjp))) | 548 (setq ans bjp))) |
586 (math-div ans (math-sub (math-mul 2 sum) bj))))))) | 549 (math-div ans (math-sub (math-mul 2 sum) bj)))))))) |
587 ) | |
588 | 550 |
589 (defun math-besJ-series (sum term k zz vk) | 551 (defun math-besJ-series (sum term k zz vk) |
590 (math-working "besJ" sum) | 552 (math-working "besJ" sum) |
591 (setq k (1+ k) | 553 (setq k (1+ k) |
592 vk (math-add 1 vk) | 554 vk (math-add 1 vk) |
593 term (math-div (math-mul term zz) (math-mul k vk))) | 555 term (math-div (math-mul term zz) (math-mul k vk))) |
594 (let ((next (math-add sum term))) | 556 (let ((next (math-add sum term))) |
595 (if (math-nearly-equal next sum) | 557 (if (math-nearly-equal next sum) |
596 next | 558 next |
597 (math-besJ-series next term k zz vk))) | 559 (math-besJ-series next term k zz vk)))) |
598 ) | |
599 | 560 |
600 (defun math-besJ0 (x &optional yflag) | 561 (defun math-besJ0 (x &optional yflag) |
601 (cond ((and (not yflag) (math-negp (calcFunc-re x))) | 562 (cond ((and (not yflag) (math-negp (calcFunc-re x))) |
602 (math-besJ0 (math-neg x))) | 563 (math-besJ0 (math-neg x))) |
603 ((Math-lessp '(float 8 0) (math-abs-approx x)) | 564 ((Math-lessp '(float 8 0) (math-abs-approx x)) |
636 '((float 1 0) | 597 '((float 1 0) |
637 (float (bigpos 712 532 678 2) -7) | 598 (float (bigpos 712 532 678 2) -7) |
638 (float (bigpos 853 264 927 5) -5) | 599 (float (bigpos 853 264 927 5) -5) |
639 (float (bigpos 718 680 494 9) -3) | 600 (float (bigpos 718 680 494 9) -3) |
640 (float (bigpos 985 532 029 1) 0) | 601 (float (bigpos 985 532 029 1) 0) |
641 (float (bigpos 411 490 568 57) 0))))))) | 602 (float (bigpos 411 490 568 57) 0)))))))) |
642 ) | |
643 | 603 |
644 (defun math-besJ1 (x &optional yflag) | 604 (defun math-besJ1 (x &optional yflag) |
645 (cond ((and (math-negp (calcFunc-re x)) (not yflag)) | 605 (cond ((and (math-negp (calcFunc-re x)) (not yflag)) |
646 (math-neg (math-besJ1 (math-neg x)))) | 606 (math-neg (math-besJ1 (math-neg x)))) |
647 ((Math-lessp '(float 8 0) (math-abs-approx x)) | 607 ((Math-lessp '(float 8 0) (math-abs-approx x)) |
684 (float (bigpos 397 991 769 3) -7) | 644 (float (bigpos 397 991 769 3) -7) |
685 (float (bigpos 394 743 944 9) -5) | 645 (float (bigpos 394 743 944 9) -5) |
686 (float (bigpos 474 330 858 1) -2) | 646 (float (bigpos 474 330 858 1) -2) |
687 (float (bigpos 178 535 300 2) 0) | 647 (float (bigpos 178 535 300 2) 0) |
688 (float (bigpos 442 228 725 144) | 648 (float (bigpos 442 228 725 144) |
689 0)))))))) | 649 0))))))))) |
690 ) | |
691 | 650 |
692 (defun calcFunc-besY (v x) | 651 (defun calcFunc-besY (v x) |
693 (math-inexact-result) | 652 (math-inexact-result) |
694 (or (math-numberp v) (math-reject-arg v 'numberp)) | 653 (or (math-numberp v) (math-reject-arg v 'numberp)) |
695 (or (math-numberp x) (math-reject-arg x 'numberp)) | 654 (or (math-numberp x) (math-reject-arg x 'numberp)) |
719 (while (< (setq j (1+ j)) v) | 678 (while (< (setq j (1+ j)) v) |
720 (setq byp (math-sub (math-mul (math-mul j two-over-x) by) | 679 (setq byp (math-sub (math-mul (math-mul j two-over-x) by) |
721 bym) | 680 bym) |
722 bym by | 681 bym by |
723 by byp)) | 682 by byp)) |
724 by))))) | 683 by)))))) |
725 ) | |
726 | 684 |
727 (defun math-besY0 (x) | 685 (defun math-besY0 (x) |
728 (cond ((Math-lessp (math-abs-approx x) '(float 8 0)) | 686 (cond ((Math-lessp (math-abs-approx x) '(float 8 0)) |
729 (let ((y (math-sqr x))) | 687 (let ((y (math-sqr x))) |
730 (math-add | 688 (math-add |
747 ((math-negp (calcFunc-re x)) | 705 ((math-negp (calcFunc-re x)) |
748 (math-add (math-besJ0 (math-neg x) t) | 706 (math-add (math-besJ0 (math-neg x) t) |
749 (math-mul '(cplx 0 2) | 707 (math-mul '(cplx 0 2) |
750 (math-besJ0 (math-neg x))))) | 708 (math-besJ0 (math-neg x))))) |
751 (t | 709 (t |
752 (math-besJ0 x t))) | 710 (math-besJ0 x t)))) |
753 ) | |
754 | 711 |
755 (defun math-besY1 (x) | 712 (defun math-besY1 (x) |
756 (cond ((Math-lessp (math-abs-approx x) '(float 8 0)) | 713 (cond ((Math-lessp (math-abs-approx x) '(float 8 0)) |
757 (let ((y (math-sqr x))) | 714 (let ((y (math-sqr x))) |
758 (math-add | 715 (math-add |
780 (math-neg | 737 (math-neg |
781 (math-add (math-besJ1 (math-neg x) t) | 738 (math-add (math-besJ1 (math-neg x) t) |
782 (math-mul '(cplx 0 2) | 739 (math-mul '(cplx 0 2) |
783 (math-besJ1 (math-neg x)))))) | 740 (math-besJ1 (math-neg x)))))) |
784 (t | 741 (t |
785 (math-besJ1 x t))) | 742 (math-besJ1 x t)))) |
786 ) | |
787 | 743 |
788 (defun math-poly-eval (x coefs) | 744 (defun math-poly-eval (x coefs) |
789 (let ((accum (car coefs))) | 745 (let ((accum (car coefs))) |
790 (while (setq coefs (cdr coefs)) | 746 (while (setq coefs (cdr coefs)) |
791 (setq accum (math-add (car coefs) (math-mul accum x)))) | 747 (setq accum (math-add (car coefs) (math-mul accum x)))) |
792 accum) | 748 accum)) |
793 ) | |
794 | 749 |
795 | 750 |
796 ;;;; Bernoulli and Euler polynomials and numbers. | 751 ;;;; Bernoulli and Euler polynomials and numbers. |
797 | 752 |
798 (defun calcFunc-bern (n &optional x) | 753 (defun calcFunc-bern (n &optional x) |
803 (or (math-num-natnump n) (math-reject-arg n 'natnump)) | 758 (or (math-num-natnump n) (math-reject-arg n 'natnump)) |
804 (if (consp n) | 759 (if (consp n) |
805 (progn | 760 (progn |
806 (math-inexact-result) | 761 (math-inexact-result) |
807 (math-float (math-bernoulli-number (math-trunc n)))) | 762 (math-float (math-bernoulli-number (math-trunc n)))) |
808 (math-bernoulli-number n))) | 763 (math-bernoulli-number n)))) |
809 ) | |
810 | 764 |
811 (defun calcFunc-euler (n &optional x) | 765 (defun calcFunc-euler (n &optional x) |
812 (or (math-num-natnump n) (math-reject-arg n 'natnump)) | 766 (or (math-num-natnump n) (math-reject-arg n 'natnump)) |
813 (if x | 767 (if x |
814 (let* ((n1 (math-add n 1)) | 768 (let* ((n1 (math-add n 1)) |
838 (math-mul (math-pow 2 n) | 792 (math-mul (math-pow 2 n) |
839 (if (consp n) | 793 (if (consp n) |
840 (progn | 794 (progn |
841 (math-inexact-result) | 795 (math-inexact-result) |
842 (calcFunc-euler n '(float 5 -1))) | 796 (calcFunc-euler n '(float 5 -1))) |
843 (calcFunc-euler n '(frac 1 2))))) | 797 (calcFunc-euler n '(frac 1 2)))))) |
844 ) | |
845 | 798 |
846 (defun math-bernoulli-coefs (n) | 799 (defun math-bernoulli-coefs (n) |
847 (let* ((coefs (list (calcFunc-bern n))) | 800 (let* ((coefs (list (calcFunc-bern n))) |
848 (nn (math-trunc n)) | 801 (nn (math-trunc n)) |
849 (k nn) | 802 (k nn) |
853 (while (>= (setq k (1- k)) 0) | 806 (while (>= (setq k (1- k)) 0) |
854 (setq term (math-div term (- nn k)) | 807 (setq term (math-div term (- nn k)) |
855 coef (math-mul term (math-bernoulli-number k)) | 808 coef (math-mul term (math-bernoulli-number k)) |
856 coefs (cons (if (consp n) (math-float coef) coef) coefs) | 809 coefs (cons (if (consp n) (math-float coef) coef) coefs) |
857 term (math-mul term k))) | 810 term (math-mul term k))) |
858 (nreverse coefs)) | 811 (nreverse coefs))) |
859 ) | |
860 | 812 |
861 (defun math-bernoulli-number (n) | 813 (defun math-bernoulli-number (n) |
862 (if (= (% n 2) 1) | 814 (if (= (% n 2) 1) |
863 (if (= n 1) | 815 (if (= n 1) |
864 '(frac -1 2) | 816 '(frac -1 2) |
882 sum (math-sub (math-div '(frac 1 2) ofact) sum) | 834 sum (math-sub (math-div '(frac 1 2) ofact) sum) |
883 math-bernoulli-b-cache (cons sum math-bernoulli-b-cache) | 835 math-bernoulli-b-cache (cons sum math-bernoulli-b-cache) |
884 math-bernoulli-B-cache (cons (math-mul sum ofact) | 836 math-bernoulli-B-cache (cons (math-mul sum ofact) |
885 math-bernoulli-B-cache) | 837 math-bernoulli-B-cache) |
886 math-bernoulli-cache-size (1+ math-bernoulli-cache-size)))) | 838 math-bernoulli-cache-size (1+ math-bernoulli-cache-size)))) |
887 (nth (- math-bernoulli-cache-size n 1) math-bernoulli-B-cache)) | 839 (nth (- math-bernoulli-cache-size n 1) math-bernoulli-B-cache))) |
888 ) | |
889 | 840 |
890 ;;; Bn = n! bn | 841 ;;; Bn = n! bn |
891 ;;; bn = - sum_k=0^n-1 bk / (n-k+1)! | 842 ;;; bn = - sum_k=0^n-1 bk / (n-k+1)! |
892 | 843 |
893 ;;; A faster method would be to use "tangent numbers", c.f., Concrete | 844 ;;; A faster method would be to use "tangent numbers", c.f., Concrete |
917 | 868 |
918 ;;; Binomial. | 869 ;;; Binomial. |
919 (defun calcFunc-utpb (x n p) | 870 (defun calcFunc-utpb (x n p) |
920 (if math-expand-formulas | 871 (if math-expand-formulas |
921 (math-normalize (list 'calcFunc-betaI p x (list '+ (list '- n x) 1))) | 872 (math-normalize (list 'calcFunc-betaI p x (list '+ (list '- n x) 1))) |
922 (calcFunc-betaI p x (math-add (math-sub n x) 1))) | 873 (calcFunc-betaI p x (math-add (math-sub n x) 1)))) |
923 ) | |
924 (put 'calcFunc-utpb 'math-expandable t) | 874 (put 'calcFunc-utpb 'math-expandable t) |
925 | 875 |
926 (defun calcFunc-ltpb (x n p) | 876 (defun calcFunc-ltpb (x n p) |
927 (math-sub 1 (calcFunc-utpb x n p)) | 877 (math-sub 1 (calcFunc-utpb x n p))) |
928 ) | |
929 (put 'calcFunc-ltpb 'math-expandable t) | 878 (put 'calcFunc-ltpb 'math-expandable t) |
930 | 879 |
931 ;;; Chi-square. | 880 ;;; Chi-square. |
932 (defun calcFunc-utpc (chisq v) | 881 (defun calcFunc-utpc (chisq v) |
933 (if math-expand-formulas | 882 (if math-expand-formulas |
934 (math-normalize (list 'calcFunc-gammaQ (list '/ v 2) (list '/ chisq 2))) | 883 (math-normalize (list 'calcFunc-gammaQ (list '/ v 2) (list '/ chisq 2))) |
935 (calcFunc-gammaQ (math-div v 2) (math-div chisq 2))) | 884 (calcFunc-gammaQ (math-div v 2) (math-div chisq 2)))) |
936 ) | |
937 (put 'calcFunc-utpc 'math-expandable t) | 885 (put 'calcFunc-utpc 'math-expandable t) |
938 | 886 |
939 (defun calcFunc-ltpc (chisq v) | 887 (defun calcFunc-ltpc (chisq v) |
940 (if math-expand-formulas | 888 (if math-expand-formulas |
941 (math-normalize (list 'calcFunc-gammaP (list '/ v 2) (list '/ chisq 2))) | 889 (math-normalize (list 'calcFunc-gammaP (list '/ v 2) (list '/ chisq 2))) |
942 (calcFunc-gammaP (math-div v 2) (math-div chisq 2))) | 890 (calcFunc-gammaP (math-div v 2) (math-div chisq 2)))) |
943 ) | |
944 (put 'calcFunc-ltpc 'math-expandable t) | 891 (put 'calcFunc-ltpc 'math-expandable t) |
945 | 892 |
946 ;;; F-distribution. | 893 ;;; F-distribution. |
947 (defun calcFunc-utpf (f v1 v2) | 894 (defun calcFunc-utpf (f v1 v2) |
948 (if math-expand-formulas | 895 (if math-expand-formulas |
950 (list '/ v2 (list '+ v2 (list '* v1 f))) | 897 (list '/ v2 (list '+ v2 (list '* v1 f))) |
951 (list '/ v2 2) | 898 (list '/ v2 2) |
952 (list '/ v1 2))) | 899 (list '/ v1 2))) |
953 (calcFunc-betaI (math-div v2 (math-add v2 (math-mul v1 f))) | 900 (calcFunc-betaI (math-div v2 (math-add v2 (math-mul v1 f))) |
954 (math-div v2 2) | 901 (math-div v2 2) |
955 (math-div v1 2))) | 902 (math-div v1 2)))) |
956 ) | |
957 (put 'calcFunc-utpf 'math-expandable t) | 903 (put 'calcFunc-utpf 'math-expandable t) |
958 | 904 |
959 (defun calcFunc-ltpf (f v1 v2) | 905 (defun calcFunc-ltpf (f v1 v2) |
960 (math-sub 1 (calcFunc-utpf f v1 v2)) | 906 (math-sub 1 (calcFunc-utpf f v1 v2))) |
961 ) | |
962 (put 'calcFunc-ltpf 'math-expandable t) | 907 (put 'calcFunc-ltpf 'math-expandable t) |
963 | 908 |
964 ;;; Normal. | 909 ;;; Normal. |
965 (defun calcFunc-utpn (x mean sdev) | 910 (defun calcFunc-utpn (x mean sdev) |
966 (if math-expand-formulas | 911 (if math-expand-formulas |
973 2)) | 918 2)) |
974 (math-mul (math-add '(float 1 0) | 919 (math-mul (math-add '(float 1 0) |
975 (calcFunc-erf | 920 (calcFunc-erf |
976 (math-div (math-sub mean x) | 921 (math-div (math-sub mean x) |
977 (math-mul sdev (math-sqrt-2))))) | 922 (math-mul sdev (math-sqrt-2))))) |
978 '(float 5 -1))) | 923 '(float 5 -1)))) |
979 ) | |
980 (put 'calcFunc-utpn 'math-expandable t) | 924 (put 'calcFunc-utpn 'math-expandable t) |
981 | 925 |
982 (defun calcFunc-ltpn (x mean sdev) | 926 (defun calcFunc-ltpn (x mean sdev) |
983 (if math-expand-formulas | 927 (if math-expand-formulas |
984 (math-normalize | 928 (math-normalize |
990 2)) | 934 2)) |
991 (math-mul (math-add '(float 1 0) | 935 (math-mul (math-add '(float 1 0) |
992 (calcFunc-erf | 936 (calcFunc-erf |
993 (math-div (math-sub x mean) | 937 (math-div (math-sub x mean) |
994 (math-mul sdev (math-sqrt-2))))) | 938 (math-mul sdev (math-sqrt-2))))) |
995 '(float 5 -1))) | 939 '(float 5 -1)))) |
996 ) | |
997 (put 'calcFunc-ltpn 'math-expandable t) | 940 (put 'calcFunc-ltpn 'math-expandable t) |
998 | 941 |
999 ;;; Poisson. | 942 ;;; Poisson. |
1000 (defun calcFunc-utpp (n x) | 943 (defun calcFunc-utpp (n x) |
1001 (if math-expand-formulas | 944 (if math-expand-formulas |
1002 (math-normalize (list 'calcFunc-gammaP x n)) | 945 (math-normalize (list 'calcFunc-gammaP x n)) |
1003 (calcFunc-gammaP x n)) | 946 (calcFunc-gammaP x n))) |
1004 ) | |
1005 (put 'calcFunc-utpp 'math-expandable t) | 947 (put 'calcFunc-utpp 'math-expandable t) |
1006 | 948 |
1007 (defun calcFunc-ltpp (n x) | 949 (defun calcFunc-ltpp (n x) |
1008 (if math-expand-formulas | 950 (if math-expand-formulas |
1009 (math-normalize (list 'calcFunc-gammaQ x n)) | 951 (math-normalize (list 'calcFunc-gammaQ x n)) |
1010 (calcFunc-gammaQ x n)) | 952 (calcFunc-gammaQ x n))) |
1011 ) | |
1012 (put 'calcFunc-ltpp 'math-expandable t) | 953 (put 'calcFunc-ltpp 'math-expandable t) |
1013 | 954 |
1014 ;;; Student's t. (As defined in Abramowitz & Stegun and Numerical Recipes.) | 955 ;;; Student's t. (As defined in Abramowitz & Stegun and Numerical Recipes.) |
1015 (defun calcFunc-utpt (tt v) | 956 (defun calcFunc-utpt (tt v) |
1016 (if math-expand-formulas | 957 (if math-expand-formulas |
1018 (list '/ v (list '+ v (list '^ tt 2))) | 959 (list '/ v (list '+ v (list '^ tt 2))) |
1019 (list '/ v 2) | 960 (list '/ v 2) |
1020 '(float 5 -1))) | 961 '(float 5 -1))) |
1021 (calcFunc-betaI (math-div v (math-add v (math-sqr tt))) | 962 (calcFunc-betaI (math-div v (math-add v (math-sqr tt))) |
1022 (math-div v 2) | 963 (math-div v 2) |
1023 '(float 5 -1))) | 964 '(float 5 -1)))) |
1024 ) | |
1025 (put 'calcFunc-utpt 'math-expandable t) | 965 (put 'calcFunc-utpt 'math-expandable t) |
1026 | 966 |
1027 (defun calcFunc-ltpt (tt v) | 967 (defun calcFunc-ltpt (tt v) |
1028 (math-sub 1 (calcFunc-utpt tt v)) | 968 (math-sub 1 (calcFunc-utpt tt v))) |
1029 ) | |
1030 (put 'calcFunc-ltpt 'math-expandable t) | 969 (put 'calcFunc-ltpt 'math-expandable t) |
1031 | 970 |
1032 | 971 |
1033 | 972 ;;; calc-funcs.el ends here |
1034 |