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