2
2
;;; ** (c) Copyright 1979 Massachusetts Institute of Technology **
8
(eval-when (compile eval)
9
(declare-top (special fun w b l alglist $true $false n c l1 l2))
11
(DECLARE-TOP (SPECIAL VAR PAR ZEROSIGNTEST PRODUCTCASE
12
FLDEG FLGKUM CHECKCOEFSIGNLIST SERIESLIST
13
$EXPONENTIALIZE $BESTRIGLIM $RADEXPAND FAIL-SYM)
9
;; #+gcl (compile eval)
10
;; #-gcl (:compile-toplevel :execute)
11
;; (declare-top (special fun w b l alglist $true $false n c l1 l2)))
13
(declare-top (special $true $false))
15
(declare-top (special var *par*
17
$exponentialize $bestriglim $radexpand))
17
20
;; (eval-when (compile eval) (load '((dsk ell) macros >)) )
19
(declare-top (special fldeg flgkum listcmdiff checkcoefsignlist serieslist
21
(SETQ FLGKUM T FLDEG T FL1F1 T CHECKCOEFSIGNLIST NIL)
23
(declare-top (special $exponentialize $bestriglim $radexpand))
25
(setq fail-sym (gensym))
26
(defvar 3//2 '((rat simp) 3 2))
27
(defvar 1//2 '((rat simp) 1 2))
28
(defvar -1//2 '((rat simp) -1 2))
30
(eval-when (eval compile)
31
(defmacro fixp (x) `(typep ,x 'fixnum))
33
(setq FLGKUM T FLDEG T FL1F1 T CHECKCOEFSIGNLIST ()
34
;; $BESTRIGLIM 3. $RADEXPAND '$ALL
37
(DEFMACRO SIMP (X) `(SIMPLIFYA ,X ()))
39
(DEFMACRO SIMP-LIST (L) `(MAPCAR #'(LAMBDA (X) (SIMP X)) ,L))
41
; The macro MABS has been renamed to HYP-MABS in order to
42
; avoid conflict with the Maxima symbol MABS. The other
43
; M* macros defined here should probably be similarly renamed
44
; for consistency. jfa 03/27/2002
46
(DEFMACRO HYP-MABS (X) `(SIMP `((MABS) ,,X)))
48
(DEFMACRO MSQRT (X) `(M^T ,X 1//2))
50
(DEFMACRO MEXPT (X) `(M^T '$%E ,X))
52
(DEFMACRO MLOG (X) `(SIMP `((%LOG) ,,X)))
54
(DEFMACRO MSIN (X) `(SIMP `((%SIN) ,,X)))
56
(DEFMACRO MCOS (X) `(SIMP `((%COS) ,,X)))
58
(DEFMACRO MASIN (X) `(SIMP `((%ASIN) ,,X)))
60
(DEFMACRO MATAN (X) `(SIMP `((%ATAN) ,,X)))
62
(DEFMACRO MGAMMA (X) `(SIMP `((%GAMMA) ,,X)))
64
(DEFMACRO MBINOM (X Y) `(SIMP `((%BINOMIAL) ,,X ,,Y)))
66
(DEFMACRO MERF (X) `(SIMP `((%ERF) ,,X)))
68
(DEFMACRO =1//2 (X) `(ALIKE1 ,X 1//2))
70
(DEFMACRO =3//2 (X) `(ALIKE1 ,X 3//2))
72
(DEFMACRO =-1//2 (X) `(ALIKE1 ,X -1//2))
77
(L1 L2 ARG &aux ($bestriglim 3) ($radexpand '$all))
79
(setq var arg par arg)
80
(return (HGFSIMP-EXEC (CDR L1)(CDR L2) ARG))))
85
(setq l1 (copy-tree l1) l2 (copy-tree l2))
86
(PROG (RES $exponentialize)
89
(COND ((OR (NUMBERP RES)(NOT (ATOM RES)))
91
(RETURN (FPQFORM L1 L2 ARG))))
22
;; Why is this needed?
23
(setq checkcoefsignlist nil)
25
;; I (rtoy) don't know what the default should be. but $hgfred sets it
26
;; to 3. But we also need to define it because some of the specint
29
(defmvar $bestriglim 3)
31
(defmvar $prefer_whittaker nil)
35
#-gcl (:execute :compile-toplevel)
36
(defmacro fixp (x) `(typep ,x 'fixnum))
38
(setq checkcoefsignlist '())
39
;; $BESTRIGLIM 3. $RADEXPAND '$ALL
41
(defmacro simp (x) `(simplifya ,x ()))
43
(defmacro simp-list (l) `(mapcar #'(lambda (x) (simp x)) ,l))
45
;; The macro MABS has been renamed to HYP-MABS in order to
46
;; avoid conflict with the Maxima symbol MABS. The other
47
;; M* macros defined here should probably be similarly renamed
48
;; for consistency. jfa 03/27/2002
50
(defmacro hyp-mabs (x) `(simp `((mabs) ,,x)))
52
(defmacro msqrt (x) `(m^t ,x 1//2))
54
(defmacro mexpt (x) `(m^t '$%e ,x))
56
(defmacro mlog (x) `(simp `((%log) ,,x)))
58
(defmacro msin (x) `(simp `((%sin) ,,x)))
60
(defmacro mcos (x) `(simp `((%cos) ,,x)))
62
(defmacro masin (x) `(simp `((%asin) ,,x)))
64
(defmacro matan (x) `(simp `((%atan) ,,x)))
66
(defmacro mgamma (x) `(simp `((%gamma) ,,x)))
68
(defmacro mbinom (x y) `(simp `((%binomial) ,,x ,,y)))
70
(defmacro merf (x) `(simp `((%erf) ,,x)))
72
(defmacro =1//2 (x) `(alike1 ,x 1//2))
74
(defmacro =3//2 (x) `(alike1 ,x 3//2))
76
(defmacro =-1//2 (x) `(alike1 ,x -1//2))
79
(defun hyp-integerp (x)
80
;; In this file, maxima-integerp was used in many places. But it
81
;; seems that this code expects maxima-integerp to return T when it
82
;; is really an integer, not something that was declared an integer.
83
;; But I'm not really sure if this is true everywhere, but it is
84
;; true in some places.
86
;; Thus, we replace all calls to maxima-integerp with hyp-integerp,
87
;; which, for now, returns T only when the arg is an integer.
88
;; Should we do something more?
89
(and (maxima-integerp x) (integerp x)))
91
;; Main entry point for simplification of hypergeometric functions.
93
;; F(a1,a2,a3,...;b1,b2,b3;z)
95
;; L1 is a (maxima) list of an's, L2 is a (maxima) list of bn's.
96
(defun $hgfred (arg-l1 arg-l2 arg)
97
;; Do we really want $radexpand set to '$all? This is probably a
98
;; bad idea in general, but we'll leave this in for now until we can
99
;; verify find all of the code that does or does not need this and
100
;; until we can verify all of the test cases are correct.
101
(let (;;($radexpand '$all)
104
(hgfsimp-exec (cdr arg-l1) (cdr arg-l2) arg)))
107
(defun hgfsimp-exec (arg-l1 arg-l2 arg)
108
(let* ((l1 (copy-tree arg-l1))
109
(l2 (copy-tree arg-l2))
110
($exponentialize nil)
111
(res (hgfsimp l1 l2 arg)))
112
(if (or (numberp res) (not (atom res)))
114
(fpqform l1 l2 arg))))
116
(defun hgfsimp (arg-l1 arg-l2 var)
117
(prog (resimp listcmdiff)
118
(setq arg-l1 (macsimp arg-l1)
119
arg-l2 (macsimp arg-l2)
120
resimp (simpg arg-l1 arg-l2))
104
(COND ((NOT (EQ (CAR RESIMP) 'FAIL))(RETURN RESIMP)))
105
(COND ((SETQ LISTCMDIFF
106
(INTDIFFL1L2 (CADR RESIMP)
108
(return (splitpfq listcmdiff
111
(RETURN (DISPATCH-SPEC-SIMP (CADR RESIMP)
120
(T (APPEND (LIST (SIMPLIFYA (CAR L) NIL)) (CDR L)))))
126
(COND ((NULL (SETQ IL (zl-INTERSECTION L1 L2)))
127
(RETURN (SIMPG-EXEC L1 L2))))
128
(RETURN (SIMPG-EXEC (DEL IL L1)(DEL IL L2)))))
133
(COND ((NULL A) B)(T (DEL (CDR A) (ZL-DELETE (CAR A) B 1)))))
139
(COND ((ZEROP-IN-L L1)(RETURN 1)))
140
(COND ((SETQ N (hyp-NEGP-IN-L L1))
141
(RETURN (CREATE-POLY L1 L2 N))))
142
(COND ((OR (ZEROP-IN-L L2)(hyp-NEGP-IN-L L2))
144
(RETURN (APPEND (LIST 'FAIL)(LIST L1)(LIST L2)))))
122
(cond ((not (eq (and (consp resimp) (car resimp))
125
(cond ((setq listcmdiff
126
(intdiffl1l2 (cadr resimp)
128
(return (splitpfq listcmdiff
131
(return (dispatch-spec-simp (cadr resimp)
136
(mapcar #'(lambda (index)
137
(simplifya ($expand index) nil))
140
;; Simplify the parameters. If L1 and L2 have common elements, remove
141
;; them from both L1 and L2.
142
(defun simpg (arg-l1 arg-l2)
143
(let ((il (zl-intersection arg-l1 arg-l2)))
145
(simpg-exec arg-l1 arg-l2))
147
(simpg-exec (del il arg-l1)
153
(del (cdr a) (zl-delete (car a) b 1)))))
155
;; Handle the simple cases where the result is either a polynomial, or
156
;; is undefined because we divide by zero.
157
(defun simpg-exec (arg-l1 arg-l2)
159
(cond ((zerop-in-l arg-l1)
160
;; A zero in the first index means the series terminates
161
;; after the first term, so the result is always 1.
163
((setq n (hyp-negp-in-l arg-l1))
164
;; A negative integer in the first series means we have a
166
(create-poly arg-l1 arg-l2 n))
167
((or (zerop-in-l arg-l2)
168
(hyp-negp-in-l arg-l2))
169
;; A zero or negative number in the second index means we
170
;; eventually divide by zero, so we're undefined.
173
;; We failed so more complicated stuff needs to be done.
174
(append (list 'fail) (list arg-l1) (list arg-l2))))))
149
(COND ((NULL L1) NIL)(T (INTDIFF L1 L2))))
154
(SETQ L L2 A (CAR L1))
156
(COND ((NULL L)(RETURN (INTDIFFL1L2 (CDR L1) L2))))
157
(COND ((NNI (SETQ DIF (SUB A (CAR L))))
158
(RETURN (LIST A DIF))))
166
(COND ((AND (EQual LEN1 2)(EQual LEN2 1))
168
((AND (EQual LEN1 1)(EQual LEN2 1))
170
((AND (EQual LEN1 2)(ZEROP LEN2))
172
(T (CREATE-ANY-POLY L1 L2 (mul -1 N)))))
189
(MUL (POWER 2 (INV 2))(POWER VAR (INV 2))))
190
(COND ((EQUAL C (DIV 1 2))
192
(INV (FACTORIAL (ADD N N)))
193
(HERMPOL (ADD N N) FACT2)))))
194
(COND ((EQUAL C (DIV 3 2))
196
(INV (FACTORIAL (ADD N N 1)))
197
(HERMPOL (ADD N N 1) FACT2)))))
198
(RETURN (MUL (FACTORIAL N)
201
(LAGPOL N (SUB C 1) VAR)))))
204
(DEFUN HERMPOL(N ARG)(LIST '(MQAPPLY)(LIST '($%HE ARRAY) N) ARG))
205
(DEFUN LAGPOL(N A ARG)(LIST '(MQAPPLY)(LIST '($%L ARRAY) N A) ARG))
211
(SETQ A (CAR L1) B (CADR L1))
212
(COND ((EQUAL (SUB B A)(DIV -1 2))
213
(SETQ TEMP A A B B TEMP)))
214
(COND ((EQUAL (SUB B A)(DIV 1 2))
215
(SETQ X (POWER (DIV 2 (MUL -1 VAR))(INV 2)))
216
(RETURN (INTERHERMPOL N A B X))))
217
(SETQ X (MUL -1 (INV VAR)) N (MUL -1 N))
218
(RETURN (MUL (FACTORIAL N)
221
(LAGPOL N (ADD B N) X)))))
226
(SETQ FACT (POWER X (MUL -1 N)))
229
(RETURN (MUL FACT (HERMPOL N X)))))
231
(SETQ N (SUB 1 (ADD N N)))
232
(RETURN (MUL FACT (HERMPOL N X)))))))
238
(COND ((NOT (EQ (CAR L1) N))(setq l1 (REVERSE L1))))
239
(SETQ L (VFVP (DIV (ADD (CADR L1) N) 2)))
240
(SETQ V (CDR (ZL-assoc 'V L)))
242
(cond ((setq lgf (legpol (car l1)(cadr l1)(car l2)))
244
(COND ((EQUAL (SUB (CAR L2) V) '((RAT SIMP) 1 2))
247
(t (mul (factorial (* -1 n))
248
(inv (factf (mul 2 v)(* -1 n))))))
251
(SUB 1 (MUL 2 PAR)))))))
252
(RETURN (mul (factorial (* -1 n))
253
(inv (factf (add 1 v) (* -1 n)))
256
(SUB (MUL 2 V)(CAR L2))
257
(SUB 1 (MUL 2 PAR)))))))
262
(LIST '(MQAPPLY)(LIST '($%P ARRAY) N A B) X))
265
(DEFUN GEGENPOL(N V X)
266
(cond ((equal v 0) (tchebypol n x))
267
(t (LIST '(MQAPPLY)(LIST '($%C ARRAY) N V) X))))
268
(defun legenpol(n x)(list '(mqapply)(list '($%P array) n) x))
269
(defun tchebypol (n x)(list '(mqapply)(list '($%T array) n) x))
270
(DEFUN CREATE-ANY-POLY
272
(PROG(RESULT EXP PRODNUM PRODEN)
273
(SETQ RESULT 1 PRODNUM 1 PRODEN 1 EXP 1)
275
(COND ((ZEROP N) (RETURN RESULT)))
277
(MUL PRODNUM (MULL L1))
279
(MUL PRODEN (MULL L2)))
285
(INV (FACTORIAL EXP)))))
297
(DEFUN MULL(L)(COND ((NULL L) 1)(T (MUL (CAR L)(MULL (CDR L))))))
303
(T (APPEND (LIST (ADD (CAR L) 1))(INCR1 (CDR L))))))
306
(DEFUN DISPATCH-SPEC-SIMP
309
(SETQ LEN1 (LENGTH L1) LEN2 (LENGTH L2))
310
(COND ((AND (LESSP LEN1 2)(LESSP LEN2 2))
311
(RETURN (SIMP2>F<2 L1 L2 LEN1 LEN2))))
312
(COND ((AND (EQUAL LEN1 2)(EQUAL LEN2 1))
313
(RETURN (SIMP2F1 L1 L2))))
314
(RETURN (FPQFORM L1 L2 VAR))))
320
(COND ((AND (ZEROP LEN1)(ZEROP LEN2))
321
(RETURN (POWER '$%E VAR))))
322
(COND ((AND (ZEROP LEN1)(EQUAL LEN2 1))
323
(RETURN (BEStrig (CAR L2) VAR))))
324
(COND ((ZEROP LEN2)(RETURN (BINOM (CAR L1)))))
325
(RETURN (CONFL L1 L2 var))))
333
(setq res (mul (gm a) (power x (div (sub 1 a) 2))))
334
(COND ((AND (MAXIMA-INTEGERP (ADD A A))
335
(NUMBERP (SETQ N (SUB A (INV 2))))
336
(LESSP N $bestriglim))
338
(MEVAL (BESREDTRIG (- N 1)
345
(cond ((equal (checksigntm x) '$negative)
347
(BES (SUB A 1) (setq X (mul -1 x)) 'J)))))
348
(return (mul res (BES (SUB A 1) X 'I)))))
355
(LIST (COND ((EQ FLG 'J) '($%J ARRAY))
358
(MUL 2 (POWER X (INV 2)))))
365
(COND ((MINUSP N)(TRIGREDMINUS (MUL -1 (ADD1 N)) Z))
366
(T (TRIGREDPLUS N Z))))
371
(ADD (MUL (SIN% (SUB Z NPINV2))
373
(MUL (COS% (SUB Z NPINV2))
375
(MUL N '$%PI (INV 2))))
382
(SUB (MUL (COS% (ADD Z NPINV2))
384
(MUL (SIN% (ADD Z NPINV2))
386
(MUL N '$%PI (INV 2))))
390
(PROG(COUNT RESULT 2R N1)
391
(SETQ N1 ($ENTIER (DIV N 2)) COUNT 0 RESULT 1)
393
(COND ((EQ COUNT N1)(RETURN RESULT)))
400
(DIV (MUL (POWER -1 COUNT)
401
(FACTORIAL (ADD N 2R)))
403
(FACTORIAL (SUB N 2R))
404
(POWER (ADD Z Z) 2R)))))
409
(PROG(COUNT RESULT 2R+1 N1)
411
($ENTIER (DIV (SUB1 N) 2))
416
(COND ((EQual N1 -1)(RETURN 0)))
418
(COND ((EQ COUNT N1)(RETURN RESULT)))
425
(DIV (MUL (POWER -1 COUNT)
426
(FACTORIAL (ADD N 2R+1)))
427
(MUL (FACTORIAL 2R+1)
428
(FACTORIAL (SUB N 2R+1))
429
(POWER (ADD Z Z) 2R+1)))))
432
(DEFUN CTR(Z)(POWER (DIV 2 (MUL '$%PI Z))(INV 2)))
437
(COND ((NULL (SETQ D (CDR (ZL-REMPROP 'D (D*U X)))))
439
(COND ((EQ (ASKSIGN (INV D)) '$POSITIVE)
444
(DEFUN BINOM(A)(POWER (SUB 1 VAR) (MUL -1 A)))
450
(MUL (LIST '(MEXPT) '$%E VAR)
451
(confl (LIST (SUB (CAR L2)(CAR L1))) L2 (MUL -1 VAR))))
458
(COND ((ZEROP (CAR L)) T)(T (ZEROP-IN-L (CDR L)))))
459
(T (ZEROP-IN-L (CDR L)))))
465
((MAXIMA-INTEGERP (CAR L))
466
(COND ((MINUSP (CAR L)) (CAR L))
467
(T (hyp-NEGP-IN-L (CDR L)))))
468
(T (hyp-NEGP-IN-L (CDR L)))))
471
(DEFUN zl-INTERSECTION
473
(cond ((null l1) nil)
474
((zl-member (car l1) l2)
476
(zl-intersection (cdr l1)
477
(zl-delete (car l1) l2 1))))
478
(t (zl-intersection (cdr l1) l2))))
485
(COND ((AND (NULL L)(GREATERP COUNT 1))(RETURN T)))
486
(COND ((NULL L)(RETURN NIL)))
487
(COND ((MAXIMA-INTEGERP (CAR L))(SETQ COUNT (ADD1 COUNT))))
497
(COND ((AND (NULL L)(GREATERP COUNT 1))(RETURN T)))
498
(COND ((NULL L)(RETURN NIL)))
499
(COND ((EQ (CAAAR L) 'RAT)(SETQ COUNT (ADD1 COUNT))))
502
;2NUMP SHOULD BE ELIMINATED. IT IS NOT EFFICIENT TO USE ANYTHING ELSE BUT JUST CONVERTING TO RAT REPRESENTATION ALL 0.X ,X IN N. ESPECIALLY LATER WHEN WE CONVERT TO OMONIMA FOR TESTING TO FIND THE RIGHT FORMULA
510
(COND ((AND (NULL L)(GREATERP COUNT 1))(RETURN T)))
511
(COND ((NULL L)(RETURN NIL)))
512
(COND ((NUMBERP (CAR L))(SETQ COUNT (ADD1 COUNT))))
517
(DEFUN WHITFUN(K M VAR)(LIST '(MQAPPLY)(LIST '($%M ARRAY) K M) VAR))
522
(SETQ A (CAR L1) B (CADR L1) C (CAR L2))
523
(cond ((and (equal a 1)
526
(return (mul (inv (mul -1 var))
527
($log (add 1 (mul -1 var)))))))
528
(cond ((or (equal c (div 3 2))
530
(cond ((setq lgf (trig-log (list a b) (list c)))
534
(equal (sub a b) (div 1 2))
535
(equal (sub b a) (div 1 2)))
536
(cond ((setq lgf (hyp-cos a b c))(return lgf)))))
537
(cond ((and (maxima-integerp a)
538
(maxima-integerp b) (maxima-integerp c))
539
(return (simpr2f1 (list a b) (list c)))))
540
(cond ((and (maxima-integerp (add c (inv 2)))
541
(maxima-integerp (add a b)))
542
(return (step4 a b c))))
543
(cond ((maxima-integerp (add (sub a b) (inv 2)))
544
(cond ((setq lgf (step7 a b c))
546
(COND ((SETQ LGF (LEGFUN A B C))(RETURN LGF)))
547
(PRINT 'SIMP2F1-WILL-CONTINUE-IN)
548
(RETURN (FPQFORM L1 L2 VAR))))
178
(defun intdiffl1l2 (arg-l1 arg-l2)
182
(intdiff arg-l1 arg-l2))))
185
(defun intdiff (arg-l1 arg-l2)
187
(setq l arg-l2 a (car arg-l1))
189
(cond ((null l)(return (intdiffl1l2 (cdr arg-l1) arg-l2))))
190
(cond ((nni (setq dif (sub a (car l))))
191
(return (list a dif))))
195
;; For each element x in arg-l1 and y in arg-l2, compute d = x - y.
196
;; Find the smallest such non-negative integer d and return (list x
198
(defun intdiff (arg-l1 arg-l2)
200
;; Compute all possible differences between elements in arg-l1 and
201
;; arg-l2. Only save the ones where the difference is a positive
205
(let ((diff (sub x y)))
207
(push (list x diff) result)))))
208
;; Find the smallest one and return it.
209
(let ((min (first result)))
210
(dolist (x (rest result))
211
(when (lessp (second x) (second min))
216
;; Create the appropriate polynomial for the hypergeometric function.
217
(defun create-poly (arg-l1 arg-l2 n)
218
(let ((len1 (length arg-l1))
219
(len2 (length arg-l2)))
220
;; n is the smallest (in magnitude) negative integer in L1. To
221
;; make everything come out right, we need to make sure this value
222
;; is first in L1. This is ok, the definition of the
223
;; hypergeometric function does not depend on the order of values
225
(setf arg-l1 (cons n (remove n arg-l1 :count 1)))
226
(cond ((and (equal len1 2)
228
(2f1polys arg-l1 arg-l2 n))
235
(t (create-any-poly arg-l1 arg-l2 (mul -1 n))))))
238
(defun 1f1polys (arg-l2 n)
239
(let* ((c (car arg-l2))
241
(fact1 (mul (power 2 n)
244
;; For all of the polynomials here, I think it's ok to
245
;; replace sqrt(z^2) with z because when everything is
246
;; expanded out they evaluate to exactly the same thing. So
247
;; $radexpand $all is ok here.
248
(fact2 (let (($radexpand '$all))
249
(mul (power 2 (inv 2))
250
(power var (inv 2))))))
251
(cond ((alike1 c (div 1 2))
253
;; hermite(2*n,x) = (-1)^n*(2*n)!/n!*M(-n,1/2,x^2)
256
;; M(-n,1/2,x) = n!/(2*n)!*(-1)^n*hermite(2*n,sqrt(x))
258
;; But hermite(m,x) = 2^(m/2)*He(sqrt(2)*sqrt(x)), so
260
;; M(-n,1/2,x) = (-1)^n*n!*2^n/(2*n)!*He(2*n,sqrt(2)*sqrt(x))
262
(inv (factorial (add n n)))
263
(hermpol (add n n) fact2)))
264
((alike1 c (div 3 2))
266
;; hermite(2*n+1,x) = (-1)^n*(2*n+1)!/n!*M(-n,3/2,x^2)*2*x
269
;; M(-n,3/2,x) = n!/(2*n+1)!*(-1)^n*hermite(2*n+1,sqrt(x))/2/sqrt(x)
271
;; and in terms of He, we get
273
;; M(-n,3/2,x) = (-1)^n*n!*2^(n-1/2)/(2*n+1)!/sqrt(x)*He(2*n+1,sqrt(2)*sqrt(x))
275
(inv (power 2 (inv 2)))
276
(inv (factorial (add n n 1)))
277
(hermpol (add n n 1) fact2)
278
;; Similarly, $radexpand here is ok to convert sqrt(z^2) to z.
279
(let (($radexpand '$all))
280
(inv (power var (inv 2))))))
284
;; gen_laguerre(n,alpha,x) =
285
;; binomial(n+alpha,n)*hgfred([-n],[alpha+1],x);
287
;; Or hgfred([-n],[alpha],x) =
288
;; gen_laguerre(n,alpha-1,x)/binomial(n+alpha-1,n)
292
(lagpol n (sub c 1) var))))))
294
;; Hermite polynomial. Note: The Hermite polynomial used here is the
295
;; He polynomial, defined as (A&S 22.5.18, 22.5.19)
297
;; He(n,x) = 2^(-n/2)*H(n,x/sqrt(2))
301
;; H(n,x) = 2^(n/2)*He(x*sqrt(2))
303
;; We want to use H, as used in specfun, so we need to convert it.
305
(defun hermpol (n arg)
306
(let ((fact (inv (power 2 (div n 2))))
307
(x (mul arg (inv (power 2 (div 1 2))))))
308
(mul fact `(($hermite) ,n ,x))))
311
;; Generalized Laguerre polynomial
312
(defun lagpol (n a arg)
313
(if (and (numberp a) (zerop a))
314
`(($laguerre) ,n ,arg)
315
`(($gen_laguerre) ,n ,a, arg)))
318
(defun 2f0polys (arg-l1 n)
319
(let ((a (car arg-l1))
321
(when (alike1 (sub b a) (div -1 2))
323
(cond ((alike1 (sub b a) (div 1 2))
324
;; 2F0(-n,-n+1/2,z) or 2F0(-n-1/2,-n,z)
325
(interhermpol n a b var))
328
(let ((x (mul -1 (inv var)))
330
(mul (factorial order)
331
(inv (power x order))
332
(inv (power -1 order))
333
(lagpol order (mul -1 (add b order)) x)))))))
335
;; Compute 2F0(-n,-n+1/2;z) and 2F0(-n-1/2,-n;z) in terms of Hermite
338
;; Ok. I couldn't find any references giving expressions for this, so
339
;; here's a quick derivation.
341
;; 2F0(-n,-n+1/2;z) = sum(pochhammer(-n,k)*pochhammer(-n+1/2,k)*z^k/k!, k, 0, n)
343
;; It's easy to show pochhammer(-n,k) = (-1)^k*n!/(n-k)!
344
;; Also, it's straightforward but tedious to show that
345
;; pochhammer(-n+1/2,k) = (-1)^k*(2*n)!*(n-k)!/2^(2*k)/n!/(2*n-2*k)!
348
;; 2F0 = (2*n)!*sum(z^k/2^(2*k)/k!/(2*n-2*k)!)
350
;; Compare this to the expression for He(2*n,x) (A&S 22.3.11):
352
;; He(2*n,x) = (2*n)! * x^(2*n) * sum((-1)^k*x^(-2*k)/2^k/k!/(2*n-2*k)!)
356
;; 2F0(-n,-n+1/2;z) = y^n * He(2*n,y)
358
;; where y = sqrt(-2/x)
360
;; For 2F0(-n-1/2,-n;z) = sum(pochhammer(-n,k)*pochhammer(-n-1/2,k)*z^k/k!)
363
;; pochhammer(-n-1/2,k) = pochhammer(-(n+1)+1/2,k)
366
;; So 2F0 = (2*n+1)!*sum(z^k/z^(2*k)/k!/(2*n+1-2*k)!)
370
;; 2F0(-n-1/2,-n;z) = y^(2*n+1) * He(2*n+1,y)
373
(defun interhermpol (n a b x)
374
(let ((arg (power (div 2 (mul -1 x)) (inv 2)))
375
(order (cond ((alike1 a n)
378
(sub 1 (add n n))))))
379
;; 2F0(-n,-n+1/2;z) = y^(-2*n)*He(2*n,y)
380
;; 2F0(-n-1/2,-n;z) = y^(-(2*n+1))*He(2*n+1,y)
382
;; where y = sqrt(-2/var);
383
(mul (power arg (mul -1 order))
384
(hermpol order arg))))
386
;; F(a,b;c;z), where either a or b is a negative integer.
387
(defun 2f1polys (arg-l1 arg-l2 n)
389
;; Since F(a,b;c;z) = F(b,a;c;z), make sure L1 has the negative
390
;; integer first, so we have F(-n,d;c;z)
391
(cond ((not (eql (car arg-l1) n))
392
(setq arg-l1 (reverse arg-l1))))
394
(format t "l1 = ~A~%" arg-l1)
395
(format t "vfvp arg = ~A~%" (div (add (cadr arg-l1) n) 2))
396
(format t "var = ~A~%" var)
397
(format t "*par* = ~A~%" *par*)
399
(setq l (vfvp (div (add (cadr arg-l1) n) 2)))
401
;;(format t "l = ~A~%" l)
402
(setq v (cdr (zl-assoc 'v l)))
404
;; Assuming we have F(-n,b;c;z), then v is (b+n)/2.
406
;;(format t "v = ~A~%" v)
408
;; See if it can be a Legendre function.
409
(cond ((setq lgf (legpol (car arg-l1) (cadr arg-l1) (car arg-l2)))
412
(cond ((alike1 (sub (car arg-l2) v) '((rat simp) 1 2))
414
;; F(-n, n + 2*a; a + 1/2; x) = n!*gegen(n, a, 1-2*x)/pochhammer(2*a,n)
416
;; So v = a, and (car arg-l2) = a + 1/2.
419
(t (mul (factorial (* -1 n))
420
(inv (factf (mul 2 v) (* -1 n))))))
423
(sub 1 (mul 2 *par*)))))))
425
;; F(-n, n + a + 1 + b; a + 1; x) = n!*jacobi_p(n,a,b,1-2*x)/pochhammer(a+1,n);
427
(return (mul (factorial (* -1 n))
428
;; I (rlt) don't think this is right, based on
429
;; 15.4.6, because v doesn't have the right value.
430
#+nil(inv (factf (add 1 v) (* -1 n)))
431
;; Based on 15.4.6, we really want the arg-l2 arg
432
(inv (factf (car arg-l2) (* -1 n)))
434
(add (car arg-l2) -1)
435
(sub (mul 2 v) (car arg-l2))
436
(sub 1 (mul 2 *par*)))))))
442
(list '(mqapply)(list '($%p array) n a b) x))
446
(defun jacobpol (n a b x)
447
`(($jacobi_p) ,n ,a ,b ,x))
451
(defun gegenpol(n v x)
452
(cond ((equal v 0) (tchebypol n x))
453
(t (list '(mqapply)(list '($%c array) n v) x))))
455
;; Gegenbauer (Ultraspherical) polynomial. We use ultraspherical to
457
(defun gegenpol (n v x)
458
(cond ((equal v 0) (tchebypol n x))
459
(t `(($ultraspherical) ,n ,v ,x))))
461
;; Legendre polynomial
462
(defun legenpol (n x)
463
`(($legendre_p) ,n ,x))
465
;; Chebyshev polynomial
466
(defun tchebypol (n x)
467
`(($chebyshev_t) ,n ,x))
469
;; Expand the hypergeometric function as a polynomial. No real checks
470
;; are made to ensure the hypergeometric function reduces to a
472
(defun $hgfpoly (arg-l1 arg-l2 arg)
475
(n (hyp-negp-in-l (cdr arg-l1))))
476
(create-any-poly (cdr arg-l1) (cdr arg-l2) (- n))))
478
(defun create-any-poly (arg-l1 arg-l2 n)
479
(prog (result exp prodnum proden)
480
(setq result 1 prodnum 1 proden 1 exp 1)
482
(cond ((zerop n) (return result)))
483
(setq prodnum (mul prodnum (mull arg-l1))
484
proden (mul proden (mull arg-l2)))
490
(inv (factorial exp)))))
493
arg-l1 (incr1 arg-l1)
494
arg-l2 (incr1 arg-l2))
498
;; Compute the product of the elements of the list L.
501
(t (mul (car l) (mull (cdr l))))))
504
;; Add 1 to each element of the list L
507
(t (append (list (add (car l) 1)) (incr1 (cdr l))))))
511
(defun dispatch-spec-simp (arg-l1 arg-l2)
513
(setq len1 (length arg-l1) len2 (length arg-l2))
514
(cond ((and (lessp len1 2)
516
(return (simp2>f<2 arg-l1 arg-l2 len1 len2))))
517
(cond ((and (equal len1 2)
519
(return (simp2f1 arg-l1 arg-l2))))
520
(return (fpqform arg-l1 arg-l2 var))))
522
;; Figure out the orders of generalized hypergeometric function we
523
;; have and call the right simplifier.
524
(defun dispatch-spec-simp (arg-l1 arg-l2)
525
(let ((len1 (length arg-l1))
526
(len2 (length arg-l2)))
527
(cond ((and (lessp len1 2)
529
;; pFq where p and q < 2.
530
(simp2>f<2 arg-l1 arg-l2 len1 len2))
534
(simp2f1 arg-l1 arg-l2))
536
;; We don't have simplifiers for any other hypergeometric
538
(fpqform arg-l1 arg-l2 var)))))
542
(defun simp2>f<2 (l1 l2 len1 len2)
544
(cond ((and (zerop len1) (zerop len2))
545
(return (power '$%e var))))
546
(cond ((and (zerop len1) (equal len2 1))
547
(return (bestrig (car l2) var))))
548
(cond ((zerop len2) (return (binom (car l1)))))
549
(return (confl l1 l2 var))))
551
;; Handle the cases where the number of indices is less than 2.
552
(defun simp2>f<2 (arg-l1 arg-l2 len1 len2)
553
(cond ((and (zerop len1) (zerop len2))
554
;; hgfred([],[],z) = e^z
556
((and (zerop len1) (equal len2 1))
559
;; The hypergeometric series is then
561
;; 1+sum(z^k/k!/[b*(b+1)*...(b+k-1)], k, 1, inf)
563
;; = 1+sum(z^k/k!*gamma(b)/gamma(b+k), k, 1, inf)
564
;; = sum(z^k/k!*gamma(b)/gamma(b+k), k, 0, inf)
565
;; = gamma(b)*sum(z^k/k!/gamma(b+k), k, 0, inf)
567
;; Note that bessel_i(b,z) has the series
569
;; (z/2)^(b)*sum((z^2/4)^k/k!/gamma(b+k+1), k, 0, inf)
571
;; bessel_i(b-1,2*sqrt(z))
572
;; = (sqrt(z))^(b-1)*sum(z^k/k!/gamma(b+k),k,0,inf)
573
;; = z^((b-1)/2)*hgfred([],[b],z)/gamma(b)
575
;; So this hypergeometric series is a Bessel I function:
577
;; hgfred([],[b],z) = bessel_i(b-1,2*sqrt(z))*z^((1-b)/2)*gamma(b)
578
(bestrig (car arg-l2) var))
580
;; hgfred([a],[],z) = 1 + sum(binomial(a+k,k)*z^k)
582
(binom (car arg-l1)))
584
;; The general case of 1F1, the confluent hypergeomtric function.
585
(confl arg-l1 arg-l2 var))))
592
;; bessel_i(a-1,2*sqrt(x))*gamma(a)*x^((1-a)/2)
598
;; bessel_j(a-1,2*sqrt(x))*gamma(a)*x^((1-a)/2)
602
;; If a is half of an odd integer and small enough, the Bessel
603
;; functions are expanded in terms of trig or hyperbolic functions.
608
;; gamma(a)*x^((1-a)/2)
609
(setq res (mul (gm a) (power x (div (sub 1 a) 2))))
611
(cond ((and (hyp-integerp (add a a))
612
(numberp (setq n (sub a (inv 2))))
613
(lessp n $bestriglim))
614
;; This is totally broken. It's got an extra (-1)^foo
615
;; factor, so let's not use it at all for now. Use the
616
;; general forms below and let expand get the right
619
(meval (besredtrig (- n 1)
626
(cond ((equal (checksigntm x) '$negative)
627
;; Not sure this is right, but the call to bes has an
628
;; extra factor (-1)^(-(a-1)/2), so we cancel that out by
629
;; multiplying by (-1)^((a-1)/2).
631
(power -1 (div (sub a 1) 2))
632
(bes (sub a 1) (setq x (mul -1 x)) 'j)))))
633
(return (mul res (bes (sub a 1) x 'i)))))
636
;; I think it's ok to have $radexpand $all here so that sqrt(z^2) is converted to z.
637
(let* (($radexpand '$all)
638
(res (mul (gm a) (power x (div (sub 1 a) 2)))))
639
;; res = gamma(a)*x^((1-a)/2)
640
(if (equal (checksigntm x) '$negative)
641
;; Not sure this is right, but the call to bes has an
642
;; extra factor (-1)^(-(a-1)/2), so we cancel that out by
643
;; multiplying by (-1)^((a-1)/2).
645
(power -1 (div (sub a 1) 2))
646
(bes (sub a 1) (setq x (mul -1 x)) 'j))
647
(mul res (bes (sub a 1) x 'i)))))
650
(let ((fun (if (eq flg 'j) '%bessel_j '%bessel_i)))
651
`((,fun) ,a ,(mul 2 (power x (inv 2))))))
654
;; Compute bessel_j(n+1/2,z) in terms of trig functions.
656
;; See A&S 10.1.8 and 10.1.9.
658
;; Note that bessel.lisp has a different implementation of this.
659
;; Should we use that instead?
660
(defun besredtrig (n z)
662
(trigredminus (mul -1 (add1 n)) z))
663
(t (trigredplus n z))))
665
(defun trigredplus (n z)
668
(add (mul (sin% (sub z npinv2))
670
(mul (cos% (sub z npinv2))
672
(mul n '$%pi (inv 2))))
679
(sub (mul (cos% (add z npinv2))
681
(mul (sin% (add z npinv2))
683
(mul n '$%pi (inv 2))))
687
(prog(count result 2r n1)
688
(setq n1 ($entier (div n 2)) count 0 result 1)
690
(cond ((eq count n1)(return result)))
697
(div (mul (power -1 count)
698
(factorial (add n 2r)))
700
(factorial (sub n 2r))
701
(power (add z z) 2r)))))
704
;; Compute Q(n+1/2,z) in A&S 10.1.9.
705
(defun secondsum (n z)
706
(prog (count result 2r+1 n1)
708
($entier (div (sub1 n) 2))
712
(mul n (add 1 n) (inv (add z z))))
713
(cond ((equal n1 -1)(return 0)))
715
(cond ((eq count n1)(return result)))
722
(div (mul (power -1 count)
723
(factorial (add n 2r+1)))
724
(mul (factorial 2r+1)
725
(factorial (sub n 2r+1))
726
(power (add z z) 2r+1)))))
731
(power (div 2 (mul '$%pi z)) (inv 2)))
735
(cond ((null (setq d (cdr (zl-remprop 'd (d*u x)))))
737
(cond ((eq (asksign (inv d)) '$positive)
743
(power (sub 1 var) (mul -1 a)))
746
;; Kummer's transformation. A&S 13.1.27
748
;; M(a,b,z) = e^z*M(b-a,b,-z)
749
(defun kummer (arg-l1 arg-l2)
750
(mul (list '(mexpt) '$%e var)
751
(confl (list (sub (car arg-l2) (car arg-l1)))
752
arg-l2 (mul -1 var))))
755
;; Return non-NIL if any element of the list L is zero.
757
(defun zerop-in-l (l)
760
(cond ((zerop (car l)) t)(t (zerop-in-l (cdr l)))))
761
(t (zerop-in-l (cdr l)))))
763
(defun zerop-in-l (l)
765
(and (numberp x) (zerop x)))
769
(defun hyp-negp-in-l (l)
771
((hyp-integerp (car l))
772
(cond ((minusp (car l)) (car l))
773
(t (hyp-negp-in-l (cdr l)))))
774
(t (hyp-negp-in-l (cdr l)))))
776
;; If the list L contains a negative integer, return the most positive
777
;; of the negative integers. Otherwise return NIL.
778
(defun hyp-negp-in-l (l)
781
(when (and (numberp x) (minusp x))
783
(setf max-neg (max max-neg x))
787
;; Compute the intersection of L1 and L2, possibly destructively
788
;; modifying L2. Perserves duplications in L1.
789
(defun zl-intersection (arg-l1 arg-l2)
790
(cond ((null arg-l1) nil)
791
((zl-member (car arg-l1) arg-l2)
793
(zl-intersection (cdr arg-l1)
794
(zl-delete (car arg-l1) arg-l2 1))))
795
(t (zl-intersection (cdr arg-l1) arg-l2))))
797
;; Whittaker M function. A&S 13.1.32 defines it as
799
;; %m[k,u](z) = exp(-z/2)*z^(u+1/2)*M(1/2+u-k,1+2*u,z)
801
;; where M is the confluent hypergeometric function.
802
(defun whitfun (k m var)
803
(list '(mqapply) (list '($%m array) k m) var))
805
(defvar $trace2f1 nil
806
"Enables simple tracing of simp2f1 so you can see how it decides
807
what approach to use to simplify hypergeometric functions")
809
(defun simp2f1 (arg-l1 arg-l2)
811
(setq a (car arg-l1) b (cadr arg-l1) c (car arg-l2))
813
(format t "Tracing SIMP2F1~%")
814
(format t " Test F(1,1,2)...~%"))
815
(cond ((and (alike1 a 1)
818
;; F(1,1;2;z), A&S 15.1.3
821
(return (mul (inv (mul -1 var))
822
(mlog (add 1 (mul -1 var)))))))
824
(format t " Test c = 1/2 or c = 3/2...~%"))
825
(cond ((or (alike1 c (div 3 2))
826
(alike1 c (div 1 2)))
827
;; F(a,b; 3/2; z) or F(a,b;1/2;z)
828
(cond ((setq lgf (trig-log (list a b) (list c)))
830
(format t " Yes: trig-log~%"))
833
(format t " Test |a-b|=1/2...~%"))
835
(alike1 (sub a b) (div 1 2))
836
(alike1 (sub b a) (div 1 2)))
837
;; F(a,b;c;z) where |a-b|=1/2
838
(cond ((setq lgf (hyp-cos a b c))
840
(format t " Yes: hyp-cos~%"))
843
(format t " Test a,b are integers, c is a numerical integer...~%"))
844
(cond ((and (maxima-integerp a)
847
;; F(a,b;c;z) when a, and b are integers (or are declared
848
;; to be integers) and c is a integral number.
849
(setf lgf (simpr2f1 (list a b) (list c)))
852
(format t " Yes: simpr2f1~%"))
855
(format t " Test a+b and c+1/2 are numerical integers...~%"))
856
(cond ((and (hyp-integerp (add c (inv 2)))
857
(hyp-integerp (add a b)))
858
;; F(a,b;c;z) where a+b is an integer and c+1/2 is an
861
(format t " Yes: step4~%"))
862
(return (step4 a b c))))
864
(format t " Test a-b+1/2 is a numerical integer...~%"))
865
(cond ((hyp-integerp (add (sub a b) (inv 2)))
866
;; F(a,b;c,z) where a-b+1/2 is an integer
867
(cond ((setq lgf (step7 a b c))
870
(format t " Yes: step7~%"))
873
(when (and (hyp-integerp (add c 1//2))
874
(or (and (hyp-integerp (add a 1//2))
876
(and (hyp-integerp (add b 1//2))
879
(format t " Test for atanh: a+1/2, b, and c+1/2 are integers~%"))
880
(return (hyp-atanh a b c)))
882
(when (hyp-integerp (add c 1//2))
884
(format t " Test for atanh: c+1/2 is an integer~%"))
885
(cond ((and (hyp-integerp (add a 1//2))
888
(format t " atanh with integers a+1/2 and b ~%"))
889
(return (hyp-atanh a b c)))
890
((and (hyp-integerp (add b 1//2))
893
(format t " atanh with integers a and b+1/2 ~%"))
894
(return (hyp-atanh b a c)))))
897
(format t " Test for Legendre function...~%"))
898
(cond ((setq lgf (legfun a b c))
900
;; LEGFUN returned something interesting, so we're done.
902
(format t " Yes: case 1~%"))
904
;; LEGFUN didn't return anything, so try it with the args
905
;; reversed, since F(a,b;c;z) is F(b,a;c;z).
906
(setf lgf (legfun b a c))
909
(format t " Yes: case 2~%"))
913
(print 'simp2f1-will-continue-in)
915
(return (fpqform arg-l1 arg-l2 var))))
917
;; As best as I (rtoy) can tell, step7 is meant to handle an extension
918
;; of hyp-cos, which handles |a-b|=1/2 and either a+b-1/2 = c or
921
;; Based on the code, step7 wants a = s + m/n and c = 2*s+k/l. For
922
;; hyp-cos, we want c-2*a to be a integer. Or k/l-2*m/n is an
550
926
(defun step7 (a b c)
551
(prog (l m n k mn kl sym sym1 r)
557
(cond ((not (equal (mul sym 2) sym1))(return nil)))
558
(setq kl (cdras 'c l)
560
r (sub (add (inv 2) (cdras 'c l)) mn)
565
(cond ((equal (* 2 l) n)
566
(cond ((maxima-integerp (// (- k m) n))
567
(return (hyp-algv k l m n a b c))))))
568
(cond ((maxima-integerp (// k (* 2 l)))
569
(cond ((maxima-integerp (// m n))
570
(return (hyp-algv k l m n a b c)))
572
((maxima-integerp (// m n))
574
((maxima-integerp (/ (- (* k n) (* 2 l m)) (* 2 l n)))
575
(return (hyp-algv k l m n a b c))))
927
(prog (l m n k mn kl sym sym1 r)
928
;; Write a = sym + mn, c = sym1 + kl
934
;; We only handle the case where 2*sym = sym1.
935
(cond ((not (equal (mul sym 2) sym1))
937
(setq kl (cdras 'c l))
938
;; a-b+1/2 is an integer.
940
r (sub (add (inv 2) (cdras 'c l)) mn)
945
;; We have a = m*s+m/n, c = 2*m*s+k/l.
946
(cond ((equal (* 2 l) n)
947
(cond ((hyp-integerp (/ (- k m) n))
948
(return (hyp-algv k l m n a b c))))))
949
(cond ((hyp-integerp (/ k (* 2 l)))
950
(cond ((hyp-integerp (/ m n))
951
(return (hyp-algv k l m n a b c)))
953
((hyp-integerp (/ m n))
955
((hyp-integerp (/ (- (* k n) (* 2 l m)) (* 2 l n)))
956
(return (hyp-algv k l m n a b c))))
583
(cond ((maxima-integerp (setq x
588
(return (list x y))))
964
(cond ((hyp-integerp (setq x (// (+ y
968
(return (list x y))))
592
972
(defun hyp-algv (k l m n a b c)
595
(setq xy (getxy k l m n)
598
(cond ((< x 0)(go out)))
599
(cond ((< x y)(cond ((< (+ a-b x (inv 2)) 0)
600
(return (f88 x y a c fun)))
601
(t (return (f87 x y a c fun)))))
602
(t (cond ((< (+ a-b x (inv 2)) 0)
603
(return (f90 x y a c fun)))
604
(t (return (f89 x y a c fun))))))
607
(cond ((< (- (+ a-b (inv 2)) w) 0)
608
(return (f92 x y a c fun)))
609
(t (return (f91 x y a c fun))))))
975
(setq xy (getxy k l m n)
978
(cond ((< x 0)(go out)))
980
(cond ((< (add a-b x (inv 2)) 0)
981
(return (f88 x y a b c fun)))
982
(t (return (f87 x y a c fun)))))
984
(cond ((< (add a-b x (inv 2)) 0)
985
(return (f90 x y a c fun)))
986
(t (return (f89 x y a c fun))))))
989
(cond ((< (- (add a-b (inv 2)) w) 0)
990
(return (f92 x y a c fun)))
991
(t (return (f91 x y a c fun))))))
611
993
(defun f87 (x y a c fun )
613
(inv (mul (factf c y)
614
(factf (sub (add c y) (add a x)) (- x y))
615
(factf (sub (add c y) (add a x (inv 2)))
616
(sub (add a x (inv 2)) (add a (inv 2))))))
617
(power 'ell (sub 1 c))
618
(power (sub 1 'ell)(sub (add y c) (add a (inv 2))))
619
($diff (mul (power 'ell (add a x))
620
(power (sub 1 'ell)(mul -1 a))
621
($diff (mul (power 'ell (sub (add (inv 2) x) y))
622
($diff (mul (power 'ell (sub (add c y) 1))
995
(inv (mul (factf c y)
996
(factf (sub (add c y) (add a x)) (- x y))
997
(factf (sub (add c y) (add a x (inv 2)))
998
(sub (add a x (inv 2)) (add a (inv 2))))))
999
(power 'ell (sub 1 c))
1000
(power (sub 1 'ell)(sub (add y c) (add a (inv 2))))
1001
($diff (mul (power 'ell (add a x))
1002
(power (sub 1 'ell)(mul -1 a))
1003
($diff (mul (power 'ell (sub (add (inv 2) x) y))
1004
($diff (mul (power 'ell (sub (add c y) 1))
633
(defun f88 (x y a c fun )
635
(inv (mul (factf c y)
636
(factf (sub (add c y) (add a x)) (- x y))
637
(factf (add a (inv 2) x)
638
(sub b (sub x (sub a (inv 2)))))))
639
(power 'ell (sub 1 c))
640
(power (sub 1 'ell)(sub (add y c) (add a (inv 2))))
641
($diff (mul (power 'ell (add a x))
642
(power (sub 1 'ell)(mul -1 a))
643
($diff (mul (power 'ell (sub c (sub x (sub (inv 2) (mul a 2))))))
644
(power (sub 1 'ell) (sub (add a x b)(sub c y)))
645
($diff (mul (power 'ell (sub b 1 ))
1015
(defun f88 (x y a b c fun )
1017
(inv (mul (factf c y)
1018
(factf (sub (add c y) (add a x)) (- x y))
1019
(factf (add a (inv 2) x)
1020
(sub b (sub x (sub a (inv 2)))))))
1021
(power 'ell (sub 1 c))
1022
(power (sub 1 'ell)(sub (add y c) (add a (inv 2))))
1023
($diff (mul (power 'ell (add a x))
1024
(power (sub 1 'ell)(mul -1 a))
1025
($diff (mul (power 'ell (sub c (sub x (sub (inv 2) (mul a 2))))))
1026
(power (sub 1 'ell) (sub (add a x b)(sub c y)))
1027
($diff (mul (power 'ell (sub b 1 ))
648
'ell (sub b (sub a (sub (x (inv 2))))))
656
((LAMBDA (INL1P INL1BP INL2P)
657
(COND (INL2P (COND ((AND INL1P INL1BP)
658
(derivint (- (car l1) 1)
664
(INL1P (GEREDno2 (CADR L1)
667
(INL1BP (GEREDno2 (CAR L1)
671
(INL1P (COND (INL1BP 'D) (T 'C)))
672
((EQ (CAAAR L1) 'RAT)
673
(COND (INL1BP 'C) (T 'D)))
675
(MAXIMA-INTEGERP (CAR L1))
676
(MAXIMA-INTEGERP (CADR L1))
677
(MAXIMA-INTEGERP (CAR L2))))
680
(COND ((AND (GREATERP (CAR L2)(CAR L1))
681
(GREATERP (CAR L2)(CADR L1)))
682
(GEREDF (CAR L1)(CADR L1)(CAR L2)))
683
(T (GERED1 L1 L2 'HGFSIMP))))
686
(COND ((GREATERP C B)(GEREDF B A C))(T (GERED2 A B C))))
688
(n m l)(subst var 'psey
690
(factorial (+ n m l 1))
693
(inv (factorial (+ n m)))
694
(inv (factorial (+ m l)))
695
($diff (mul (power (sub 1 'psey) (+ m l))
696
($diff (mul (power 'psey -1)
698
($log (sub 1 'psey)))
709
(setq a1 (div (sub (add a b) (div 1 2)) 2))
710
(setq z1 (sub 1 var))
711
(setq a2 (mul c (inv 2)))
712
(cond ((equal (sub (add a b) (div 1 2)) c)
713
(return (mul (power 2 (sub (mul a1 2) 1))
714
(inv (power z1 (div 1 2)))
719
(sub 1 (mul 2 a1)))))))
720
(cond ((equal (add 1 (mul 2 a1)) c)
721
(return (mul (power 2 (sub c 1))
726
(mul -1 (sub c 1)))))))
731
(COND ((EQ (QUEST (SUB C B)) '$NEGATIVE)
732
(COND ((EQ (QUEST (SUB C A)) '$NEGATIVE)
733
(GERED1 (LIST A B)(LIST C) 'HGFSIMP))
735
((EQ (QUEST (SUB C A)) '$NEGATIVE)(GERED2 B A C))
736
(T (REST-DEGEN A B C))))
742
(COND ((NNI (SETQ M (SUB A 1)))
743
(RETURN (REST-DEGEN-1 A B C M))))
744
(COND ((NI B)(RETURN (REST-DEGEN-2 A B C))))
745
(COND ((AND (NNI (SETQ N (SUB C 1)))
746
(NNI (SETQ M (SUB (SUB A N) 1)))
747
(NNI (SETQ L (SUB B A)))
748
(EQ (SUB (SUB C A) B)
749
(MUL -1 (ADD M M N L 1))))
750
(RETURN (GERED1 (LIST A B)
753
(RETURN (hyp-DEG B A C))))
760
(NI (SUB (SUB C A) B))
761
(NNI (SUB (SUB C A) 1)))
762
(RETURN (DEG299 A B C))))
763
(COND ((AND (NNI (SETQ N (SUB (SUB C M) 2)))
764
(NNI (SETQ L (SUB B C)))
765
(EQUAL (SUB (SUB C A) B)
766
(MUL -1 (ADD L M 1))))
767
(RETURN (GERED1 (LIST A B)
770
(COND ((NNI (SETQ L (SUB (SUB B M) 1)))
771
(RETURN (REST-DEGEN-1A A B C M L))))
772
(RETURN (hyp-DEG B A C))))
778
(COND ((AND (NNI (SETQ N
779
(SUB (SUB (SUB C M) L) 2)))
780
(EQUAL (SUB N M)(SUB (SUB C A) B)))
781
(RETURN (DEG2913 A B C))))
782
(COND ((AND (EQUAL C (MUL -1 N))
783
(EQUAL (SUB (SUB C A) B)
784
(MUL -1 (ADD M M L N 2))))
785
(RETURN (DEG2918 A B C))))
786
(RETURN (hyp-DEG B A C))))
792
(COND ((AND (NI C)(NI (SUB (SUB C A) B)))
793
(RETURN (REST-DEGEN-2A A B C))))
794
(COND ((AND (NNI (SETQ M (SUB C 1)))
795
(NNI (SETQ L (SUB A C)))
796
(NI (SUB (SUB C A) B)))
797
(RETURN (DEG292 A B C))))
798
(RETURN (hyp-DEG B A C))))
804
(COND ((NNI (SUB A C))
805
(RETURN (GERED1 (LIST A B)
808
(COND ((NNI (SUB (SUB C A) 1))
809
(RETURN (DEG2917 A B C))))
810
(RETURN (hyp-DEG B A C))))
814
(COND ((NUMBERP A)(CHECKSIGNTM A))
815
((EQUAL (CAaR A) 'RAT)(CHECKSIGNTM A))
820
(DEFUN NNI(A)(COND ((MAXIMA-INTEGERP A)(NOT (MINUSP A)))))
823
(DEFUN NI(A)(NOT (MAXIMA-INTEGERP A)))
829
(COND (FLDEG (SETQ FLDEG NIL)
830
(RETURN (HGFSIMP (LIST A B)
834
(RETURN (FPQFORM (LIST A B)(LIST C) VAR))))
839
(MUL (POWER (MUL -1 VAR)(MUL -1 B))
840
(HGFSIMP (LIST (ADD B 1 (MUL -1 C)) B)
841
(LIST (ADD B 1 (MUL -1 A)))
847
(MUL (POWER VAR (SUB 1 C))
848
(POWER (SUB 1 VAR)(ADD C (MUL -1 A)(MUL -1 B)))
849
(HGFSIMP (LIST (SUB 1 A)(SUB 1 B))
856
(MUL (POWER VAR (SUB 1 C))
857
(HGFSIMP (LIST (ADD A 1 (MUL -1 C))
858
(ADD B 1 (MUL -1 C)))
865
(MUL (POWER (MUL -1 VAR)(MUL -1 A))
866
(HGFSIMP (LIST A (ADD A 1 (MUL -1 C)))
867
(LIST (ADD A 1 (MUL -1 B)))
873
(PROG(1-C A-B C-A-B INV2)
882
(COND ((EQUAL A-B INV2)
883
(RETURN (GERED1 (LIST A B)(LIST C) 'LEGF24))))
884
(COND ((EQUAL A-B (MUL -1 INV2))
885
(RETURN (LEGF24 (LIST A B)(LIST C) VAR))))
886
(COND ((EQUAL C-A-B INV2)
887
(RETURN (LEGF20 (LIST A B)(LIST C) VAR))))
888
(COND ((EQUAL C-A-B (MUL -1 INV2))
889
(RETURN (GERED1 (LIST A B)(LIST C) 'LEGF20))))
890
(COND ((EQUAL 1-C A-B)
891
(RETURN (LEGF16 (LIST A B)(LIST C) VAR))))
892
(COND ((EQUAL 1-C (MUL -1 A-B))
893
(RETURN (GERED1 (LIST A B)(LIST C) 'LEGF16))))
894
(COND ((EQUAL 1-C C-A-B)
895
(RETURN (GERED1 (LIST A B)(LIST C) 'LEGF14))))
896
(COND ((EQUAL 1-C (MUL -1 C-A-B))
897
(RETURN (LEGF14 (LIST A B)(LIST C) VAR))))
898
(COND ((EQUAL A-B (MUL -1 C-A-B))
899
(RETURN (LEGF36 (LIST A B)(LIST C) VAR))))
900
(COND ((OR (EQUAL 1-C INV2)
901
(EQUAL 1-C (MUL -1 INV2)))
902
(RETURN (LEGPOL A B C))))
903
(COND ((EQUAL A-B C-A-B)
904
(RETURN 'LEGENDRE-FUNCT-TO-BE-DISCOVERED)))
912
(SETQ B (CADR L1) C (CAR L2))
913
(SETQ M (SUB 1 C) N (MUL -1 (ADD B B M)))
914
(RETURN (MUL (LF N M)
917
(POWER (SUB 1 VAR) (INV 2))
931
(MUL -1 (ADD A A M)))
932
(RETURN (MUL (LF N M)
933
(POWER VAR (ADD N M))
936
(INV (POWER (SUB 1 VAR)
944
(SETQ A (CAR L1) C (CAR L2) M (SUB 1 C) N (MUL -1 A))
945
(RETURN (MUL (POWER 2 (MUL -1 N))
946
(POWER (SUB VAR 1)(DIV M -2))
948
(POWER (ADD VAR 1)(ADD (DIV M 2) N))
951
(DIV (ADD 1 VAR)(SUB 1 VAR))
958
(INV (POWER (SUB (POWER VAR 2) 1)(DIV M 2)))
959
(INV (GM (SUB 1 M)))))
965
(SETQ l (s+c (car l1))
966
a (cond ((eq (cdras 'c l) 0) (cdras 'f l))
967
(t (mul -1 (cdras 'f l))))
968
C (CAR L2) M (SUB 1 C)
970
(RETURN (MUL (POWER (ADD VAR 1)(DIV M 2))
971
(POWER (SUB VAR 1)(DIV M -2))
973
(LEGEN N M (SUB 1 (MUL 2 VAR)) '$P)))))
979
(SETQ A (CAR L1) B (CADR L1) N (SUB B 1) M (SUB B A))
980
(RETURN (MUL (POWER 2 N)
984
(ADD (DIV M 2)(MUL -1 N) -1))
985
(POWER (SUB VAR 1)(DIV M -2))
986
(INV (GM (ADD 2 N N)))
987
(POWER '$%E (MUL -1 '$%I M '$%PI))
988
(LEGEN N M (DIV (SUB 2 VAR) VAR) '$Q)))))
994
(LIST (COND ((EQ PQ '$Q) '($%Q ARRAY))
1004
(COND ((NOT (hyp-NEGP-IN-L (LIST A)))
1005
(RETURN 'FAIL-1-IN-C-1-CASE)))
1006
(SETQ L (VFVP (DIV (ADD B A) 2)))
1007
(SETQ V (CDR (ZL-ASSOC 'V L)))
1008
(COND ((AND (EQUAL V '((RAT SIMP) 1 2))(EQUAL C 1))
1009
(RETURN (LEGENPOL (MUL -1 A)
1010
(SUB 1 (MUL 2 VAR))))))
1011
(COND ((AND (EQUAL C '((RAT SIMP) 1 2))
1012
(EQUAL (SUB B A) '((RAT SIMP) 1 2)))
1013
(RETURN (MUL (FACTORIAL (MUL -1 A))
1015
(MULTAUG (INV 2) (MUL -1 A))
1016
(LEGENPOL (MUL -1 A)
1026
(COND ((ZEROP N) 1)(T (MUL A (MULTAUG (ADD A 1)(SUB1 N))))))
1031
(MUL (POWER (SUB 1 VAR)
1034
(MUL -1 (CADR L1))))
1036
(LIST (SUB (CAR L2) (CAR L1))
1037
(SUB (CAR L2) (CADR L1)))
1047
(MUL (POWER (SUB 1 VAR)(MUL -1 A))
1048
(HGFSIMP (LIST A (SUB C B))
1050
(DIV VAR (SUB VAR 1)))))
1055
(ADD (DIV (MUL (GM C)
1056
(GM (ADD C (MUL -1 A)(MUL -1 B)))
1057
(POWER VAR (MUL -1 A))
1058
(HGFSIMP (LIST A (ADD A 1 (MUL -1 C)))
1059
(LIST (ADD A B (MUL -1 C) 1))
1060
(SUB 1 (DIV 1 VAR))))
1061
(MUL (GM (SUB C A))(GM (SUB C B))))
1063
(GM (ADD A B (MUL -1 C)))
1065
(ADD C (MUL -1 A)(MUL -1 B)))
1066
(POWER VAR (SUB A C))
1067
(HGFSIMP (LIST (SUB C A)(SUB 1 A))
1030
'ell (sub b (sub a (sub x (inv 2)))))
1035
;; A new version of step7.
1036
(defun step7 (a b c)
1037
;; To get here, we know that a-b+1/2 is an integer. To make further
1038
;; progress, we want a+b-1/2-c to be an integer too.
1040
;; Let a-b+1/2 = p and a+b+1/2-c = q where p and q are (numerical)
1043
;; A&S 15.2.3 and 15.2.5 allows us to increase or decrease a. Then
1044
;; F(a,b;c;z) can be written in terms of F(a',b;c;z) where a' = a-p
1045
;; and a'-b = a-p-b = 1/2. Also, a'+b+1/2-c = a-p+b+1/2-c = q-p =
1046
;; r, which is also an integer.
1048
;; A&S 15.2.4 and 15.2.6 allows us to increase or decrese c. Thus,
1049
;; we can write F(a',b;c;z) in terms of F(a',b;c',z) where c' =
1050
;; c+r. Now a'-b=1/2 and a'+b+1/2-c' = a-p+b+1/2-c-r =
1051
;; a+b+1/2-c-p-r = q-p-(q-p)=0.
1053
;; Thus F(a',b;c';z) is exactly the form we want for hyp-cos. In
1054
;; fact, it's A&S 15.1.14: F(a,a+1/2,;1+2a;z) =
1055
;; 2^(2*a)*(1+sqrt(1-z))^(-2*a).
1056
(declare (special var))
1057
(let ((q (sub (add a b (inv 2))
1059
(unless (hyp-integerp q)
1060
;; Wrong form, so return NIL
1061
(return-from step7 nil))
1062
;; Since F(a,b;c;z) = F(b,a;c;z), we need to figure out which form
1063
;; to use. The criterion will be the fewest number of derivatives
1065
(let* ((p1 (add (sub a b) (inv 2)))
1067
(p2 (add (sub b a) (inv 2)))
1070
(format t "step 7:~%")
1071
(format t " q, p1, r1 = ~A ~A ~A~%" q p1 r1)
1072
(format t " p2, r2 = ~A ~A~%" p2 r2))
1073
(cond ((<= (+ (abs p1) (abs r1))
1074
(+ (abs p2) (abs r2)))
1077
(step7-core b a c))))))
1079
(defun step7-core (a b c)
1080
(let* ((p (add (sub a b) (inv 2)))
1081
(q (sub (add a b (inv 2))
1085
(c-prime (add 1 (mul 2 a-prime))))
1086
;; Ok, p and q are integers. We can compute something. There are
1087
;; four cases to handle depending on the sign of p and r.
1089
(let ((fun (hyp-cos a-prime (add a-prime 1//2) c-prime)))
1090
;; fun is F(a',a'+1/2;2*a'+1;z)
1092
(format t "step7-core~%")
1093
(format t " a,b,c = ~A ~A ~A~%" a b c)
1094
(format t " p,q,r = ~A ~A ~A~%" p q r)
1095
(format t " a', c' = ~A ~A~%" a-prime c-prime)
1096
(format t " F(a',a'+1/2; 1+2*a';z) =~%")
1097
(maxima-display fun))
1100
(step-7-pp a-prime b c-prime p r var fun))
1102
(step-7-pm a-prime b c-prime p r var fun))))
1105
(step-7-mp a-prime b c-prime p r var fun))
1107
(step-7-mm a-prime b c-prime p r var fun))))))))
1109
;; F(a,b;c;z) in terms of F(a',b;c';z)
1111
;; F(a'+p,b;c'-r;z) where p >= 0, r >= 0.
1112
(defun step-7-pp (a b c p r var fun)
1113
;; Apply A&S 15.2.4 and 15.2.3
1114
(let ((res (as-15.2.4 a b c r var fun)))
1115
(as-15.2.3 a b (sub c r) p var res)))
1120
;; F(a'+p,b;c'-r;z) = F(a'+p,b;c'+r';z)
1121
(defun step-7-pm (a b c p r var fun)
1122
;; Apply A&S 15.2.6 and 15.2.3
1123
(let ((res (as-15.2.6 a b c (- r) var fun)))
1124
(as-15.2.3 a b (sub c r) p var res)))
1129
;; F(a'+p,b;c'-r;z) = F(a'-p',b;c'-r;z)
1130
(defun step-7-mp (a b c p r var fun)
1131
;; Apply A&S 15.2.4 and 15.2.5
1132
(let ((res (as-15.2.4 a b c r var fun)))
1133
(as-15.2.5 a b (sub c r) (- p) var res)))
1137
;; Let p' = - p, r' = -r
1139
;; F(a'+p,b;c'-r;z) = F(a'-p',b;c'+r';z)
1140
(defun step-7-mm (a b c p r var fun)
1141
;; Apply A&S 15.2.6 and A&S 15.2.5
1142
(let ((res (as-15.2.6 a b c (- r) var fun)))
1143
(as-15.2.5 a b (sub c r) (- p) var res)))
1145
;; F(a,b;c;z) when a and b are integers (or declared to be integers)
1146
;; and c is an integral number.
1147
(defun simpr2f1 (arg-l1 arg-l2)
1148
(destructuring-bind (a b)
1150
(destructuring-bind (c)
1152
(let ((inl1p (hyp-integerp a))
1153
(inl1bp (hyp-integerp b))
1154
(inl2p (hyp-integerp c)))
1157
(cond ((and inl1p inl1bp)
1158
;; a, b, c are (numerical) integers
1161
;; a and c are integers
1164
;; b and c are integers.
1167
;; Can't really do anything else if c is not an integer.
1173
((eq (caaar arg-l1) 'rat)
1174
;; How do we ever get here?
1184
(cond ((and (greaterp (car arg-l2)(car arg-l1))
1185
(greaterp (car arg-l2)(cadr arg-l1)))
1186
(geredf (car arg-l1)(cadr arg-l1)(car arg-l2)))
1187
(t (gered1 arg-l1 arg-l2 #'hgfsimp))))
1189
(defun geredno2 (a b c)
1190
(cond ((greaterp c b) (geredf b a c))
1191
(t (gered2 a b c))))
1193
;; Consider F(1,1;2;z). A&S 15.1.3 says this is equal to -log(1-z)/z.
1195
;; Apply A&S 15.2.2:
1197
;; diff(F(1,1;2;z),z,ell) = poch(1,ell)*poch(1,ell)/poch(2,ell)*F(1+ell,1+ell;2+ell;z)
1201
;; diff((1-z)^(m+ell)*F(1+ell;1+ell;2+ell;z),z,m)
1202
;; = (-1)^m*poch(1+ell,m)*poch(1,m)/poch(2+ell,m)*(1-z)^ell*F(1+ell+m,1+ell;2+ell+m;z)
1206
;; diff((1-z)^ell*F(1+ell+m,1+ell;2+ell+m;z),z,n)
1207
;; = poch(1,n)*poch(1+m,n)/poch(2+ell+m,n)*(1-z)^(ell-n)*F(1+ell+m,1+ell;2+ell+m+n;z)
1209
;; The derivation above assumes that ell, m, and n are all
1210
;; non-negative integers. Thus, F(a,b;c;z), where a, b, and c are
1211
;; integers and a <= b <= c, can be written in terms of F(1,1;2;z).
1212
;; The result also holds for b <= a <= c, of course.
1214
;; Also note that the last two differentiations can be combined into
1215
;; one differention since the result of the first is in exactly the
1216
;; form required for the second. The code below does one
1219
;; So if a = 1+ell, b = 1+ell+m, and c = 2+ell+m+n, we have ell = a-1,
1220
;; m = b - a, and n = c - ell - m - 2 = c - b - 1.
1222
(defun derivint (a b c)
1230
(factorial (+ n m l 1))
1233
(inv (factorial (+ n m)))
1234
(inv (factorial (+ m l)))
1235
(power (sub 1 'psey) (sub n l))
1236
($diff (mul (power (sub 1 'psey) (+ m l))
1237
($diff (mul (power 'psey -1)
1239
(mlog (sub 1 'psey)))
1247
(defun hyp-cos (a b c)
1252
(setq a1 (div (sub (add a b) (div 1 2)) 2))
1253
(setq z1 (sub 1 var))
1254
(setq a2 (mul c (inv 2)))
1255
(cond ((equal (sub (add a b) (div 1 2)) c)
1258
;; 2^(2*a1 - 1)/sqrt(z1)*(1+sqrt(z1))^(1-2*a1)
1259
(return (mul (power 2 (sub (mul a1 2) 1))
1260
(inv (power z1 (div 1 2)))
1265
(sub 1 (mul 2 a1)))))))
1266
(cond ((equal (add 1 (mul 2 a1)) c)
1267
;; c = 1+2*a1 = a+b+1/2
1269
;; 2^(c-1)*(1+sqrt(z1))^(-(c-1))
1270
(return (mul (power 2 (sub c 1))
1275
(mul -1 (sub c 1)))))))
1278
(defun hyp-cos (a b c)
1279
(let ((a1 (div (sub (add a b) (div 1 2)) 2))
1283
(cond ((alike1 (sub (add a b)
1291
;; = 2^(2*a-1)*(1-z)^(-1/2)*(1+sqrt(1-z))^(1-2*a)
1293
;; But if 1-2*a is a negative integer, let's rationalize the answer to give
1296
;; = 2^(2*a-1)*(1-z)^(-1/2)*(1-sqrt(1-z))^(2*a-1)/z^(2*a-1)
1298
(format t " Case a+b-1/2=c~%"))
1299
(let ((2a-1 (sub (mul a1 2) 1)))
1300
(cond ((and (integerp 2a-1) (plusp 2a-1))
1301
;; 2^(2*a-1)*(1-z)^(-1/2)*(1-sqrt(1-z))^(2*a-1)/z^(2*a-1)
1303
(inv (power z1 1//2))
1304
(power (sub 1 (power z1 1//2)) 2a-1)
1305
(inv (power var 2a-1))))
1307
;; 2^(2*a-1)*(1-z)^(-1/2)*(1+sqrt(1-z))^(1-2*a)
1308
(mul (power 2 (sub (mul a1 2) 1))
1309
(inv (power z1 (div 1 2)))
1314
(sub 1 (mul 2 a1))))))))
1315
((alike1 (add 1 (mul 2 a1)) c)
1316
;; c = 1+2*a1 = a+b+1/2
1320
;; F(a,1/2+a;1+2*a;z) = 2^(2*a)*(1+sqrt(1-z))^(-2*a)
1322
;; But if 2*a is a positive integer, let's rationalize the answer to give
1324
;; F(a,1/2+a;1+2*a;z) = 2^(2*a)*(1-sqrt(1-z))^(2*a)/z^(2*a)
1326
(format t " Case c=1+2*a=a+b+1/2~%"))
1327
(let ((2a (sub c 1)))
1328
(cond ((and (integerp 2a) (plusp 2a))
1329
;; 2^(2*a)*(1-sqrt(1-z))^(2*a)/z^(2*a)
1331
(power (sub 1 (power z1 1//2))
1333
(power var (mul -1 2a))))
1335
;; 2^(2*a)*(1+sqrt(1-z))^(-2*a)
1337
(power (add 1 (power z1 1//2))
1338
(mul -1 2a))))))))))
1340
;; Is A a non-negative integer?
1342
(cond ((hyp-integerp a)
1346
;;; The following code doesn't appear to be used at all. Comment it all out for now.
1350
(cond ((eq (quest (sub c b)) '$negative)
1351
(cond ((eq (quest (sub c a)) '$negative)
1352
(gered1 (list a b)(list c) #'hgfsimp))
1353
(t (gered2 a b c))))
1354
((eq (quest (sub c a)) '$negative)(gered2 b a c))
1355
(t (rest-degen a b c))))
1361
(cond ((nni (setq m (sub a 1)))
1362
(return (rest-degen-1 a b c m))))
1363
(cond ((ni b)(return (rest-degen-2 a b c))))
1364
(cond ((and (nni (setq n (sub c 1)))
1365
(nni (setq m (sub (sub a n) 1)))
1366
(nni (setq l (sub b a)))
1367
(eq (sub (sub c a) b)
1368
(mul -1 (add m m n l 1))))
1369
(return (gered1 (list a b)
1372
(return (hyp-deg b a c))))
1379
(ni (sub (sub c a) b))
1380
(nni (sub (sub c a) 1)))
1381
(return (deg299 a b c))))
1382
(cond ((and (nni (setq n (sub (sub c m) 2)))
1383
(nni (setq l (sub b c)))
1384
(equal (sub (sub c a) b)
1385
(mul -1 (add l m 1))))
1386
(return (gered1 (list a b)
1389
(cond ((nni (setq l (sub (sub b m) 1)))
1390
(return (rest-degen-1a a b c m l))))
1391
(return (hyp-deg b a c))))
1394
(defun rest-degen-1a
1397
(cond ((and (nni (setq n
1398
(sub (sub (sub c m) l) 2)))
1399
(equal (sub n m)(sub (sub c a) b)))
1400
(return (deg2913 a b c))))
1401
(cond ((and (equal c (mul -1 n))
1402
(equal (sub (sub c a) b)
1403
(mul -1 (add m m l n 2))))
1404
(return (deg2918 a b c))))
1405
(return (hyp-deg b a c))))
1411
(cond ((and (ni c)(ni (sub (sub c a) b)))
1412
(return (rest-degen-2a a b c))))
1413
(cond ((and (nni (setq m (sub c 1)))
1414
(nni (setq l (sub a c)))
1415
(ni (sub (sub c a) b)))
1416
(return (deg292 a b c))))
1417
(return (hyp-deg b a c))))
1420
(defun rest-degen-2a
1423
(cond ((nni (sub a c))
1424
(return (gered1 (list a b)
1427
(cond ((nni (sub (sub c a) 1))
1428
(return (deg2917 a b c))))
1429
(return (hyp-deg b a c))))
1433
(cond ((numberp a)(checksigntm a))
1434
((equal (caar a) 'rat)(checksigntm a))
1437
(defun ni(a)(not (hyp-integerp a)))
1443
(cond (fldeg (setq fldeg nil)
1444
(return (hgfsimp (list a b)
1448
(return (fpqform (list a b)(list c) var))))
1453
(mul (power (mul -1 var)(mul -1 b))
1454
(hgfsimp (list (add b 1 (mul -1 c)) b)
1455
(list (add b 1 (mul -1 a)))
1461
(mul (power var (sub 1 c))
1462
(power (sub 1 var)(add c (mul -1 a)(mul -1 b)))
1463
(hgfsimp (list (sub 1 a)(sub 1 b))
1470
(mul (power var (sub 1 c))
1471
(hgfsimp (list (add a 1 (mul -1 c))
1472
(add b 1 (mul -1 c)))
1479
(mul (power (mul -1 var)(mul -1 a))
1480
(hgfsimp (list a (add a 1 (mul -1 c)))
1481
(list (add a 1 (mul -1 b)))
1487
;; Is F(a, b; c; z) is Legendre function?
1488
(defun legfun (a b c)
1489
(let ((1-c (sub 1 c))
1491
(c-a-b (sub (sub c a) b))
1493
(cond ((alike1 a-b inv2)
1496
(format t "Legendre a-b = 1/2~%"))
1497
(gered1 (list a b) (list c) #'legf24))
1498
((alike1 a-b (mul -1 inv2))
1501
;; For example F(a,a+1/2;c;x)
1503
(format t "Legendre a-b = -1/2~%"))
1504
(legf24 (list a b) (list c) var))
1505
((alike1 c-a-b inv2)
1508
;; For example F(a,b;a+b+1/2;z)
1510
(format t "Legendre c-a-b = 1/2~%"))
1511
(legf20 (list a b) (list c) var))
1512
((alike1 c-a-b (mul -1 inv2))
1515
(format t "Legendre c-a-b = -1/2~%"))
1516
(gered1 (list a b) (list c) #'legf20))
1520
(format t "Legendre 1-c = a-b~%"))
1521
(gered1 (list a b) (list c) #'legf16))
1522
((alike1 1-c (mul -1 a-b))
1525
(format t "Legendre 1-c = b-a~%"))
1526
(legf16 (list a b) (list c) var))
1530
(format t "Legendre 1-c = c-a-b~%"))
1531
(gered1 (list a b) (list c) #'legf14))
1532
((alike1 1-c (mul -1 c-a-b))
1535
;; For example F(a,1-a;c;x)
1537
(format t "Legendre 1-c = a+b-c~%"))
1538
(legf14 (list a b) (list c) var))
1539
((alike1 a-b (mul -1 c-a-b))
1542
;; For example F(a,b;2*b;z)
1544
(format t "Legendre a-b = a+b-c~%"))
1545
(legf36 (list a b) (list c) var))
1546
((or (alike1 1-c inv2)
1547
(alike1 1-c (mul -1 inv2)))
1548
;; 1-c = 1/2 or 1-c = -1/2
1550
;; For example F(a,b;1/2;z) or F(a,b;3/2;z)
1552
(format t "Legendre |1-c| = 1/2~%"))
1557
(format t "Legendre a-b = c-a-b~%"))
1558
'legendre-funct-to-be-discovered)
1563
;;; The following legf<n> functions correspond to formulas in Higher
1564
;;; Transcendental Functions. See the chapter on Legendre functions,
1565
;;; in particular the table on page 124ff,
1567
;; Handle the case c-a-b = 1/2
1571
;; P(n,m,z) = 2^m*(z^2-1)^(-m/2)/gamma(1-m)*F(1/2+n/2-m/2, -n/2-m/2; 1-m; 1-z^2)
1573
;; See also A&S 15.4.12 and 15.4.13.
1575
;; Let a = 1/2+n/2-m/2, b = -n/2-m/2, c = 1-m. Then, m = 1-c. And we
1576
;; have two equivalent expressions for n:
1578
;; n = c - 2*b - 1 or n = 2*a - c
1580
;; The code below chooses the first solution. A&S chooses second.
1582
;; F(a,b;c;w) = 2^(c-1)*gamma(c)*(-w)^((1-c)/2)*P(c-2*b-1,1-c,sqrt(1-w))
1586
(defun legf20 (arg-l1 arg-l2 var)
1587
(let* ((b (cadr arg-l1))
1590
(n (mul -1 (add b b m))))
1592
;; n = -(2*b+1-c) = c - 1 - 2*b
1594
(power 2 (mul -1 m))
1595
(power (mul -1 var) (div m 2))
1598
(power (sub 1 var) (inv 2))
1601
(defun legf20 (arg-l1 arg-l2 var)
1602
(let* (($radexpand nil)
1605
(a (sub (sub c b) (inv 2)))
1607
(n (mul -1 (add b b m))))
1609
;; n = -(2*b+1-c) = c - 1 - 2*b
1610
(cond ((and (eq (asksign var) '$positive)
1611
(eq (asksign (sub 1 var)) '$positive))
1614
;; F(a,b;a+b+1/2;x) = 2^(a+b-1/2)*gamma(a+b+1/2)*x^((1/2-a-b)/2)
1615
;; *assoc_legendre_p(a-b-1/2,1/2-a-b,sqrt(1-x))
1617
(mul (power 2 (add a b (inv -2)))
1618
(gm (add a b (inv 2)))
1620
(div (sub (inv 2) (add a b))
1624
(power (sub 1 var) (inv 2))
1627
(mul (power 2 (add a b (inv -2)))
1628
(gm (add a b (inv 2)))
1630
(div (sub (inv 2) (add a b))
1634
(power (sub 1 var) (inv 2))
1637
;; Handle the case a-b = -1/2.
1641
;; P(n,m,z) = 2^m*(z^2-1)^(-m/2)*z^(n+m)/gamma(1-m)*F(-n/2-m/2,1/2-n/2-m/2;1-m;1-1/z^2)
1643
;; See also A&S 15.4.10 and 15.4.11.
1645
;; Let a = -n/2-m/2, b = 1/2-n/2-m/2, c = 1-m. Then m = 1-c. Again,
1646
;; we have 2 possible (equivalent) values for n:
1648
;; n = -(2*a + 1 - c) or n = c-2*b
1650
;; The code below chooses the first solution.
1652
;; F(a,b;c;w) = 2^(c-1)*w^(1/2-c/2)*(1-w)^(c/2-a-1/2)*P(c-2*a-1,1-c,1/sqrt(1-w))
1654
;; F(a,b;c;w) = 2^(c-1)*w^(1/2-c/2)*(1-w)^(c/2-b)*P(c-2*b,1-c,sqrt(1-w))
1656
;; Is there a mistake in 15.4.10 and 15.4.11?
1659
(defun legf24 (arg-l1 arg-l2 var)
1660
(let* ((a (car arg-l1))
1663
(n (mul -1 (add a a m)))
1664
(z (inv (power (sub 1 var) (inv 2)))))
1665
(mul (inv (power 2 m))
1666
(power (sub (power z 2) 1)
1668
(power z (mul -1 (add n m)))
1675
(defun legf24 (arg-l1 arg-l2 var)
1676
(let* (($radexpand nil)
1680
(n (mul -1 (add a a m)))
1681
(z (inv (power (sub 1 var) (inv 2)))))
1682
;; A&S 15.4.10, 15.4.11
1683
(cond ((eq (asksign var) '$negative)
1686
;; F(a,a+1/2;c;x) = 2^(c-1)*gamma(c)(-x)^(1/2-c/2)*(1-x)^(c/2-a-1/2)
1687
;; *assoc_legendre_p(2*a-c,1-c,1/sqrt(1-x))
1688
(mul (inv (power 2 m))
1690
(power (mul -1 var) (div m 2))
1691
(power (sub 1 var) (sub (div m -2) a))
1697
(mul (inv (power 2 m))
1699
(power var (div m 2))
1700
(power (sub 1 var) (sub (div m -2) a))
1712
;; P(n,m,z) = 2^(-n)*(z+1)^(m/2+n)*(z-1)^(-m/2)/gamma(1-m)*F(-n,-n-m;1-m;(z-1)/(z+1))
1714
;; See also A&S 15.4.14 and 15.4.15.
1716
;; Let a = -n, b = -n-m, c = 1-m. Then m = 1-c. We have 2 solutions
1719
;; n = -a or n = c-b-1.
1721
;; The code below chooses the first solution.
1723
;; F(a,b;c;w) = gamma(c)*w^(1/2-c/2)*(1-w)^(-a)*P(-a,1-c,(1+w)/(1-w));
1725
;; FIXME: We don't correctly handle the branch cut here!
1727
(defun legf16 (arg-l1 arg-l2 var)
1728
(let* ((a (car arg-l1))
1735
(power (sub z 1) (div m 2))
1737
(inv (power (add z 1) (add (div m 2) n)))
1743
(defun legf16 (arg-l1 arg-l2 var)
1744
(let* (($radexpand nil)
1750
;; m = b - a, so b = a + m
1756
(format t "a, c = ~A ~A~%" a c)
1757
(format t "b = ~A~%" b))
1758
;; A&S 15.4.14, 15.4.15
1759
(cond ((eq (asksign var) '$negative)
1762
;; F(a,b;a-b+1,x) = gamma(a-b+1)*(1-x)^(-b)*(-x)^(b/2-a/2)
1763
;; * assoc_legendre_p(-b,b-a,(1+x)/(1-x))
1767
(power (sub 1 var) (mul -1 b))
1768
(power (mul -1 var) (div m 2))
1775
(power (sub 1 var) (mul -1 b))
1776
(power var (div m 2))
1783
;; Handle the case 1-c = a+b-c.
1785
;; See, for example, A&S 8.1.2 (which
1786
;; might have a bug?) or
1787
;; http://functions.wolfram.com/HypergeometricFunctions/LegendreP2General/26/01/02/
1791
;; P(n,m,z) = (z+1)^(m/2)*(z-1)^(-m/2)/gamma(1-m)*F(-n,1+n;1-m;(1-z)/2)
1793
;; See also A&S 8.1.2, 15.4.16, 15.4.17
1795
;; Let a=-n, b = 1+n, c = 1-m. Then m = 1-c and n has 2 solutions:
1797
;; n = -a or n = b - 1.
1799
;; The code belows chooses the first solution.
1801
;; F(a,b;c;w) = gamma(c)*(-w)^(1/2-c/2)*(1-w)^(c/2-1/2)*P(-a,1-c,1-2*w)
1803
(defun legf14 (arg-l1 arg-l2 var)
1804
(let* ((a (car arg-l1))
1808
(z (sub 1 (mul 2 var))))
1809
(mul (power (add z 1) (div m -2))
1810
(power (sub z 1) (div m 2))
1812
(legen n m (sub 1 (mul 2 var)) '$p))))
1814
(defun legf14 (arg-l1 arg-l2 var)
1815
;; Set $radexpand to NIL, because we don't want (-z)^x getting
1816
;; expanded to (-1)^x*z^x because that's wrong for this.
1817
(let* (($radexpand nil)
1822
(z (sub 1 (mul 2 var))))
1823
;; A&S 15.4.16, 15.4.17
1824
(cond ((and (eq (asksign var) '$positive)
1825
(eq (asksign (sub 1 var)) '$positive))
1828
;; F(a,1-a;c;x) = gamma(c)*x^(1/2-c/2)*(1-x)^(c/2-1/2)*
1829
;; assoc_legendre_p(-a,1-c,1-2*x)
1833
(power var (div (sub 1 c) 2))
1834
(power (sub 1 var) (div (sub c 1) 2))
1839
;; F(a,1-a;c;z) = gamma(c)*(-z)^(1/2-c/2)*(1-z)^(c/2-1/2)*
1840
;; assoc_legendre_p(-a,1-c,1-2*z)
1842
(power (mul -1 var) (div (sub 1 c) 2))
1843
(power (sub 1 var) (div (sub c 1) 2))
1844
(legen n m z '$p))))))
1846
;; Handle a-b = a+b-c
1850
;; exp(-%i*m*%pi)*Q(n,m,z) =
1851
;; 2^n*gamma(1+n)*gamma(1+n+m)*(z+1)^(m/2-n-1)*(z-1)^(-m/2)/gamma(2+2*n)
1852
;; * hgfred([1+n-m,1+n],[2+2*n],2/(1+z))
1854
;; Let a = 1+n-m, b = 1+n, c = 2+2*n. then n = b-1 and m = b - a.
1855
;; (There are other solutions.)
1857
;; F(a,b;c;z) = 2*gamma(2*b)/gamma(b)/gamma(2*b-a)*w^(-b)*(1-w)^((b-a)/2)
1858
;; *Q(b-1,b-a,2/w-1)*exp(-%i*%pi*(b-a))
1860
(defun legf36 (arg-l1 arg-l2 var)
1861
(declare (ignore arg-l2))
1862
(let* ((a (car arg-l1))
1866
;;z (div (sub 2 var) var)
1867
(z (sub (div 2 var) 1)))
1868
(mul (inv (power 2 n))
1869
(inv (gm (add 1 n)))
1870
(inv (gm (add 1 n m)))
1871
(inv (power (add z 1)
1875
(inv (power (sub z 1) (div m -2)))
1877
(power '$%e (mul -1 '$%i m '$%pi))
1878
(legen n m z '$q))))
1880
(defun legen (n m x pq)
1881
;; A&S 8.2.1: P(-n-1,m,z) = P(n,m,z)
1883
;; Currently only applied if n is a number. (Should this be
1884
;; extended to any expression? We'll have to ask the user for the
1885
;; sign if we can' figure it out ourselves. Should we?)
1886
(let ((n (if (and (mnump n)
1887
(eq (checksigntm n) '$negative))
1891
`((,(if (eq pq '$q) '$legendre_q '$legendre_p)) ,n ,x))
1893
`((,(if (eq pq '$q) '$assoc_legendre_q '$assoc_legendre_p))
1897
(defun legpol (a b c)
1898
;; Why do we insist that a be a negative (numerical) integer?
1899
(when (not (hyp-negp-in-l (list a)))
1900
(print 'fail-1-in-c-1-case)
1901
(return-from legpol nil))
1902
(let* ((l (vfvp (div (add b a) 2)))
1903
(v (cdr (zl-assoc 'v l))))
1906
((and (alike1 v '((rat simp) 1 2))
1909
;; P(n,x) = F(-n,n+1;1;(1-x)/2)
1910
(legenpol (mul -1 a)
1911
(sub 1 (mul 2 var))))
1913
((and (alike1 c '((rat simp) 1 2))
1914
(alike1 (add b a) '((rat simp) 1 2)))
1916
;; P(2*n,x) = (-1)^n*(2*n)!/2^(2*n)/(n!)^2*F(-n,n+1/2;1/2;x^2)
1918
;; F(-n,n+1/2;1/2;x^2) = P(2*n,x)*(-1)^n*(n!)^2/(2*n)!*2^(2*n)
1920
(let ((n (mul -1 a)))
1922
(power (factorial n) 2)
1923
(inv (factorial (mul 2 n)))
1926
(power var (div 1 2))))))
1928
((and (alike1 c '((rat simp) 3 2))
1929
(alike1 (add b a) '((rat simp) 3 2)))
1931
;; P(2*n+1,x) = (-1)^n*(2*n+1)!/2^(2*n)/(n!)^2*F(-n,n+3/2;3/2;x^2)*x
1933
;; F(-n,n+3/2;3/2;x^2) = P(2*n+1,x)*(-1)^n*(n!)^2/(2*n+1)!*2^(2*n)/x
1935
(let ((n (mul -1 a)))
1937
(power (factorial n) 2)
1938
(inv (factorial (add 1 (mul 2 n))))
1940
(legenpol (add 1 (mul 2 n))
1941
(power var (div 1 2)))
1942
(inv (power var (div 1 2))))))
1944
((and (zerp (sub b a))
1945
(zerp (sub c (add a b))))
1947
;; P(n,x) = binomial(2*n,n)*((x-1)/2)^n*F(-n,-n;-2*n;2/(1-x))
1949
;; F(-n,-n;-2*n;x) = P(n,1-2/x)/binomial(2*n,n)(-1/x)^(-n)
1950
(mul (power (factorial (mul -1 a)) 2)
1951
(inv (factorial (mul -2 a)))
1952
(power (mul -1 var) (mul -1 a))
1953
(legenpol (mul -1 a)
1954
(add 1 (div -2 var)))))
1955
((and (alike1 (sub a b) '((rat simp) 1 2))
1956
(alike1 (sub c (mul 2 b)) '((rat simp) 1 2)))
1958
;; P(n,x) = binomial(2*n,n)*(x/2)^n*F(-n/2,(1-n)/2;1/2-n;1/x^2)
1960
;; F(-n/2,(1-n)/2;1/2-n,1/x^2) = P(n,x)/binomial(2*n,n)*(x/2)^(-n)
1961
(mul (power (factorial (mul -2 b)) 2)
1962
(inv (factorial (mul -4 b)))
1963
(power (mul 2 (power var (div 1 2))) (mul -2 b))
1964
(legenpol (mul -2 b)
1965
(power var (div -1 2)))))
1966
((and (alike1 (sub b a) '((rat simp) 1 2))
1967
(alike1 (sub c (mul 2 a)) '((rat simp) 1 2)))
1969
;; P(n,x) = binomial(2*n,n)*(x/2)^n*F(-n/2,(1-n)/2;1/2-n;1/x^2)
1971
;; F(-n/2,(1-n)/2;1/2-n,1/x^2) = P(n,x)/binomial(2*n,n)*(x/2)^(-n)
1972
(mul (power (factorial (mul -2 a)) 2)
1973
(inv (factorial (mul -4 a)))
1974
(power (mul 2 (power var (div 1 2))) (mul -2 a))
1975
(legenpol (mul -2 a)
1976
(power var (div -1 2)))))
1984
;; F(a,b;c;z) = (1-z)^(c-a-b)*F(c-a,c-b;c;z)
1985
(defun gered1 (arg-l1 arg-l2 simpflg)
1986
(destructuring-bind (a b)
1988
(destructuring-bind (c)
1990
(mul (power (sub 1 var)
2001
;; F(a,b;c;z) = (1-z)^(-a)*F(a,c-b;c;z/(z-1))
2002
(defun gered2 (a b c)
2003
(mul (power (sub 1 var) (mul -1 a))
2004
(hgfsimp (list a (sub c b))
2006
(div var (sub var 1)))))
2010
;; F(a,b;c;z) = A*z^(-a)*F(a,a-c+1;a+b-c+1;1-1/z)
2011
;; + B*(1-z)^(c-a-b)*z^(a-c)*F(c-a,1-a;c-a-b+1,1-1/z)
2013
;; where A = gamma(c)*gamma(c-a-b)/gamma(c-a)/gamma(c-b)
2014
;; B = gamma(c)*gamma(a+b-c)/gamma(a)/gamma(b)
2016
(defun geredf (a b c)
2017
(add (div (mul (gm c)
2018
(gm (add c (mul -1 a)(mul -1 b)))
2019
(power var (mul -1 a))
2020
($hgfred `((mlist) ,a ,(add a 1 (mul -1 c)))
2021
`((mlist) ,(add a b (mul -1 c) 1))
2022
(sub 1 (div 1 var))))
2023
(mul (gm (sub c a))(gm (sub c b))))
2025
(gm (add a b (mul -1 c)))
2027
(add c (mul -1 a)(mul -1 b)))
2028
(power var (sub a c))
2029
($hgfred `((mlist) ,(sub c a) ,(sub 1 a))
1072
(SUB 1 (DIV 1 VAR))))
1073
(MUL (GM A)(GM B)))))
1079
(COND ((EQUAL (SIMPLIFYA (CAR L2) NIL) '((RAT SIMP) 3 2))
1081
((EQUAL (SIMPLIFYA (CAR L2) NIL) '((RAT SIMP) 1 2))
1088
(COND ((AND (OR (equal (car l1) 1) (equal (cadr l1) 1))
1089
(OR (equal (car l1) (div 1 2))
1090
(equal (cadr l1) (div 1 2))))
1091
(TRIG-LOG-3-EXEC L1 L2))
1092
((and (equal (car l1) (cadr l1))
1093
(or (equal 1 (car l1))
1094
(equal (div 1 2) (car l1))))
1095
(trig-log-3a-exec l1 l2))
1096
((or(equal (add (car l1) (cadr l1)) 1)
1097
(equal (add (car l1) (cadr l1)) 2))
1099
((or (equal (sub (car l1) (cadr l1)) (div 1 2))
1100
(equal (sub (cadr l1) (car l1)) (div 1 2)))
1107
(return (mul (inv (setq z (power var (div 1 2))))
1114
(sub (power (add 1 z) a)
1115
(power (sub 1 z) a))))))
1119
(setq a (car l1) b (cadr l1) c (car l2))
1120
(cond ((equal (add a b) 1)
1121
(return (mul (inv (mul (mul -1 (sub a b))
1122
($sin ($asin ($sqrt var)))))
1125
($asin ($sqrt var)))))))
1127
(return (mul ($sin (mul (setq z1
1139
;Generates atan if arg positive else log
1140
(defun trig-log-3-exec
1143
(cond ((equal (checksigntm var) '$positive)
1144
(return (mul (power (setq z
1150
($log (div (add 1 z)
1152
((equal (checksigntm var) '$negative)
1153
(return (mul (power (setq z
1163
(prog (a b c z1 $exponentialize)
1165
(setq a (car l1) b (cadr l1) c (car l2))
1166
(cond ((equal (add a b) 0)
1167
(cond ((equal (checksigntm var) '$positive)
1168
(return ($cos (mul (mul 2 a)
1171
(t (return (div (add (power (add (setq
1185
((equal (add a b) 1)
1186
(return (mul (inv ($cos (setq z1
1190
($cos (mul z1 (sub a b))))))
1191
((or (equal (sub a b) (inv 2))
1192
(equal (sub a b) (inv -2)))
1193
(return (add (div (power (add 1
1202
(div (power (sub 1 z1)
1210
(DEFUN TRIG-LOG-1 (A B) ;; 2F1's with C = 1/2
1211
(LET (X Z $EXPONENTIALIZE) ;; 15.1.17, 11, 18, 12, 9, and 19
1212
(setq a (car l1) b (cadr l1))
1213
(COND ((=0 (M+T A B))
1214
(COND ((EQUAL (CHECKSIGNTM VAR) '$POSITIVE)
1215
(MCOS (M*T 2. A (MASIN (MSQRT VAR)))))
1216
((EQUAL (CHECKSIGNTM VAR) '$NEGATIVE)
1218
(M+T (M^T (M+T (SETQ X (MSQRT (M-T 1. VAR)))
1219
(SETQ Z (MSQRT (M-T VAR))))
1220
(SETQ B (M*T 2. B)))
1221
(M^T (M-T X Z) B))))
1223
((EQUAL (M+T A B) 1.)
1224
(COND ((EQUAL (CHECKSIGNTM VAR) '$POSITIVE)
1225
(M//T (MCOS (M*T (M-T A B) (SETQ Z (MASIN (MSQRT VAR)))))
1227
((EQUAL (CHECKSIGNTM VAR) '$NEGATIVE)
1228
(M*T 1//2 (M//T (SETQ X (MSQRT (M-T 1. VAR))))
1229
(M+T (M^T (M+T X (SETQ Z (MSQRT (M-T VAR))))
1231
(M^T (M-T X Z) B))))
1233
((=1//2 (HYP-MABS (M-T B A)))
1234
(COND ((EQUAL (CHECKSIGNTM VAR) '$POSITIVE)
1236
(M+T (M^T (M1+T (SETQ Z (MSQRT VAR)))
1237
(SETQ B (M-T 1//2 (M+T A B))))
1238
(M^T (M-T 1. Z) B))))
1239
((EQUAL (CHECKSIGNTM VAR) '$NEGATIVE)
1240
(M*T (M^T (MCOS (SETQ Z (MATAN (MSQRT (M-T VAR)))))
1241
(SETQ B (M+T A B -1//2)))
1247
; List L contains two elements first the numerator parameter that
1248
;exceeds the denumerator one and is called "C", second
1249
;the difference of the two parameters which is called "M".
1251
(DEFUN DIFFINTPROP-GEN-EXEC (L L1 L2)
1252
(PROG (C M POLY CONSTFACT )
1255
L1 (ZL-DELETE C L1 1.)
1257
L2 (ZL-DELETE C L2 1.)
1258
POLY ($EXPAND (CONSTRPOLY C M 'AVGOUSTIS))
1259
CONSTFACT (CREATECONSTFACT C M))
1260
(RETURN (YANMULT CONSTFACT
1261
(DIFFINTPROP-EXEC POLY L1 L2)))))
1263
(DEFUN CONSTRPOLY (C M K)
1264
(COND ((ZEROP M) 1.)
1265
(T (MUL (ADD C K (SUB1 M)) (CONSTRPOLY C (SUB1 M) K)))))
1267
(DEFUN CREATECONSTFACT (C M)
1268
(COND ((ZEROP M) 1.)
1269
(T (MUL (INV (ADD C (SUB1 M)))
1270
(CREATECONSTFACT C (SUB1 M))))))
1272
(DEFUN DIFFINTPROP-EXEC (POLY L1 L2)
1273
(DISTRDIFFINTPROP (CREATECOEFPOWLIST-EXEC POLY) L1 L2))
1275
(DEFUN DISTRDIFFINTPROP (L L1 L2)
1277
(T (ADD (YANMULT ($FACTOR (CADAR L))
1278
(DIFFINTPROP (CAAR L) L1 L2))
1279
(DISTRDIFFINTPROP (CDR L) L1 L2)))))
1281
(DEFUN DIFFINTPROP (POW L1 L2)
1282
(COND ((ZEROP POW) (HGFSIMP L1 L2 VAR))
1284
(YANMULT (MUL (DIV (MULTPL L1) (MULTPL L2)) VAR)
1285
(HGFSIMP (INCR1 L1) (INCR1 L2) VAR)))
1286
(T (SEARCHADDSERIESLIST POW L1 L2))))
1288
(DEFUN SEARCHADDSERIESLIST (POW L1 L2)
1290
(COND ((SETQ SERIES (SEARCHSERIESLISTP SERIESLIST POW))
1291
(RETURN (EVAL SERIES))))
1300
'(YANMULT (MUL (DIV (MULTPL L1) (MULTPL L2))
1302
(DIFFINTPROPRECURSE (SUB1 POW)
1305
(RETURN (EVAL RES))))
1307
(DEFUN DIFFINTPROPRECURSE (POW L1 L2)
1310
($EXPAND (POWER (ADD 'AVGOUSTIS 1.) POW)))
1311
(RETURN (DIFFINTPROP-EXEC POLY L1 L2))))
1314
(COND ((NULL L) 1.) (T (MUL (CAR L) (MULTPL (CDR L))))))
1316
(DEFUN CREATECOEFPOWLIST-EXEC (POLY)
1318
(SETQ CONSTER (CONSTERMINIT POLY 'AVGOUSTIS)
1319
HP ($HIPOW POLY 'AVGOUSTIS))
1320
(RETURN (APPEND (LIST (LIST 0. CONSTER))
1321
(CREATECOEFPOWLIST POLY HP)))))
1323
(DEFUN CREATECOEFPOWLIST (POLY HP)
1324
(COND ((EQUAL HP 1.)
1325
(LIST (LIST 1. ($COEFF POLY 'AVGOUSTIS))))
1326
(T (APPEND (CREATECOEFPOWLIST POLY (SUB1 HP))
1332
(DEFUN CONSTERMINIT (FUN VAR)
1333
(COND ((EQ (CAAR FUN) 'MPLUS) (CONSTERM (CDR FUN) VAR))
1334
(T (COND ((FREEVAR FUN) FUN) (T 0.)))))
1336
(DEFUN SEARCHSERIESLISTP (SERIESLIST POW)
1337
(COND ((NULL SERIESLIST) NIL)
1338
((EQUAL (CAAR SERIESLIST) POW) (CADAR SERIESLIST))
1339
(T (SEARCHSERIESLISTP (CDR SERIESLIST) POW))))
1341
(DEFUN YANMULT (A B)
1342
(COND ((EQ (CAAR B) 'MPLUS) (YANMUL A (CDR B)))
1347
(T (ADD (MUL A (CAR B)) (YANMUL A (CDR B))))))
1350
(DEFUN FREEVARPAR(EXP)(COND ((FREEVAR EXP)(FREEPAR EXP))(T NIL)))
1352
(DECLARE-top (SPECIAL serieslist VAR PAR ZEROSIGNTEST PRODUCTCASE))
1355
(COND ((ATOM A) (NOT (EQ A VAR)))
1357
((AND (NOT (ATOM (CAR A)))
1358
(MEMQ 'ARRAY (CDAR A)))
1359
(COND ((FREEVAR (CDR A)) T)
1360
(T (PRINC 'VARIABLE-OF-INTEGRATION-APPEARED-IN-SUBSCRIPT)
1362
(T (AND (FREEVAR (CAR A)) (FREEVAR (CDR A))))))
1366
(COND ((ATOM EXP)(NOT (EQ EXP PAR)))
1367
(T (AND (FREEPAR (CAR EXP))(FREEPAR (CDR EXP))))))
1369
(DEFUN HASPAR(EXP)(COND ((FREEPAR EXP) NIL)(T T)))
1373
(PROG(A C A-C K M z)
1374
(SETQ A (CAR L1) C (CAR L2))
1375
(COND ((EQUAL C (ADD A A))
1377
(RETURN (MUL (POWER '$%E (setq z (DIV VAR 2)))
1378
(bestrig (add a (inv 2))
1379
(div (mul z z) 4))))))
1382
(COND ((NOT (MAXIMA-INTEGERP (SETQ A-C (SUB A C))))
1384
(COND ((MINUSP A-C)(RETURN (ERFGAMMARED A C VAR))))
1385
(RETURN (KUMMER L1 L2))
1387
(COND ((MAXIMA-INTEGERP A)(RETURN (KUMMER L1 L2))))
1391
(ADD (INV 2) M (MUL -1 A)))
1392
(RETURN (MUL (POWER VAR (MUL -1 (ADD (INV 2) M)))
1393
(POWER '$%E (DIV VAR 2))
1394
(WHITFUN K M VAR)))))
1398
(SETQ X (MUL '$%I (POWER X (INV 2))))
1399
(RETURN (MUL (POWER '$%PI (INV 2))
1402
(LIST '(%ERF) X)))))
1405
(COND ((AND (NUMP A)(NUMP C))(ERFGAMNUMRED A C Z))
1406
(T (GAMMAREDS A C Z))))
1409
(PROG(M NUMPROD RESULT COUNT ATEMP)
1411
(COND ((EQ M 1)(RETURN (HYPREDINCGM A Z))))
1422
(HYPREDINCGM ATEMP Z))
1425
(INV (SETQ ATEMP (ADD ATEMP 1)))
1426
(HYPREDINCGM ATEMP Z))))
1428
(COND ((EQ COUNT M)(RETURN RESULT)))
1435
(MUL (POWER -1 COUNT)
1436
(INV (FACTORIAL (SUB M
1440
(HYPREDINCGM ATEMP Z))))
1447
(POWER Z (MUL -1 A))
1448
(LIST '($%GAMMAGREEK) A Z)))))
1451
(COND ((EQ M 2) (MUL A (ADD A 1)))
1452
(T (MUL (ADD A (SUB1 M))(PROD A (SUB1 M))))))
1455
(COND ((MAXIMA-INTEGERP (SUB C (INV 2)))(ERFRED A C Z))
1456
(T (GAMMAREDS A C Z))))
1460
(SETQ N (SUB A (INV 2)) M (SUB C (DIV 3 2)))
1461
(COND ((NOT (OR (GREATERP N M)(MINUSP N)))
1462
(RETURN (THno33 N M Z))))
1463
(COND ((AND (MINUSP N)(MINUSP M))
1464
(RETURN (THno35 (MUL -1 N)(MUL -1 M) Z))))
1465
(COND ((AND (MINUSP N)(PLUSP M))
1466
(RETURN (THno34 (MUL -1 N) M Z))))
1467
(RETURN (GAMMAREDS (ADD N (INV 2))
1475
(MUL (DIV (MUL (POWER -1 M-N)
1476
(FCTRL (DIV 3 2) M-N)
1482
(MEVAL (LIST '($DIFF)
1485
(MEVAL (LIST '($DIFF)
1504
(DIV (MUL (FCTRL (DIV 3 2) M)
1505
(POWER '$%E 'YANNIS))
1509
(MEVAL (LIST '($DIFF)
1512
(MEVAL (LIST '($DIFF)
1517
(HYPREDERF 'YANNIS))
1526
(MUL (DIV (POWER 'YANNIS (SUB M (inv 2)))
1527
(MUL (POWER -1 (TIMES -1 M))
1529
(FCTRL (INV -2) M)))
1530
(MEVAL (LIST '($DIFF)
1531
(MUL (POWER 'YANNIS (inv 2))
1532
(POWER '$%E 'YANNIS)
1533
(MEVAL (LIST '($DIFF)
1540
(HYPREDERF 'YANNIS))
1549
(T (MUL (ADD A (SUB1 N))(FCTRL A (SUB1 N))))))
1551
(defun one (x)(equal x 1))
1557
(PROG (ASLIST QUEST ZEROSIGNTEST PRODUCTCASE)
1558
(SETQ ASLIST CHECKCOEFSIGNLIST)
1559
(COND ((ATOM EXPR) (GO LOOP)))
1560
(COND ((EQ (CAAR EXPR) 'MTIMES)
1561
(SETQ PRODUCTCASE T)))
1563
(COND ((NULL ASLIST)
1564
(SETQ CHECKCOEFSIGNLIST
1565
(APPEND CHECKCOEFSIGNLIST
1574
(COND ((EQUAL (CAAR ASLIST) EXPR)
1575
(RETURN (CADAR ASLIST))))
1576
(SETQ ASLIST (CDR ASLIST))
1579
(DEFUN CHECKFLAGANDACT
1581
(COND (PRODUCTCASE (SETQ PRODUCTCASE NIL)
1582
(FINDSIGNOFTHEIRPRODUCT (FINDSIGNOFACTORS
1584
(T (ASKSIGN ($REALPART EXPR)))))
1586
(DEFUN FINDSIGNOFACTORS
1588
(COND ((NULL LISTOFACTORS) NIL)
1589
((EQ ZEROSIGNTEST '$ZERO) '$ZERO)
1590
(T (APPEND (LIST (SETQ ZEROSIGNTEST
1593
(FINDSIGNOFACTORS (CDR LISTOFACTORS))))))
1595
(DEFUN FINDSIGNOFTHEIRPRODUCT
1598
(COND ((EQ LIST '$ZERO) (RETURN '$ZERO)))
1599
(SETQ SIGN '$POSITIVE)
1601
(COND ((NULL LIST) (RETURN SIGN)))
1602
(COND ((EQ (CAR LIST) '$POSITIVE)
1603
(SETQ LIST (CDR LIST))
1605
(COND ((EQ (CAR LIST) '$NEGATIVE)
1615
(COND ((EQ SIGN '$POSITIVE) '$NEGATIVE) (T '$POSITIVE)))
1620
(DEFUN VFVP(EXP)(M2 EXP '(V FREEVARPAR) NIL))
1626
'((MTIMES)((COEFFTT)(D FREEPAR))((COEFFTT)(U HASPAR)))
1632
(LIST '($%F ARRAY)(LENGTH L1)(LENGTH L2))
1633
(APPEND (LIST '(MLIST)) L1)
1634
(APPEND (LIST '(MLIST)) L2)
1641
(prog(result prodnum proden count k a1 b1)
1661
(hgfsimp l1 l2 var))
1663
(cond ((eq count k) (return result)))
1667
(mul prodnum (mull l1))
1669
(mul proden (mull l2))
1672
(mul (combin k count)
1673
(div prodnum proden)
1675
(hgfsimp (setq l1 (incr1 l1))
1676
(setq l2 (incr1 l2))
1683
(mul (factorial count)(factorial (sub k count)))))
1686
;Algor. II from thesis:minimizes differentiations
1688
(prog (m n ap con sym m+n)
1689
(cond ((not (setq sym (cdras 'f (s+c a))))
1691
(setq con (sub a sym))
1693
(setq m+n (add a b))
1694
(setq m ($entier con))
1695
(cond ((minusp m)(add1 m)))
1696
(setq ap (add (sub con m) ap))
1698
(cond ((and (minusp (mul n m))(greaterp (abs m) (abs n)))
1699
(return (list ap (sub ap n) m+n))))
1700
(return (list ap (add ap m) m+n))))
1706
;Algor. 2F1-RL from thesis:step 4:dispatch on a+m,-a+n,1/2+l cases
1709
(prog (aprime m n $ratsimpexponens $ratprint newf)
1718
(setq $ratsimpexponens $true $ratprint $false)
1720
($ratsimp (subst aprime
1731
(return (subst var 'ell
1732
(algiii (subst 'ell var newf)
1735
;Pattern match for s(ymbolic) + c(onstant) in parameter
1739
'((MPLUS)((COEFFPT)(F nonnump))((COEFFPP)(C $numberp)))
2034
(sub 1 (div 1 var))))
2035
(mul (gm a)(gm b)))))
2039
(defun trig-log (arg-l1 arg-l2)
2040
(cond ((equal (simplifya (car arg-l2) nil) '((rat simp) 3 2))
2043
(format t " trig-log: Test c=3/2~%"))
2044
(trig-log-3 arg-l1 arg-l2))
2045
((equal (simplifya (car arg-l2) nil) '((rat simp) 1 2))
2048
(format t " trig-log: Test c=1/2~%"))
2049
(trig-log-1 arg-l1 arg-l2))
2053
(defun trig-log-3 (arg-l1 arg-l2)
2054
(cond ((and (or (equal (car arg-l1) 1) (equal (cadr arg-l1) 1))
2055
(or (equal (car arg-l1) (div 1 2))
2056
(equal (cadr arg-l1) (div 1 2))))
2057
;; (a = 1 or b = 1) and (a = 1/2 or b = 1/2)
2059
(format t " Case a or b is 1 and the other is 1/2~%"))
2060
(trig-log-3-exec arg-l1 arg-l2))
2061
((and (equal (car arg-l1) (cadr arg-l1))
2062
(or (equal 1 (car arg-l1))
2063
(equal (div 1 2) (car arg-l1))))
2064
;; a = b and (a = 1 or a = 1/2)
2066
(format t " Case a=b and a is 1 or 1/2~%"))
2067
(trig-log-3a-exec arg-l1 arg-l2))
2068
((or (equal (add (car arg-l1) (cadr arg-l1)) 1)
2069
(equal (add (car arg-l1) (cadr arg-l1)) 2))
2070
;; a + b = 1 or a + b = 2
2072
(format t " Case a+b is 1 or 2~%"))
2073
(trig-sin arg-l1 arg-l2))
2074
((or (equal (sub (car arg-l1) (cadr arg-l1)) (div 1 2))
2075
(equal (sub (cadr arg-l1) (car arg-l1)) (div 1 2)))
2076
;; a - b = 1/2 or b - a = 1/2
2078
(format t " Case a-b=1/2 or b-a=1/2~%"))
2079
(trig-3 arg-l1 arg-l2))
2082
(defun trig-3 (arg-l1 arg-l2)
2083
(declare (ignore arg-l2))
2086
;; F(a,a+1/2,3/2,z^2) =
2087
;; ((1+z)^(1-2*a) - (1-z)^(1-2*a))/2/z/(1-2*a)
2088
(let* (($radexpand '$all)
2090
(sub (add (car arg-l1)
2093
(z (power var (div 1 2))))
2097
(sub (power (add 1 z) a)
2098
(power (sub 1 z) a)))))
2100
(defun trig-sin (arg-l1 arg-l2)
2101
(declare (ignore arg-l2))
2102
;; A&S 15.1.15, 15.1.16
2103
(destructuring-bind (a b)
2105
;; I think it's ok to convert sqrt(z^2) to z here, so $radexpand
2107
(let (($radexpand '$all)
2109
(cond ((equal (add a b) 1)
2112
;; F(a,1-a;3/2;sin(z)^2) =
2114
;; sin((2*a-1)*z)/(2*a-1)/sin(z)
2115
(mul (inv (mul (mul -1 (sub a b))
2116
(msin (masin (msqrt var)))))
2119
(masin (msqrt var))))))
2120
((equal (add a b) 2)
2123
;; F(a, 2-a; 3/2; sin(z)^2) =
2125
;; sin((2*a-2)*z)/(a-1)/sin(2*z)
2126
(mul (msin (mul (setq z1
2140
;;Generates atan if arg positive else log
2141
(defun trig-log-3-exec (arg-l1 arg-l2)
2142
(declare (ignore arg-l1 arg-l2))
2143
;; See A&S 15.1.4 and 15.1.5
2145
;; F(a,b;3/2;z) where a = 1/2 and b = 1 (or vice versa).
2147
;; I think it's ok to convert sqrt(z^2) to z here, so $radexpand is
2149
(let (($radexpand '$all))
2150
(cond ((equal (checksigntm var) '$positive)
2153
;; F(1/2,1;3/2,z^2) =
2155
;; log((1+z)/(1-z))/z/2
2157
;; This is the same as atanh(z)/z. Should we return that
2158
;; instead? This would make this match what hyp-atanh
2160
(let ((z (power var (div 1 2))))
2163
(mlog (div (add 1 z)
2165
((equal (checksigntm var) '$negative)
2168
;; F(1/2,1;3/2,z^2) =
2170
(let ((z (power (mul -1 var)
2175
(defun trig-log-3a-exec (arg-l1 arg-l2)
2176
;; See A&S 15.1.6 and 15.1.7
2178
;; F(a,b;3/2,z) where a = b and a = 1/2 or a = 1.
2180
;; I think it's ok to convert sqrt(z^2) to z here, so $radexpand is
2182
(let ((a (first arg-l1))
2184
(cond ((equal (checksigntm var) '$positive)
2187
;; F(1/2,1/2; 3/2; z^2) = sqrt(1-z^2)*F(1,1;3/2;z^2) =
2189
(let ((z (power var (div 1 2))))
2191
(div (trig-log-3a-exec (list (div 1 2) (div 1 2)) arg-l2)
2192
(power (sub 1 (power z 2)) (div 1 2)))
2193
(div (masin z) z))))
2194
((equal (checksigntm var) '$negative)
2197
;; F(1/2,1/2; 3/2; -z^2) = sqrt(1+z^2)*F(1,1,3/2; -z^2) =
2198
;;log(z + sqrt(1+z^2))/z
2199
(let* ((z (power (mul -1 var)
2201
(1+z^2 (add 1 (power z 2))))
2203
(div (trig-log-3a-exec (list (div 1 2) (div 1 2))
2207
(div (mlog (add z (power 1+z^2
2212
(defun trig-log-1 (arg-l1 arg-l2) ;; 2F1's with C = 1/2
2213
(declare (ignore arg-l2))
2214
;; 15.1.17, 11, 18, 12, 9, and 19
2216
;; I think it's ok to convert sqrt(z^2) to z here, so $radexpand is
2218
(let (($radexpand '$all)
2219
x z $exponentialize a b)
2220
(setq a (car arg-l1) b (cadr arg-l1))
2221
(cond ((=0 (m+t a b))
2224
(cond ((equal (checksigntm var) '$positive)
2226
;; F(-a,a;1/2;sin(z)^2) = cos(2*a*z)
2227
(trig-log-1-pos a var))
2228
((equal (checksigntm var) '$negative)
2230
;; F(-a,a;1/2;-z^2) = 1/2*((sqrt(1+z^2)+z)^(2*a)
2231
;; +(sqrt(1+z^2)-z)^(2*a))
2233
(trig-log-1-neg a b var))
2235
((equal (m+t a b) 1.)
2237
(cond ((equal (checksigntm var) '$positive)
2239
;; F(a,1-a;1/2;sin(z)^2) = cos((2*a-1)*z)/cos(z)
2240
(m//t (mcos (m*t (m-t a b) (setq z (masin (msqrt var)))))
2242
((equal (checksigntm var) '$negative)
2244
;; F(a,1-a;1/2;-z^2) = 1/2*(1+z^2)^(-1/2)*
2245
;; {[(sqrt(1+z^2)+z]^(2*a-1)
2246
;; +[sqrt(1+z^2)-z]^(2*a-1)}
2247
(m*t 1//2 (m//t (setq x (msqrt (m-t 1. var))))
2248
(m+t (m^t (m+t x (setq z (msqrt (m-t var))))
2250
(m^t (m-t x z) b))))
2252
((=1//2 (hyp-mabs (m-t b a)))
2253
;; F(a, a+1/2; 1/2; z)
2254
(cond ((equal (checksigntm var) '$positive)
2256
;; F(a,1/2+a;1/2;z^2) = ((1+z)^(-2*a)+(1-z)^(-2*a))/2
2258
(m+t (m^t (m1+t (setq z (msqrt var)))
2259
(setq b (m-t 1//2 (m+t a b))))
2260
(m^t (m-t 1. z) b))))
2261
((equal (checksigntm var) '$negative)
2263
;; F(a,1/2+a;1/2;-tan(z)^2) = cos(z)^(2*a)*cos(2*a*z)
2264
(m*t (m^t (mcos (setq z (matan (msqrt (m-t var)))))
2265
(setq b (m+t a b -1//2)))
2270
(defun trig-log-1-pos (a z)
2271
;; I think it's ok to convert sqrt(z^2) to z here, so $radexpand is
2273
(let (($radexpand '$all))
2274
(mcos (m*t 2. a (masin (msqrt z))))))
2276
(defun trig-log-1-neg (a b v)
2277
;; Look to see a is of the form m*s+c where m and c
2278
;; are numbers. If m is positive, swap a and b.
2279
;; Basically we want F(-a,a;1/2;-z^2) =
2280
;; F(a,-a;1/2;-z^2), as they should be.
2281
(let* (($radexpand '$all)
2283
(m (cdras 'm match))
2284
(s (cdras 's match))
2286
(if (and m (plusp m))
2289
(if (eq (checksigntm a) '$negative)
2292
(x (msqrt (m-t 1. v)))
2293
(z (msqrt (m-t v))))
2296
(setq b (m*t 2. b)))
2297
(m^t (m-t x z) b)))))
2299
;; Pattern match for m*s+c where a is a number, x is symbolic, and c
2303
'((mplus) ((coeffpt) (m $numberp) (s nonnump))
2304
((coeffpp) (c $numberp)))
2307
;; List L contains two elements first the numerator parameter that
2308
;;exceeds the denumerator one and is called "C", second
2309
;;the difference of the two parameters which is called "M".
2312
(defun diffintprop-gen-exec (l l1 l2)
2313
(prog (c m poly constfact )
2316
l1 (zl-delete c l1 1.)
2318
l2 (zl-delete c l2 1.)
2319
poly ($expand (constrpoly c m 'avgoustis))
2320
constfact (createconstfact c m))
2321
(return (yanmult constfact
2322
(diffintprop-exec poly l1 l2)))))
2324
(defun constrpoly (c m k)
2325
(cond ((zerop m) 1.)
2326
(t (mul (add c k (sub1 m)) (constrpoly c (sub1 m) k)))))
2328
(defun createconstfact (c m)
2329
(cond ((zerop m) 1.)
2330
(t (mul (inv (add c (sub1 m)))
2331
(createconstfact c (sub1 m))))))
2333
(defun diffintprop-exec (poly l1 l2)
2334
(distrdiffintprop (createcoefpowlist-exec poly) l1 l2))
2336
(defun distrdiffintprop (l l1 l2)
2338
(t (add (yanmult ($factor (cadar l))
2339
(diffintprop (caar l) l1 l2))
2340
(distrdiffintprop (cdr l) l1 l2)))))
2342
(defun diffintprop (pow l1 l2)
2343
(cond ((zerop pow) (hgfsimp l1 l2 var))
2345
(yanmult (mul (div (multpl l1) (multpl l2)) var)
2346
(hgfsimp (incr1 l1) (incr1 l2) var)))
2347
(t (searchaddserieslist pow l1 l2))))
2349
(defun searchaddserieslist (pow l1 l2)
2351
(cond ((setq series (searchserieslistp serieslist pow))
2352
(return (eval series))))
2361
'(yanmult (mul (div (multpl l1) (multpl l2))
2363
(diffintproprecurse (sub1 pow)
2366
(return (eval res))))
2368
(defun diffintproprecurse (pow l1 l2)
2371
($expand (power (add 'avgoustis 1.) pow)))
2372
(return (diffintprop-exec poly l1 l2))))
2375
(cond ((null l) 1.) (t (mul (car l) (multpl (cdr l))))))
2377
(defun createcoefpowlist-exec (poly)
2379
(setq conster (consterminit poly 'avgoustis)
2380
hp ($hipow poly 'avgoustis))
2381
(return (append (list (list 0. conster))
2382
(createcoefpowlist poly hp)))))
2384
(defun createcoefpowlist (poly hp)
2385
(cond ((equal hp 1.)
2386
(list (list 1. ($coeff poly 'avgoustis))))
2387
(t (append (createcoefpowlist poly (sub1 hp))
2393
(defun consterminit (fun var)
2394
(cond ((eq (caar fun) 'mplus) (consterm (cdr fun) var))
2395
(t (cond ((freevar fun) fun) (t 0.)))))
2397
(defun searchserieslistp (serieslist pow)
2398
(cond ((null serieslist) nil)
2399
((equal (caar serieslist) pow) (cadar serieslist))
2400
(t (searchserieslistp (cdr serieslist) pow))))
2402
(defun yanmult (a b)
2403
(cond ((eq (caar b) 'mplus) (yanmul a (cdr b)))
2408
(t (add (mul a (car b)) (yanmul a (cdr b))))))
2412
(defun freevarpar (exp)
2413
(cond ((freevar exp) (freepar exp))
2416
;; Why is this needed?
2419
;;(DEFUN FREEVAR (A)
2420
;; (COND ((ATOM A) (NOT (EQ A VAR)))
2421
;; ((ALIKE1 A VAR)NIL)
2422
;; ((AND (NOT (ATOM (CAR A)))
2423
;; (MEMQ 'ARRAY (CDAR A)))
2424
;; (if (FREEVAR (CDR A))
2426
;; (merror "`variable-of-integration-appeared-in-subscript'")))
2427
;; (T (AND (FREEVAR (CAR A)) (FREEVAR (CDR A))))))
2429
(defun freepar (exp)
2431
(not (eq exp *par*)))
2432
(t (and (freepar (car exp))
2433
(freepar (cdr exp))))))
2436
(cond ((freepar exp) nil)
2439
;; Confluent hypergeometric function.
2442
(defun confl (arg-l1 arg-l2 var)
2443
(let* ((a (car arg-l1))
2446
(cond ((alike1 c (add a a))
2450
;; F(n+1;2*n+1;2*z) =
2451
;; gamma(3/2+n)*exp(z)*(z/2)^(-n-1/2)*bessel_i(n+1/2,z).
2456
;; gamma(n+1/2)*exp(z/2)*(z/4)^(-n-3/2)*bessel_i(n-1/2,z/2);
2457
(let ((z (div var 2)))
2459
(bestrig (add a (inv 2))
2460
(div (mul z z) 4)))))
2461
((not (hyp-integerp a-c))
2462
(cond ((hyp-integerp a)
2463
(kummer arg-l1 arg-l2))
2467
;; %m[k,u](z) = exp(-z/2)*z^(u+1/2)*M(1/2+u-k,1+2*u,z)
2471
;; M(a,c,z) = exp(z/2)*z^(-c/2)*%m[c/2-a,c/2-1/2](z)
2473
;; But if we apply Kummer's transformation, we can also
2474
;; derive the expression
2476
;; %m[k,u](z) = exp(z/2)*z^(u+1/2)*M(1/2+u+k,1+2*u,-z)
2480
;; M(a,c,z) = exp(-z/2)*(-z)^(-c/2)*%m[a-c/2,c/2-1/2](-z)
2482
;; For Laplace transforms it might be more beneficial to
2483
;; return this second form instead.
2484
(let* ((m (div (sub c 1) 2))
2485
(k (add (inv 2) m (mul -1 a))))
2486
(mul (power var (mul -1 (add (inv 2) m)))
2487
(power '$%e (div var 2))
2488
(whitfun k m var))))
2491
(fpqform arg-l1 arg-l2 var))))
2493
(erfgammared a c var))
2495
(kummer arg-l1 arg-l2)))))
2498
;; M(1/2,3/2,-z^2) = sqrt(%pi)*erf(z)/2/sqrt(z)
2500
;; So M(1/2,3/2,z) = sqrt(%pi)*erf(sqrt(-z))/2/sqrt(-z)
2501
;; = sqrt(%pi)*erf(%i*sqrt(z))/2/(%i*sqrt(z))
2502
(defun hyprederf (x)
2503
(let ((x (mul '$%i (power x (inv 2)))))
2504
(mul (power '$%pi (inv 2))
2509
;; M(a,c,z), where a-c is a negative integer.
2510
(defun erfgammared (a c z)
2511
(cond ((and (nump a) (nump c))
2512
(erfgamnumred a c z))
2513
(t (gammareds a c z))))
2515
;; M(a,c,z) where a-c is a negative integer, and at least one of a or
2516
;; c is not a number.
2518
(defun gammareds (a c z)
2519
(prog (m numprod result count atemp)
2523
;; We have M(a,a+1,z)
2524
(return (hypredincgm a z))))
2525
(setq numprod (prod a m)
2531
(hypredincgm atemp z))
2534
(inv (setq atemp (add atemp 1)))
2535
(hypredincgm atemp z))))
2537
(cond ((eq count m)(return result)))
2538
(setq count (add1 count)
2541
(mul (power -1 count)
2542
(inv (factorial (sub m
2546
(hypredincgm atemp z))))
2549
;; I (rtoy) think this is what the function above is doing, but I'm
2550
;; not sure. Plus, I think it's wrong.
2552
;; For hgfred([n],[2+n],-z), the above returns
2554
;; 2*n*(n+1)*z^(-n-1)*(%gammagreek(n,z)*z-%gammagreek(n+1,z))
2556
;; But from A&S 13.4.3
2558
;; -M(n,2+n,z) - n*M(n+1,n+2,z) + (n+1)*M(n,n+1,z) = 0
2560
;; so M(n,2+n,z) = (n+1)*M(n,n+1,z)-n*M(n+1,n+2,z)
2562
;; And M(n,n+1,-z) = n*z^(-n)*%gammagreek(n,z)
2566
;; M(n,2+n,z) = (n+1)*n*z^(-n)*%gammagreek(n,z) - n*(n+1)*z^(-n-1)*%gammagreek(n+1,z)
2567
;; = n*(n+1)*z^(-n-1)*(%gammagreek(n,z)*n-%gammagreek(n+1,z))
2569
;; So the version above is off by a factor of 2. But I think it's more than that.
2570
;; Using A&S 13.4.3 again,
2572
;; M(n,n+3,-z) = [n*M(n+1,n+3,-z) - (n+2)*M(n,n+2,-z)]/(-2);
2574
;; The version above doesn't produce anything like this equation would
2575
;; produce, given the value of M(n,n+2,-z) derived above.
2576
(defun gammareds (a c z)
2577
;; M(a,c,z) where a-c is a negative integer.
2578
(let ((diff (sub c a)))
2580
;; We have M(a,a+1,z).
2584
;; Apply Kummer's tranformation to get the form M(a-1,a,z)
2586
;; (I don't think we ever get here, but just in case, we leave it.)
2588
(kummer (list a) (list c))))
2590
;; We have M(a, a+n, z)
2593
;; (1+a-b)*M(a,b,z) - a*M(a+1,b,z)+(b-1)*M(a,b-1,z) = 0
2597
;; M(a,b,z) = [a*M(a+1,b,z) - (b-1)*M(a,b-1,z)]/(1+a-b);
2599
;; Thus, the difference between b and a is reduced, until
2600
;; b-a=1, which we handle above.
2602
(gammareds (add 1 a) c z))
2604
(gammareds a (sub c 1) z)))
2605
(inv (sub (add 1 a) c)))))))
2607
;; Incomplete gamma function
2609
;; gammagreek(a,x) = integrate(t^(a-1)*exp(-t),t,0,x)
2610
(defun gammagreek (a z)
2611
(cond ((and (integerp a) (eql a 1))
2612
;; gammagreek(0, x) = 1-exp(x)
2613
(sub 1 (mexpt (neg z))))
2614
((and (integerp a) (plusp a))
2615
;; gammagreek(a,z) can be simplified if a is a positive
2620
;; gammagreek(a+1,x) = a*gammagreek(a,x) - x^a*exp(-x)
2624
;; gammagreek(a,x) = (a-1)*gammagreek(a-1,x)-x^(a-1)*exp(-x)
2625
(let ((a-1 (sub a 1)))
2626
(sub (mul a-1 (gammagreek a-1 z))
2632
;; gammagreak(1/2,x^2) = sqrt(%pi)*erf(x)
2634
;; gammagreek(1/2,z) = sqrt(%pi)*erf(sqrt(x))
2637
`((%erf) ,(msqrt z))))
2638
((and (integerp (add a 1//2)))
2639
;; gammagreek(n+1/2,z) can be simplified using A&S 6.5.22 to
2640
;; reduce the problem to gammagreek(1/2,x), which we know,
2642
(if (ratgreaterp a 0)
2643
(let ((a-1 (sub a 1)))
2644
(sub (mul a-1 (gammagreek a-1 z))
2647
(let ((a+1 (add a 1)))
2648
(div (add (gammagreek a+1 z)
2654
`(($%gammagreek) ,a ,z))))
2657
;; %gammagreek(a,x) = x^a/a*M(a,1+a,-x)
2658
;; = x^a/a*exp(-x)*M(1,1+a,x)
2660
;; where %gammagreek(a,x) is the incomplete gamma function.
2662
;; M(a,1+a,x) = a*(-x)^(-a)*%gammagreek(a,-x)
2669
(power z (mul -1 a))
2670
(list '($%gammagreek) a z)))))
2672
(defun hypredincgm (a z)
2673
(let ((-z (mul -1 z)))
2674
(mul a (power -z (mul -1 a))
2675
(gammagreek a -z))))
2677
;; M(a,c,z), when a and c are numbers, and a-c is a negative integer
2678
(defun erfgamnumred (a c z)
2679
(cond ((hyp-integerp (sub c (inv 2)))
2681
(t (gammareds a c z))))
2683
;; M(a,c,z) when a and c are numbers and c-1/2 is an integer and a-c
2684
;; is a negative integer. Thus, we have M(p+1/2, q+1/2,z)
2685
(defun erfred (a c z)
2687
(setq n (sub a (inv 2))
2688
m (sub c (div 3 2)))
2691
;; a - c < 0 so n - m - 1 < 0
2692
(cond ((not (or (greaterp n m) (minusp n)))
2694
(return (thno33 n m z))))
2695
(cond ((and (minusp n) (minusp m))
2697
(return (thno35 (mul -1 n) (mul -1 m) z))))
2698
(cond ((and (minusp n) (plusp m))
2700
(return (thno34 (mul -1 n) m z))))
2702
(return (gammareds (add n (inv 2))
2705
;; Compute M(n+1/2, m+3/2, z) with 0 <= n <= m.
2707
;; I (rtoy) think this is what this routine is doing. (I'm guessing
2708
;; that thno33 means theorem number 33 from Yannis Avgoustis' thesis.)
2710
;; I don't have his thesis, but I see there are similar ways to derive
2711
;; the result we want.
2714
;; Use Kummer's transformation (A&S ) to get
2716
;; M(n+1/2,m+3/2,z) = exp(z)*M(m-n+1,m+3/2,-z)
2718
;; From A&S, we have
2720
;; diff(M(1,n+3/2,z),z,m-n) = poch(1,m-n)/poch(n+3/2,m-n)*M(m-n+1,m+3/2,z)
2722
;; Apply Kummer's transformation again:
2724
;; M(1,n+3/2,z) = exp(z)*M(n+1/2,n+3/2,-z)
2726
;; Apply the differentiation formula again:
2728
;; diff(M(1/2,3/2,z),z,n) = poch(1/2,n)/poch(3/2,n)*M(n+1/2,n+3/2,z)
2730
;; And we know that M(1/2,3/2,z) can be expressed in terms of erf.
2734
;; Since n <= m, apply the differentiation formula:
2736
;; diff(M(1/2,m-n+3/2,z),z,n) = poch(1/2,n)/poch(m-n+3/2,n)*M(n+1/2,m+3/2,z)
2738
;; Apply Kummer's transformation:
2740
;; M(1/2,m-n+3/2,z) = exp(z)*M(m-n+1,m-n+3/2,z)
2742
;; Apply the differentiation formula again:
2744
;; diff(M(1,3/2,z),z,m-n) = poch(1,m-n)/poch(3/2,m-n)*M(m-n+1,m-n+3/2,z)
2746
;; I think this routine uses Method 2.
2753
(mul (div (mul (power -1 m-n)
2754
(fctrl (div 3 2) m-n)
2760
;; diff(M(1/2,m-n+3/2,z),z,n)
2761
(meval (list '($diff)
2762
;; Kummer's transformation
2765
;; diff(M(1,3/2,z),z,m-n)
2766
(meval (list '($diff)
2767
;; M(1,3/2,-z) = e^(-z)*M(1/2,3/2,z)
2782
(defun thno33 (n m x)
2783
;; M(n+1/2,m+3/2,z) = diff(M(1/2,m-n+3/2,z),z,n)*poch(m-n+3/2,n)/poch(1/2,n)
2784
;; M(1/2,m-n+3/2,z) = exp(z)*M(m-n+1,m-n+3/2,-z)
2785
;; M(m-n+1,m-n+3/2,z) = diff(M(1,3/2,z),z,m-n)*poch(3/2,m-n)/poch(1,m-n)
2786
;; diff(M(1,3/2,z),z,m-n) = (-1)^(m-n)*diff(M(1,3/2,-z),z,m-n)
2787
;; M(1,3/2,-z) = exp(-z)*M(1/2,3/2,z)
2788
(let* ((m-n (sub m n))
2789
;; poch(m-n+3/2,n)/poch(1/2,n)
2790
(factor1 (div (fctrl (add m-n (div 3 2)) n)
2792
;; poch(3/2,m-n)/poch(1,m-n)
2793
(factor2 (div (fctrl (div 3 2) m-n)
2795
;; M(1,3/2,-z) = exp(-z)*M(1/2,3/2,z)
2796
(hgferf (mul (power '$%e (mul -1 'yannis))
2797
(hyprederf 'yannis)))
2798
;; diff(M(1,3/2,z),z,m-n)
2799
(diff1 (meval `(($diff) ,hgferf 'yannis ,m-n)))
2800
;; exp(z)*M(m-n+1,m-n+3/2,-z)
2801
(kummer (mul (power '$%e 'yannis)
2803
;; diff(M(1/2,m-n+3/2,z),z,n)
2804
(diff2 (meval `(($diff) ,kummer 'yannis ,n))))
2805
;; Multiply all the terms together.
2809
(subst x 'yannis diff2))))
2811
;; M(n+1/2,m+3/2,z), with n < 0 and m > 0
2813
;; Let's write it more explicitly as M(-n+1/2,m+3/2,z) with n > 0 and
2816
;; First, use Kummer's transformation to get
2818
;; M(-n+1/2,m+3/2,z) = exp(z)*M(m+n+1,m+3/2,-z)
2822
;; diff(z^(n+m)*M(m+1,m+3/2,z),z,n) = poch(m+1,n)*z^m*M(m+n+1,m+3/2,z)
2826
;; diff(M(1,3/2,z),z,m) = poch(1,m)/poch(3/2,m)*M(m+1,m+3/2,z)
2828
;; Thus, we can compute M(-n+1/2,m+3/2,z) from M(1,3/2,z).
2830
;; The second formula above can be derived easily by multiplying the
2831
;; series for M(m+1,m+3/2,z) by z^(n+m) and differentiating n times.
2833
(defun thno34 (n m x)
2837
(div (mul (fctrl (div 3 2) m)
2838
(power '$%e 'yannis))
2842
(meval (list '($diff)
2845
(meval (list '($diff)
2850
(hyprederf 'yannis))
2856
;; M(n+1/2,m+3/2,z), with n < 0 and m < 0
2858
;; Write it more explicitly as M(-n+1/2,-m+3/2,z) with n > 0 and m >
2863
;; diff(sqrt(z)*M(-n+1/2,3/2,z),z,m) = poch(3/2-m,m)*M(-n+1/2,-m+3/2,z).
2865
;; Apply Kummer's transformation:
2867
;; M(-n+1/2,3/2,z) = exp(z) * M(n+1,3/2,-z)
2871
;; diff(z^n*M(1,3/2,z),z,n) = n!*M(n+1,3/2,z)
2873
;; So we can express M(-n+1/2,-m+3/2,z) in terms of M(1,3/2,z).
2875
;; The first formula above follows from the more general formula
2877
;; diff(z^(b-1)*M(a,b,z),z,n) = poch(b-n,n)*z^(b-n-1)*M(a,b-n,z)
2879
;; The last formula follows from the general result
2881
;; diff(z^(a+n-1)*M(a,b,z),z,n) = poch(a,n)*z^(a-1)*M(a+n,b,z)
2883
;; Both of these are easily derived by using the series for M and
2885
(defun thno35 (n m x)
2888
(mul (div (power 'yannis (sub m (inv 2)))
2889
(mul (power -1 (times -1 m))
2891
(fctrl (inv -2) m)))
2892
(meval (list '($diff)
2893
(mul (power 'yannis (inv 2))
2894
(power '$%e 'yannis)
2895
(meval (list '($diff)
2902
(hyprederf 'yannis))
2908
;; Pochhammer symbol. fctrl(a,n) = a*(a+1)*(a+2)*...*(a+n-1).
2910
;; N must be a positive integer!
2912
;; FIXME: This appears to be identical to factf below.
2919
(mul (add a (sub1 n))
2920
(fctrl a (sub1 n))))))
2927
(m2 exp '(v freevarpar) nil))
2933
'((mtimes)((coefftt)(d freepar))((coefftt)(u haspar)))
2936
(defun fpqform (arg-l1 arg-l2 arg)
2938
(list '($%f array) (length arg-l1)(length arg-l2))
2939
(append (list '(mlist)) arg-l1)
2940
(append (list '(mlist)) arg-l2)
2945
;; Consider pFq([a_k]; [c_j]; z). If a_k = c_j + m for some k and j
2946
;; and m >= 0, we can express pFq in terms of (p-1)F(q-1).
2948
;; Here is a derivation for F(a,b;c;z), but it generalizes to the
2949
;; generalized hypergeometric very easily.
2953
;; diff(z^(a+n-1)*F(a,b;c;z), z, n) = poch(a,n)*z^(a-1)*F(a+n,b;c;z)
2955
;; F(a+n,b;c;z) = diff(z^(a+n-1)*F(a,b;c;z), z, n)/poch(a,n)/z^(a-1)
2958
;; So this expresses F(a+n,b;c;z) in terms of F(a,b;c;z). Let a = c +
2959
;; n. This therefore gives F(c+n,b;c;z) in terms of F(c,b;c;z) =
2960
;; 1F0(b;;z), which we know.
2962
;; For simplicity, we will write F(z) for F(a,b;c;z).
2967
;; diff(z^x*F(z),z,n) = sum binomial(n,k)*diff(z^x,z,n-k)*diff(F(z),z,k)
2970
;; But diff(z^x,z,n-k) = x*(x-1)*...*(x-n+k+1)*z^(x-n+k)
2971
;; = poch(x-n+k+1,n-k)*z^(x-n+k)
2975
;; z^(-a+1)/poch(a,n)*diff(z^(a+n-1),z,n-k)
2976
;; = poch(a+n-1-n+k+1,n-k)/poch(a,n)*z^(a+n-1-n+k)*z^(-a+1)
2977
;; = poch(a+k,n-k)/poch(a,n)*z^k
2980
;; Combining these we have
2983
;; F(a+n,b;c;z) = sum z^k/poch(a,k)*binomial(n,k)*diff(F(a,b;c;z),z,k)
2986
;; Since a = c, we have
2989
;; F(a+n,b;a;z) = sum z^k/poch(a,k)*binomial(n,k)*diff(F(a,b;a;z),z,k)
2992
;; But F(a,b;a;z) = F(b;;z) and it's easy to see that A&S 15.2.2 can
2993
;; be specialized to this case to give
2995
;; diff(F(b;;z),z,k) = poch(b,k)*F(b+k;;z)
2997
;; Finally, combining all of these, we have
3000
;; F(a+n,b;c;z) = sum z^k/poch(a,k)*binomial(n,k)*poch(b,k)*F(b+k;;z)
3003
;; Thus, F(a+n,b;c;z) is expressed in terms of 1F0(b+k;;z), as desired.
3004
(defun splitpfq (l arg-l1 arg-l2)
3005
(destructuring-bind (a1 k)
3012
(arg-l1 (zl-delete a1 arg-l1 1))
3013
(arg-l2 (zl-delete b1 arg-l2 1)))
3014
(loop for count from 0 upto k
3017
(format t "splitpfg term:~%")
3018
(maxima-display (mul (combin k count)
3019
(div prodnum proden)
3022
(format t "F(~A, ~A)~%" arg-l1 arg-l2))
3023
(setq result (add result
3024
(mul (combin k count)
3025
(div prodnum proden)
3028
(hgfsimp arg-l1 arg-l2 var))))
3029
(setq prod-b (mul prod-b (add b1 count)))
3030
(setq prodnum (mul prodnum (mull arg-l1))
3031
proden (mul proden (mull arg-l2)))
3032
(setq arg-l1 (incr1 arg-l1))
3033
(setq arg-l2 (incr1 arg-l2)))
3036
;; binomial(k,count)
3037
(defun combin (k count)
3039
(mul (factorial count)
3040
(factorial (sub k count)))))
3043
;; Algor. II from thesis:minimizes differentiations
3045
;; We're looking at F(a+m,-a+n;1/2+L;z)
3047
(defun algii (a b c)
3048
(declare (ignore c))
3049
(prog (m n ap con sym m+n)
3050
;; We know a+b is an integer. In the most general form, we can
3051
;; have a = r*s+f+m and b = -(r*s+f)+n.
3052
(cond ((not (setq sym (cdras 'f (s+c a))))
3054
(setq con (sub a sym))
3056
(setq m+n (add a b))
3057
(setq m ($entier con))
3060
;; At this point sym = r*s, con is f+m, and m is m.
3061
(setq ap (add (sub con m) ap))
3064
;; Return r*a+f, r*a+f+p, and m+n, where p is chosen to minimize
3065
;; the number of derivatives we need to take. Basically
3066
;; p=min(abs(m),abs(n)).
3068
;; So we have changed F(a+m,-a+n;c;z) to F(a',-a'+m+n;c;z).
3069
(cond ((and (minusp (mul n m))
3070
(greaterp (abs m) (abs n)))
3071
(return (list ap (sub ap n) m+n))))
3072
(return (list ap (add ap m) m+n))))
3074
;; F(a,b;c;z), where a + b is a (numerical) integer.
3076
;; If we're here, a and b are not integers. In general, a = s+f+m,
3077
;; where s is symbolic stuff and |f|<1 and m is an integer. We can
3078
;; also write b = -(s+f)+n.
3080
;; Let a' = s+f+m. Then a=a' and b = -a'+m+n, and we have converted
3081
;; F(a,b;c;z) to F(a',-a'+m+n;c;z), or, equivalently,
3082
;; F(-a'+m+n,a';c;z).
3084
(defun algii (a b c)
3085
(declare (ignore c))
3086
(prog (m n ap con sym m+n)
3087
;; We know a+b is an integer. In the most general form, we can
3088
;; have a = x+f+m and b = -(x-f)+n.
3090
;; Express a in the form sym + con.
3091
(cond ((not (setq sym (cdras 'f (s+c a))))
3094
;; a is of the form s+c. Look at the coefficient of s.
3095
;; If it's negative, swap a and b.
3096
(let ((res (m2 sym '((mtimes) ((coefft) (m $numberp)) ((coefft) (s nonnump)))
3098
(when (and res (minusp (cdras 'm res)))
3100
(setf sym (cdras 'f (s+c a)))))))
3101
(setq con (sub a sym))
3103
(setq m+n (add a b))
3104
;; Truncate con to an integer m.
3105
(setq m ($entier con))
3108
;; Make ap = s+con-m
3109
(setq ap (add (sub con m) ap))
3111
;; Return r*a+f, r*a+f+p, and m+n, where p is chosen to minimize
3112
;; the number of derivatives we need to take. Basically
3113
;; p=min(abs(m),abs(n)).
3115
;; Hmm, not sure why we do this. In any case, we need to do m+n
3116
;; differentiations. So the only simplification we could do is
3118
(cond ((and (minusp (mul n m))
3119
(greaterp (abs m) (abs n)))
3120
(return (list ap (sub ap n) m+n m n))))
3121
(return (list ap (add ap m) m+n m n))))
3123
;; We have something like F(s+m,-s+n;c;z)
3124
;; Rewrite it like F(a'+d,-a';c;z) where a'=s-n=-b and d=m+n.
3127
(let* ((sym (cdras 'f (s+c a)))
3128
(sign (cdras 'm (m2 sym '((mtimes) ((coefft) (m $numberp)) ((coefft) (s nonnump)))
3130
(when (and sign (minusp sign))
3132
(list nil (mul -1 b) (add a b))))
3135
;;Algor. 2F1-RL from thesis:step 4:dispatch on a+m,-a+n,1/2+l cases
3136
(defun step4 (a b c)
3137
;; F(a,b;c;z) where a+b is an integer and c+1/2 is an integer. If a
3138
;; and b are not integers themselves, we can derive the result from
3139
;; F(a1,-a1;1/2;z). However, if a and b are integers, we can't use
3140
;; that because F(a1,-a1;1/2;z) is a polynomial. We need to derive
3141
;; the result from F(1,1;3/2;z).
3142
(if (and (hyp-integerp a)
3147
(defun step4-a (a b c)
3148
(let* ((alglist (algii a b))
3149
(aprime (cadr alglist))
3152
($ratsimpexpons $true)
3154
;; At this point, we have F(a'+m,-a';1/2+n;z) where m and n are
3156
(cond ((hyp-integerp (add aprime (inv 2)))
3157
;; Ok. We have a problem if aprime + 1/2 is an integer.
3158
;; We can't always use the algorithm below because we have
3159
;; F(1/2,-1/2;1/2;z) which is 1F0(-1/2;;z) so the
3160
;; derivation isn't quite right. Also, sometimes we'll end
3161
;; up with a division by zero.
3163
;; Thus, We need to do something else. So, use A&S 15.3.3
3164
;; to change the problem:
3166
;; F(a,b;c;z) = (1-z)^(c-a-b)*F(c-a, c-b; c; z)
3170
;; F('a+m,-a';1/2+n;z) = (1-z)^(1/2+n-m)*F(1/2+n-a'-m,1/2+n+a';1/2+n;z)
3172
;; Recall that a' + 1/2 is an integer. Thus we have
3173
;; F(<int>,<int>,1/2+n;z), which we know how to handle in
3175
(gered1 (list a b) (list c) #'hgfsimp))
3178
(cond ((equal (checksigntm var) '$positive)
3179
(trig-log-1-pos aprime 'ell))
3180
((equal (checksigntm var) '$negative)
3181
(trig-log-1-neg (mul -1 aprime) aprime 'ell)))))
3182
;; Ok, this uses F(a,-a;1/2;z). Since there are 2 possible
3183
;; representations (A&S 15.1.11 and 15.1.17), we check the sign
3184
;; of the var (as done in trig-log-1) to select which form we
3185
;; want to use. The original didn't and seemed to want to use
3186
;; the negative form.
3188
;; With this change, F(a,-a;3/2;z) matches what A&S 15.2.6 would
3189
;; produce starting from F(a,-a;1/2;z), assuming z < 0.
3195
;; F(a,b;c;z), where a and b are (positive) integers and c = 1/2+l.
3196
;; This can be computed from F(1,1;3/2;z).
3198
;; Assume a < b, without loss of generality.
3200
;; F(m,n;3/2+L;z), m < n.
3202
;; We start from F(1,1;3/2;z). Use A&S 15.2.3, differentiating m
3203
;; times to get F(m,1;3/2;z). Swap a and b to get F(m,1;3/2;z) =
3204
;; F(1,m;3/2;z) and use A&S 15.2.3 again to get F(n,m;3/2;z) by
3205
;; differentiating n times. Finally, if L < 0, use A&S 15.2.4.
3206
;; Otherwise use A&S 15.2.6.
3208
;; I (rtoy) can't think of any way to do this with less than 3
3209
;; differentiations.
3211
;; Note that if d = (n-m)/2 is not an integer, we can write F(m,n;c;z)
3212
;; as F(-d+u,d+u;c;z) where u = (n+m)/2. In this case, we could use
3213
;; step4-a to compute the result.
3216
;; Transform F(a,b;c;z) to F(a+n,b;c;z), given F(a,b;c;z)
3217
(defun as-15.2.3 (a b c n arg fun)
3218
(declare (ignore b c))
3221
;; F(a+n,b;c;z) = z^(1-a)/poch(a,n)*diff(z^(a+n-1)*fun,z,n)
3222
(mul (inv (factf a n))
3223
(power arg (sub 1 a))
3224
($diff (mul (power arg (sub (add a n) 1))
3228
;; Transform F(a,b;c;z) to F(a,b;c-n;z), given F(a,b;c;z)
3229
(defun as-15.2.4 (a b c n arg fun)
3230
(declare (ignore a b))
3233
;; F(a,b;c-n;z) = 1/poch(c-n,n)/z^(c-n-1)*diff(z^(c-1)*fun,z,n)
3234
(mul (inv (factf (sub c n) n))
3235
(inv (power arg (sub (sub c n) 1)))
3236
($diff (mul (power arg (sub c 1))
3240
;; Transform F(a,b;c;z) to F(a-n,b;c;z), given F(a,b;c;z)
3241
(defun as-15.2.5 (a b c n arg fun)
3243
;; F(a-n,b;c;z) = 1/poch(c-a,n)*z^(1+a-c)*(1-z)^(c+n-a-b)
3244
;; *diff(z^(c-a+n-1)*(1-z)^(a+b-c)*F(a,b;c;z),z,n)
3246
(mul (inv (factf (sub c a) n))
3247
(power arg (sub (add a 1) c))
3249
(sub (add c n) (add a b)))
3250
($diff (mul (power arg (sub (add c n)
3257
;; Transform F(a,b;c;z) to F(a,b;c+n;z), given F(a,b;c;z)
3258
(defun as-15.2.6 (a b c n arg fun)
3260
;; F(a,b;c+n;z) = poch(c,n)/poch(c-a,n)/poch(c-b,n)*(1-z)^(c+n-a-b)
3261
;; *diff((1-z)^(a+b-c)*fun,z,n)
3264
(inv (factf (sub c a) n))
3265
(inv (factf (sub c b) n))
3266
(inv (power (sub 1 arg) (sub (add a b)
3268
($diff (mul (power (sub 1 arg) (sub (add a b) c))
3272
(defun step4-int (a b c)
3275
(let* ((s (gensym (symbol-name '#:step4-var-)))
3279
(res (cond ((eq (checksigntm var) '$negative)
3281
;; -%i*log(%i*sqrt(zn)+sqrt(1-zn))/(sqrt(1-zn)*sqrt(zn))
3283
(let ((root1-z (power (sub 1 s) (inv 2)))
3284
(rootz (power s (inv 2))))
3286
(mlog (add (mul '$%i rootz)
3291
;; F(1,1;3/2;z) = asin(sqrt(zp))/(sqrt(1-zp)*sqrt(zp))
3293
(let ((root1-z (power (sub 1 s) (inv 2)))
3294
(rootz (power s (inv 2))))
3298
;; Start with res = F(1,1;3/2;z). Compute F(m,1;3/2;z)
3299
(setf res (as-15.2.3 1 1 3//2 m s res))
3300
;; We now have res = C*F(m,1;3/2;z). Compute F(m,n;3/2;z)
3301
(setf res (as-15.2.3 1 a 3//2 n s res))
3302
;; We now have res = C*F(m,n;3/2;z). Now compute F(m,n;3/2+ell;z):
3305
(as-15.2.4 a b 3//2 (- ell) s res))
3307
(as-15.2.6 a b 3//2 ell s res)))))))
3309
;;Pattern match for s(ymbolic) + c(onstant) in parameter
3312
'((mplus) ((coeffpt)(f nonnump)) ((coeffpp)(c $numberp)))
1742
3315
(defun nonnump (z)
1743
(cond ((not ($numberp z)) t)
3316
(cond ((not ($numberp z)) t)
1746
;Algor. III from thesis:determines which Differ. Formula to use
3319
;;Algor. III from thesis:determines which Differ. Formula to use
1747
3320
(defun algiii (fun m n aprime)
1749
(setq mm (abs m) nn (abs n))
1750
(cond ((and (nni m) (nni n))
1751
(cond ((lessp m n) (return (f81 fun m n aprime)))
1752
(t (return (f85 fun mm nn aprime)))))
1753
((and (hyp-negp n) (hyp-negp m))
1754
(cond ((greaterp (abs n) (abs m))
1755
(return (f86 fun mm nn aprime)))
1756
(t (return (f82 fun mm nn aprime)))))
1757
((and (hyp-negp m) (nni n))(return (f83 fun mm nn aprime)))
1758
(t (return (f84 fun mm nn aprime))))))
3323
(cond ((and (nni m) (nni n))
3325
(f81 fun m n aprime))
3327
(f85 fun mm nn aprime))))
3328
((and (hyp-negp n) (hyp-negp m))
3329
(cond ((greaterp (abs m) (abs n))
3330
(f86 fun mm nn aprime))
3332
(f82 fun mm nn aprime))))
3333
((and (hyp-negp m) (nni n))
3334
(f83 fun mm nn aprime))
3336
(f84 fun mm nn aprime)))))
1760
;Factorial function:x*(x+1)*(x+2)...(x+n-1)
3338
;; Factorial function:x*(x+1)*(x+2)...(x+n-1)
3340
;; FIXME: This appears to be identical to fctrl above
1761
3341
(defun factf (x n)
1763
(t (mul x (factf (add x 1) (sub n 1))))))
1765
;Formula #85 from Yannis thesis:finds by differentiating F[2,1](a,b,c,z)
1766
; given F[2,1](a+m,b,c+n,z) where b=-a and c=1/2, n,m integers
3343
(t (mul x (factf (add x 1) (sub n 1))))))
3345
;;Formula #85 from Yannis thesis:finds by differentiating F[2,1](a,b,c,z)
3346
;; given F[2,1](a+m,b,c+n,z) where b=-a and c=1/2, n,m integers
3348
;; Like F81, except m > n.
3350
;; F(a+m,-a;c+n;z), m > n, c = 1/2, m and n are non-negative integers
3353
;; diff(z^(a+m-n-1)*F(a,-a;1/2;z),z,m-n) = poch(a,m-n)*z^(a-1)*F(a+m-n,-a;1/2;z)
3356
;; diff((1-z)^(a+m-1)*F(a+m-n,-a;1/2;z),z,n)
3357
;; = (-1)^n*poch(a+m-n,n)*poch(1/2+a,n)/poch(1/2,n)*(1-z)^(a+m-n-1)
3358
;; * F(a+m,-a;1/2+n;z)
1767
3360
(defun f85 (fun m n a)
1768
(mul (factf (inv 2) n)
1770
(inv (factf (sub (add a m) n) n))
1771
(inv (factf (sub (inv 2) (mul a -1)) n))
1772
(inv (factf a (- m n)))
1773
(power (sub 1 'ell) (sub (sub (add 1 n) m) a))
1775
(power (sub 1 'ell) (sub (add a m) 1))
1776
(power 'ell (sub 1 a))
1778
(power 'ell (sub (add a m -1) n))
1779
fun) 'ell (- m n))) 'ell n)))
1781
;Used to find negative things that are not integers,eg RAT's
1782
(defun hyp-negp(x) (cond ((equal (asksign x) '$negative) t)(t nil)))
3361
(mul (factf (inv 2) n)
3363
(inv (factf (sub (add a m)
3366
(inv (factf (sub (inv 2)
3369
(inv (factf a (- m n)))
3370
(power (sub 1 'ell) (sub (sub (add 1 n) m) a))
3371
($diff (mul (power (sub 1 'ell) (sub (add a m) 1))
3372
(power 'ell (sub 1 a))
3373
($diff (mul (power 'ell (sub (add a m -1) n))
3378
;;Used to find negative things that are not integers,eg RAT's
3380
(cond ((equal (asksign x) '$negative)
3384
;; F(a,-a+m; c+n; z) where m,n are non-negative integers, m < n, c = 1/2.
3387
;; diff((1-z)^(a+b-c)*F(a,b;c;z),z,n)
3388
;; = poch(c-a,n)*poch(c-b,n)/poch(c,n)*(1-z)^(a+b-c-n)*F(a,b;c+n;z)
3391
;; diff((1-z)^(a+m-1))*F(a,b;c;z),z,m)
3392
;; = (-1)^m*poch(a,m)*poch(c-b,m)/poch(c,m)*(1-z)^(a-1)*F(a+m,b;c+m;z)
3394
;; Rewrite F(a,-a+m; c+n;z) as F(-a+m, a; c+n; z). Then apply 15.2.6
3395
;; to F(-a,a;1/2;z), differentiating n-m times:
3397
;; diff((1-z)^(-1/2)*F(-a,a;1/2;z),z,n-m)
3398
;; = poch(1/2+a,n-m)*poch(1/2-a,n-m)/poch(1/2,n-m)*(1-z)^(-1/2-n+m)*F(-a,a;1/2+n-m;z)
3400
;; Multiply this result by (1-z)^(n-a-1/2) and apply 15.2.7, differentiating m times:
3402
;; diff((1-z)^(m-a-1)*F(-a,a;1/2+n-m;z),z,m)
3403
;; = (-1)^m*poch(-a,m)*poch(1/2+n-m-a,m)/poch(1/2+n-m)*(1-z)^(-a-1)*F(-a+m,a;1/2+n;z)
3405
;; Which gives F(-a+m,a;1/2+n;z), which is what we wanted.
1784
3406
(defun f81 (fun m n a)
1785
(mul (factf (add (inv 2) (- n m)) m)
1786
(factf (inv 2) (- n m))
1789
(inv (factf (add (inv 2) n (sub a m)) m))
1790
(inv (factf (sub (inv 2) a) (- n m)))
1791
(inv (factf (add (inv 2) a) (- n m)))
1792
(power (sub 1 'ell) (sub 1 a))
1794
(power (sub 1 'ell) (add a n (inv -2)))
1796
(power (sub 1 'ell) (inv -2))
1797
fun) 'ell (- n m))) 'ell m)))
1801
(mul (inv (factf (sub (inv 2) n) m))
1802
;; Was this both inverse?
1803
(inv (factf (sub (add (inv 2) m) n) (- n m)))
1804
(power 'ell (add n (inv 2)))
1805
(power (sub 1 'ell) (sub (add m (inv 2) a) n))
1806
($diff (mul (power (sub 1 'ell)
1807
(sub (sub n a) (inv 2)))
1808
($diff (mul (power 'ell (inv -2)) fun)
1816
(mul (factf (inv 2) n)
1817
(inv (factf (sub (inv 2) a) n))
1818
(inv (factf (add (sub (inv 2) a) n) m))
1819
(inv (factf (add (inv 2) a) n))
1820
(power (sub 1 'ell) (add m n (inv 2)))
1821
(power 'ell (sub (add (inv 2) a) n))
1822
($diff (mul (power 'ell (sub (sub (+ m n) a)(inv 2)))
1823
($diff (mul (power (sub 1 'ell)
1833
(mul (inv (mul (factf a m) (factf (sub (inv 2) n) n)))
1834
(power 'ell (sub 1 a))
1835
($diff (mul (power 'ell (sub (add a m n) (inv 2)))
1836
($diff (mul (power 'ell (inv -2)) fun)
1844
(mul (inv (mul (factf (sub (inv 2) n) n)
1845
(factf (sub (inv 2) a) (- m n))))
1846
(power 'ell (add n (inv 2)))
1847
(power (sub 1 'ell)(add (inv 2) a))
1848
($diff (mul (power 'ell a)
1849
(power (sub 1 'ell)(sub m a))
1850
($diff (mul (power 'ell
1851
(sub (sub (sub m n) (inv 2)) a))
1858
(eval-when (compile)
1859
(DECLARE-top (unSPECIAL serieslist VAR PAR ZEROSIGNTEST PRODUCTCASE
1860
fldeg flgkum listcmdiff checkcoefsignlist ))
1862
(declare-top (unspecial fun w b l alglist n c ))
b'\\ No newline at end of file'
3407
(mul (factf (add (inv 2) (- n m)) m)
3408
(factf (inv 2) (- n m))
3411
(inv (factf (add (inv 2) n (sub a m)) m))
3412
(inv (factf (sub (inv 2) a) (- n m)))
3413
(inv (factf (add (inv 2) a) (- n m)))
3414
(power (sub 1 'ell) (sub 1 a))
3415
($diff (mul (power (sub 1 'ell) (add a n (inv -2)))
3416
($diff (mul (power (sub 1 'ell) (inv -2))
3421
;; Like f86, but |n|>=|m|
3423
;; F(a-m,-a;1/2-n;z) where n >= m >0
3426
;; diff(z^(c-1)*F(a,b;c;z),z,n)
3427
;; = poch(c-n,n)*z^(c-n-1)*F(a;b;c-n;z)
3430
;; diff(z^(c-1)*(1-z)^(b-c+n)*F(a,b;c;z),z,n)
3431
;; = poch(c-n,n)*z^(c-n-1)*(1-z)^(b-c)*F(a-n,b;c-n;z)
3435
;; diff(z^(-1/2)*F(a,-a;1/2;z),z,n-m)
3436
;; = poch(1/2-n+m,n-m)*z^(m-n-1/2)*F(a,-a;1/2-n+m;z)
3438
;; diff(z^(m-n-1/2)*(1-z)^(n-a-1/2)*F(a,-a;1/2-n+m;z),z,m)
3439
;; = poch(1/2-n,m)*z^(-1/2-n)*(1-z)^(n-m-a-1/2)*F(a-m,-a;1/2-n;z)
3443
;; G(z) = z^(m-n-1/2)*F(a,-a;1/2-n+m;z)
3444
;; = z^(n-m+1/2)/poch(1/2-n+m,n-m)*diff(z^(-1/2)*F(a,-a;1/2;z),z,n-m)
3446
;; F(a-m,-a;1/2-n;z)
3447
;; = z^(n+1/2)*(1-z)^(m+a-1/2-n)/poch(1/2-n,m)*diff((1-z)^(n-a-1/2)*G(z),z,m)
3448
(defun f82 (fun m n a)
3449
(mul (inv (factf (sub (inv 2) n) m))
3450
(inv (factf (sub (add (inv 2) m) n) (- n m)))
3451
(power 'ell (add n (inv 2)))
3452
(power (sub 1 'ell) (sub (add m (inv 2) a) n))
3453
($diff (mul (power (sub 1 'ell)
3454
(sub (sub n a) (inv 2)))
3455
($diff (mul (power 'ell (inv -2)) fun)
3461
;; F(a+m,-a;1/2+n;z) with m,n integers and m < 0, n >= 0
3463
;; Write this more clearly as F(a-m,-a;1/2+n;z), m > 0, n >= 0
3464
;; or equivalently F(a-m,-a;c+n;z)
3467
;; diff((1-z)^(-1/2)*F(a,-a;1/2;z),z,n)
3468
;; = poch((1/2+a,n)*poch(1/2-a,n)/poch(1/2,n)*(1-z)^(-1/2-n)
3469
;; * F(a,-a;1/2+n;z)
3472
;; diff(z^(n+m-a-1/2)*(1-z)^(-1/2-n)*F(a,-a;1/2+n;z),z,m)
3473
;; = poch(1/2+n-a,m)*z^(1/2+n-a-1)*(1-z)^(-1/2-n-m)*F(a-m,-a;1/2+n;z)
3474
;; = poch(1/2+n-a,m)*z^(n-a-1/2)*(1-z)^(-1/2-n-m)*F(a-m,-a;1/2+n;z)
3476
;; (1-z)^(-1/2-n)*F(a,-a;1/2+n;z)
3477
;; = poch(1/2,n)/poch(1/2-a,n)/poch(1/2+a,n)*diff((1-z)^(-1/2)*F(a,-a;1/2;z),z,n)
3479
;; F(a-m,-a;1/2+n;z)
3480
;; = (1-z)^(n+m+1/2)*z^(a-n+1/2)/poch(1/2+n-a,m)
3481
;; *diff(z^(n+m-a-1/2)*(1-z)^(-1/2-n)*F(a,-a;1/2+n;z),z,m)
3482
(defun f83 (fun m n a)
3483
(mul (factf (inv 2) n)
3484
(inv (factf (sub (inv 2) a) n))
3485
(inv (factf (sub (add n (inv 2)) a) m))
3486
(inv (factf (add (inv 2) a) n))
3487
(power (sub 1 'ell) (add m n (inv 2)))
3488
(power 'ell (add (sub a n) (inv 2)))
3489
($diff (mul (power 'ell (sub (sub (+ m n) a) (inv 2)))
3490
($diff (mul (power (sub 1 'ell)
3498
;; The last case F(a+m,-a;c+n;z), m,n integers, m >= 0, n < 0
3500
;; F(a+m,-a;1/2-n;z)
3503
;; diff(z^(c-1)*F(a,b;c;z),z,n) = poch(c-n,n)*z^(c-n-1)*F(a,b;c-n;z)
3506
;; diff(z^(a+m-1)*F(a,b;c;z),z,m) = poch(a,n)*z^(a-1)*F(a+n,b;c;z)
3510
;; diff(z^(-1/2)*F(a,-a;1/2;z),z,n) = poch(1/2-n,n)*z^(-n-1/2)*F(a,-a;1/2-n;z)
3512
;; diff(z^(a+m-1)*F(a,-a;1/2-n;z),z,m) = poch(a,m)*z^(a-1)*F(a+m,-a;1/2-n;z)
3513
(defun f84 (fun m n a)
3514
(mul (inv (mul (factf a m)
3515
(factf (sub (inv 2) n) n)))
3516
(power 'ell (sub 1 a))
3517
($diff (mul (power 'ell (sub (add a m n) (inv 2)))
3518
($diff (mul (power 'ell (inv -2)) fun)
3524
;; Like f82, but |n|<|m|
3526
;; F(a-m,-a;1/2-n;z), 0 < n < m
3529
;; diff(z^(c-a+n-1)*(1-z)^(a+b-c)*F(a,b;c;z),z,n)
3530
;; = poch(c-a,n)*z^(c-a-1)*(1-z)^(a+b-c-n)*F(a-n,b;c;z)
3533
;; diff(z^(c-1)*(1-z)^(b-c+n)*F(a,b;c;z),z,n)
3534
;; = poch(c-n,n)*z^(c-n-1)*(1-z)^(b-c)*F(a-n,b;c-n;z)
3538
;; diff(z^(-a+m-n-1/2)*(1-z)^(-1/2)*F(a,-a;1/2;z),z,m-n)
3539
;; = poch(1/2-a,m-n)*z^(-a-1/2)*(1-z)^(-1/2-m+n)*F(a-m+n,-a;1/2;z)
3541
;; diff(z^(-1/2)*(1-z)^(-a-1/2+n)*F(a-m+n,-a;1/2;z),z,n)
3542
;; = poch(1/2-n,n)*z^(-n-1/2)*(1-z)^(-a-1/2)*F(a-m,-a;1/2-n;z)
3544
;; G(z) = z^(-a-1/2)*(1-z)^(-1/2-m+n)*F(a-m+n,-a;1/2;z)
3545
;; = 1/poch(1/2-a,m-n)*diff(z^(-a+m-n-1/2)*(1-z)^(-1/2)*F(a,-a;1/2;z),z,m-n)
3547
;; F(a-m,-a;1/2-n;z)
3548
;; = z^(n+1/2)*(1-z)^(a+1/2)/poch(1/2-n,n)
3549
;; *diff(z^(-1/2)*(1-z)^(-a-1/2+n)*F(a-m+n,-a;1/2;z),z,n)
3550
;; = z^(n+1/2)*(1-z)^(a+1/2)/poch(1/2-n,n)
3551
;; *diff(z^a*(1-z)^(m-a)*G(z),z,n)
3553
(defun f86 (fun m n a)
3554
(mul (inv (mul (factf (sub (inv 2) n) n)
3555
(factf (sub (inv 2) a) (- m n))))
3556
(power 'ell (add n (inv 2)))
3557
(power (sub 1 'ell) (add (inv 2) a))
3558
($diff (mul (power 'ell a)
3559
(power (sub 1 'ell) (sub m a))
3560
($diff (mul (power 'ell
3561
(sub (sub (sub m n) (inv 2)) a))
3568
;; F(-1/2+n, 1+m; 1/2+l; z)
3569
(defun hyp-atanh (a b c)
3570
;; We start with F(-1/2,1;1/2;z) = 1-sqrt(z)*atanh(sqrt(z)). We can
3571
;; derive the remaining forms by differentiating this enough times.
3573
;; FIXME: Do we need to assume z > 0? We do that anyway, here.
3574
(let* ((s (gensym (symbol-name '#:hyp-atanh-)))
3578
(res (sub 1 (mul (power s 1//2)
3579
`((%atanh) ,(power s 1//2)))))
3583
;; The total number of derivates we compute is n + m + ell. We
3584
;; should do something to reduce the number of derivatives.
3587
(format t "a ,b ,c = ~a ~a ~a~%" a b c)
3588
(format t "n, m, ell = ~a ~a ~a~%" n m ell)
3589
(format t "init a b c = ~a ~a ~a~%" new-a new-b new-c))
3591
(setf res (as-15.2.6 new-a new-b new-c ell s res))
3592
(setf new-c (add new-c ell)))
3594
(setf res (as-15.2.4 new-a new-b new-c (- ell) s res))
3595
(setf new-c (add new-c ell))))
3598
(maxima-display res)
3599
(format t "new a b c = ~a ~a ~a~%" new-a new-b new-c))
3602
(setf res (as-15.2.3 new-a new-b new-c n s res))
3603
(setf new-a (add new-a n)))
3605
(setf res (as-15.2.5 new-a new-b new-c (- n) s res))
3606
(setf new-a (add new-a n))))
3609
(format t "new a b c = ~a ~a ~a~%" new-a new-b new-c)
3610
(maxima-display res))
3612
(setf res (as-15.2.3 new-b new-a new-c m s res))
3613
(setf new-b (add new-b m)))
3615
(setf res (as-15.2.5 new-b new-a new-c (- m) s res))
3616
(setf new-b (add new-b m))))
3619
(format t "new a b c = ~a ~a ~a~%" new-a new-b new-c)
3620
(maxima-display res))
3625
#-gcl (:compile-toplevel)
3626
(declare-top (unspecial var *par* checkcoefsignlist))
added
removed
src/hyp.lisp
