-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathappell.lisp
145 lines (120 loc) · 4.33 KB
/
appell.lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
;; not done:
;; (a) conjugate function
;; (b) TeX support
;; (c) Power series
;; (d) wxMaxima display code
;; (e) sign-function (unlikely anything useful)
;; (f) rbp & rbp ?
;; (g) limit code
;; (h) bigfloat numerical evaluation
;; (i) msimpind (Done)
(in-package :maxima)
(defmvar $use_appell_f1_integral_rep nil
"When 'true', the general simplifier returns an integral representation of 'appell_f1"
:properties ((assign #'(lambda (s a)
(if (or (eq a nil) (eq a t))
t
(mtell "The value of ~M must be either 't' or 'false, but found ~M . ~% " s a))))))
(def-simplifier appell_f1 (a b1 b2 c x y)
(cond
((eql 0 y)
(ftake '%hypergeometric (ftake 'mlist a b1) (ftake 'mlist c) x))
((eql 0 x)
(ftake '%hypergeometric (ftake 'mlist a b2) (ftake 'mlist c) y))
((alike1 x y)
(ftake '%hypergeometric (ftake 'mlist a (add b1 b2)) (ftake 'mlist c) x))
((alike1 (add a b1 b2) c)
(ftake 'mexpt (sub 1 (add x y)) (mul -1 a)))
;; When all arguments are numbers and at least one is a binary64, use numerical evaluation.
((and (every #'mnump (list a b1 b2 c x y))
(some #'floatp (list a b1 b2 c x y))
(appell-f1-float a b1 b2 c x y)))
((great x y)
(ftake '%appell_f1 a b2 b1 c y x))
((alike1 y 1)
(mul (ftake '%hypergeometric (ftake 'mlist a b2) (ftake 'mlist c) 1)
(ftake '%hypergeometric (ftake 'mlist a b1) (ftake 'mlist (sub c b1)) x)))
((alike1 x (mul -1 x))
(ftake '%hypergeometric
(ftake (div (add a 1) 2) (div a 2) b1)
(ftake 'mlist (add (add c 1) 2) (div c 2))
(mul x x)))
;; http://dlmf.nist.gov/16.16.E1
((alike1 c (add b1 b2))
(mul (ftake 'mexpt (sub 1 y) (mul -1 a))
(ftake '%hypergeometric (ftake 'mlist a b1) (ftake 'mlist (add b1 b2))
(div (sub x y) (sub 1 y)))))
((and $use_appell_f1_integral_rep
(or (not (integerp c)) (< 0 c))
(or (not (integerp a)) (< 0 a))
(or (not (integerp (sub c a))) (< 0 (- c a))))
(integral-rep-appell-f1 a b1 b2 c x y))
(t (give-up))))
(defprop %appell_f1 (%appell_f1 simp) msimpind)
;; Integral representation
(defun integral-rep-appell-f1 (a b1 b2 c x y)
(let ((z ($gensym)))
(mul (div (ftake '%gamma c)
(mul (ftake '%gamma a) (ftake '%gamma (sub c a))))
(ftake '%integrate
(mul (ftake 'mexpt z (sub a 1))
(ftake 'mexpt (sub 1 z) (add c (mul -1 a) -1))
(ftake 'mexpt (sub 1 (mul x z)) (mul -1 b1))
(ftake 'mexpt (sub 1 (mul y z)) (mul -1 b2)))
z 0 1))))
;; Derivatives--we'll not attempt the derivatives with respect to the first four
;; arguments.
(defun diff-appell-f1-x (a b1 b2 c x y)
(if (eql c 0)
nil
(div (mul a b1 (ftake '%appell_f1 (add 1 a) (add 1 b1) b2 (add 1 c) x y)) c)))
(defun diff-appell-f1-y (a b1 b2 c x y)
(if (eql c 0)
nil
(div (mul a b2 (ftake '%appell_f1 (add 1 a) b1 (add 1 b2) (add 1 c) x y)) c)))
(defprop %appell_f1
((a b1 b2 c x y)
nil
nil
nil
nil
diff-appell-f1-x
diff-appell-f1-y)
grad)
;; Antiderivatives
(defun integrate-appell-f1-x (a b1 b2 c x y)
(div (mul (sub c 1) (ftake '%appell_f1 (sub a 1) (sub b1 1) b2 (sub c 1) x y))
(mul (sub a 1) (sub b1 1))))
(defun integrate-appell-f1-y (a b1 b2 c x y)
(div (mul (sub c 1) (ftake '%appell_f1 (sub a 1) b1 (sub b2 1) (sub c 1) x y))
(mul (sub a 1) (sub b2 1))))
(defprop %appell_f1
((a b1 b2 c x y)
nil
nil
nil
nil
integrate-appell-f1-x
integrate-appell-f1-y)
integral)
;; Use quad_qaws for floating point evaluation.
(defun appell-f1-float (a b1 b2 c x y)
(let ((s (gensym)) (g))
;; Convert inputs to float
(setq a ($float a)
b1 ($float b1)
b2 ($float b2)
c ($float c)
x ($float x)
y ($float y))
;; Check conditions and compute
(cond
((and (< x 1) (< y 1))
(setq g ($quad_qaws
(div 1
(mul (ftake 'mexpt (sub 1 (mul x s)) b1)
(ftake 'mexpt (sub 1 (mul y s)) b2)))
s 0 1 (- a 1) (- c (+ a 1)) 1))
(setq g (cdr g))
(if (eql 0 (fourth g)) (div (first g) (ftake '%beta a (- c a))) nil))
(t nil))))