calc-funcs.el

Revision 4220, 29.8 kB (checked in by miyoshi, 9 months ago)
Sync up with Emacs22.2.

Line
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, 2008 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 3, 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