8
8
;;; (c) copyright 1982 massachusetts institute of technology ;;;
9
9
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
12
12
(macsyma-module defint)
14
14
;;; this is the definite integration package.
15
; defint does definite integration by trying to find an
16
;appropriate method for the integral in question. the first thing that
17
;is looked at is the endpoints of the problem.
19
; i(grand,var,a,b) will be used for integrate(grand,var,a,b)
15
;; defint does definite integration by trying to find an
16
;;appropriate method for the integral in question. the first thing that
17
;;is looked at is the endpoints of the problem.
19
;; i(grand,var,a,b) will be used for integrate(grand,var,a,b)
21
;;references are to evaluation of definite integrals by symbolic
22
;;manipulation by paul s. wang.
21
;; References are to "Evaluation of Definite Integrals by Symbolic
22
;; Manipulation", by Paul S. Wang,
23
;; (http://www.lcs.mit.edu/publications/pubs/pdf/MIT-LCS-TR-092.pdf)
24
25
;; nointegrate is a macsyma level flag which inhibits indefinite
26
; abconv is a macsyma level flag which inhibits the absolute
29
; $defint is the top level function that takes the user input
30
;and does minor changes to make the integrand ready for the package.
32
; next comes defint, which is the function that does the
33
;integration. it is often called recursively from the bowels of the
34
;package. defint does some of the easy cases and dispatches to:
36
; dintegrate. this program first sees if the limits of
37
;integration are 0,inf or minf,inf. if so it sends the problem to
38
;ztoinf or mtoinf, respectivly.
39
; else, dintegrate tries:
41
; intsc1 - does integrals of sin's or cos's or exp(%i var)'s
42
; when the interval is 0,2 %pi or 0,%pi.
43
; method is conversion to rational function and find
44
; residues in the unit circle. [wang, pp 107-109]
46
; ratfnt - does rational functions over finite interval by
47
; doing polynomial part directly, and converting
48
; the rational part to an integral on 0,inf and finding
49
; the answer by residues.
51
; zto1 - i(x^(k-1)*(1-x)^(l-1),x,0,1) = beta(k,l) or
52
; i(log(x)*x^(x-1)*(1-x)^(l-1),x,0,1) = psi...
55
; dintrad- i(x^m/(a*x^2+b*x+c)^(n+3/2),x,0,inf) [wang, p 74]
57
; dintlog- i(log(g(x))*f(x),x,0,inf) = 0 (by symetry) or
58
; tries an integration by parts. (only routine to
59
; try integration by parts) [wang, pp 93-95]
61
; dintexp- i(f(exp(x)),x,a,b) = i(f(x+1),x,a',b')
63
;dintegrate also tries indefinite integration based on certain
64
;predicates (such as abconv) and tries breaking up the integrand
65
;over a sum or tries a change of variable.
67
; ztoinf is the routine for doing integrals over the range 0,inf.
68
; it goes over a series of routines and sees if any will work:
70
; scaxn - sc(b*x^n) (sc stands for sin or cos) [wang, pp 81-83]
72
; ssp - a*sc^n(r*x)/x^m [wang, pp 83,84]
74
; zmtorat- rational function. done by multiplication by plog(-x)
75
; and finding the residues over the keyhole contour
78
; log*rat- r(x)*log^n(x) [wang, pp 89-92]
80
; logquad0 log(x)/(a*x^2+b*x+c) uses formula
81
; i(log(x)/(x^2+2*x*a*cos(t)+a^2),x,0,inf) =
82
; t*log(a)/sin(t). a better formula might be
83
; i(log(x)/(x+b)/(x+c),x,0,inf) =
84
; (log^2(b)-log^2(c))/(2*(b-c))
86
; batapp - x^(p-1)/(b*x^n+a)^m uses formula related to the beta
87
; function [wang, p 71]
88
; there is also a special case when m=1 and a*b<0
91
; sinnu - x^-a*n(x)/d(x) [wang, pp 69-70]
93
; ggr - x^r*exp(a*x^n+b)
95
; dintexp- see dintegrate
97
; ztoinf also tries 1/2*mtoinf if the integrand is an even function
99
; mtoinf is the routine for doing integrals on minf,inf.
100
; it too tries a series of routines and sees if any succeed.
102
; scaxn - when the integrand is an even function, see ztoinf
104
; mtosc - exp(%i*m*x)*r(x) by residues on either the upper half
105
; plane or the lower half plane, depending on whether
106
; m is positive or negative.
108
; zmtorat- does rational function by finding residues in upper
111
; dintexp- see dintegrate
113
; rectzto%pi2 - poly(x)*rat(exp(x)) by finding residues in
114
; rectangle [wang, pp98-100]
116
; ggrm - x^r*exp((x+a)^n+b)
118
; mtoinf also tries 2*ztoinf if the integrand is an even function.
27
;; abconv is a macsyma level flag which inhibits the absolute
30
;; $defint is the top level function that takes the user input
31
;;and does minor changes to make the integrand ready for the package.
33
;; next comes defint, which is the function that does the
34
;;integration. it is often called recursively from the bowels of the
35
;;package. defint does some of the easy cases and dispatches to:
37
;; dintegrate. this program first sees if the limits of
38
;;integration are 0,inf or minf,inf. if so it sends the problem to
39
;;ztoinf or mtoinf, respectivly.
40
;; else, dintegrate tries:
42
;; intsc1 - does integrals of sin's or cos's or exp(%i var)'s
43
;; when the interval is 0,2 %pi or 0,%pi.
44
;; method is conversion to rational function and find
45
;; residues in the unit circle. [wang, pp 107-109]
47
;; ratfnt - does rational functions over finite interval by
48
;; doing polynomial part directly, and converting
49
;; the rational part to an integral on 0,inf and finding
50
;; the answer by residues.
52
;; zto1 - i(x^(k-1)*(1-x)^(l-1),x,0,1) = beta(k,l) or
53
;; i(log(x)*x^(x-1)*(1-x)^(l-1),x,0,1) = psi...
56
;; dintrad- i(x^m/(a*x^2+b*x+c)^(n+3/2),x,0,inf) [wang, p 74]
58
;; dintlog- i(log(g(x))*f(x),x,0,inf) = 0 (by symetry) or
59
;; tries an integration by parts. (only routine to
60
;; try integration by parts) [wang, pp 93-95]
62
;; dintexp- i(f(exp(k*x)),x,a,inf) = i(f(x+1)/(x+1),x,0,inf)
63
;; or i(f(x)/x,x,0,inf)/k. First case hold for a=0;
64
;; the second for a=minf. [wang 96-97]
66
;;dintegrate also tries indefinite integration based on certain
67
;;predicates (such as abconv) and tries breaking up the integrand
68
;;over a sum or tries a change of variable.
70
;; ztoinf is the routine for doing integrals over the range 0,inf.
71
;; it goes over a series of routines and sees if any will work:
73
;; scaxn - sc(b*x^n) (sc stands for sin or cos) [wang, pp 81-83]
75
;; ssp - a*sc^n(r*x)/x^m [wang, pp 86,87]
77
;; zmtorat- rational function. done by multiplication by plog(-x)
78
;; and finding the residues over the keyhole contour
81
;; log*rat- r(x)*log^n(x) [wang, pp 89-92]
83
;; logquad0 log(x)/(a*x^2+b*x+c) uses formula
84
;; i(log(x)/(x^2+2*x*a*cos(t)+a^2),x,0,inf) =
85
;; t*log(a)/sin(t). a better formula might be
86
;; i(log(x)/(x+b)/(x+c),x,0,inf) =
87
;; (log^2(b)-log^2(c))/(2*(b-c))
89
;; batapp - x^(p-1)/(b*x^n+a)^m uses formula related to the beta
90
;; function [wang, p 71]
91
;; there is also a special case when m=1 and a*b<0
94
;; sinnu - x^-a*n(x)/d(x) [wang, pp 69-70]
96
;; ggr - x^r*exp(a*x^n+b)
98
;; dintexp- see dintegrate
100
;; ztoinf also tries 1/2*mtoinf if the integrand is an even function
102
;; mtoinf is the routine for doing integrals on minf,inf.
103
;; it too tries a series of routines and sees if any succeed.
105
;; scaxn - when the integrand is an even function, see ztoinf
107
;; mtosc - exp(%i*m*x)*r(x) by residues on either the upper half
108
;; plane or the lower half plane, depending on whether
109
;; m is positive or negative.
111
;; zmtorat- does rational function by finding residues in upper
114
;; dintexp- see dintegrate
116
;; rectzto%pi2 - poly(x)*rat(exp(x)) by finding residues in
117
;; rectangle [wang, pp98-100]
119
;; ggrm - x^r*exp((x+a)^n+b)
121
;; mtoinf also tries 2*ztoinf if the integrand is an even function.
120
123
(load-macsyma-macros rzmac)
122
(DECLARE-TOP(*lexpr $DIFF $LIMIT $SUBSTITUTE $EZGCD $RATSIMP context)
123
(*expr subfunmake $coeff $logcontract $radcan $makegamma
124
$constantp $subvarp MAXIMA-SUBSTITUTE freeof ith
125
$oddp $hipow $multthru $xthru $num $denom
126
stripdollar MAXIMA-FIND sdiff partition
127
constant free mapatom
129
$ratdisrep ratdisrep $ratp ratp ratnumerator
130
sratsimp ratdenominator $ratsubst ratnump ratcoef
131
pterm rdis pdis ratrep newvar pdivide pointergp
133
$factor factor $sqfr oddelm zerop1
135
$asksign asksign $sign ask-integer assume forget
137
$residue residue res res1 polelist partnum
141
$rectform $realpart $imagpart trisplit cabs
143
among involve notinvolve
145
subin polyinx genfind xor fmt polyp numden andmapcar
146
abless1 even1 rddeg tansc radicalp deg simplerd
147
no-err-sub oscip %einvolve sin-sq-cos-sq-sub)
125
(declare-top(*lexpr $diff $limit $substitute $ezgcd $ratsimp context)
126
(*expr subfunmake $coeff $logcontract $radcan $makegamma
127
$constantp $subvarp maxima-substitute freeof ith
128
$oddp $hipow $multthru $xthru $num $denom
129
stripdollar maxima-find sdiff partition
130
constant free mapatom
132
$ratdisrep ratdisrep $ratp ratp ratnumerator
133
sratsimp ratdenominator $ratsubst ratnump ratcoef
134
pterm rdis pdis ratrep newvar pdivide pointergp
136
$factor factor $sqfr oddelm zerop1
138
$asksign asksign $sign ask-integer assume forget
140
$residue residue res res1 polelist partnum
144
$rectform $realpart $imagpart trisplit cabs
146
among involve notinvolve
148
subin polyinx genfind xor fmt polyp numden andmapcar
149
abless1 even1 rddeg tansc radicalp deg simplerd
150
no-err-sub oscip %einvolve sin-sq-cos-sq-sub)
149
152
;;;rsn* is in comdenom. does a ratsimp of numerator.
150
(SPECIAL *DEF2* PCPRNTD MTOINF* RSN* SEMIRAT*
151
SN* SD* LEADCOEF CHECKFACTORS
153
UL1 LL1 *DFLAG BPTU BPTD PLM* ZN ZD
154
*UPDN UL LL EXP PE* PL* RL* PL*1 RL*1
155
LOOPSTOP* VAR NN* ND* DN* p*
157
PLOGABS *ZEXPTSIMP? SCFLAG
158
sin-cos-recur rad-poly-recur dintlog-recur
159
dintexp-recur defintdebug defint-assumptions
161
global-defint-assumptions)
153
(special *def2* pcprntd *mtoinf* rsn* semirat*
154
sn* sd* leadcoef checkfactors
156
*ul1* *ll1* *dflag bptu bptd plm* zn
158
*updn ul ll exp pe* pl* rl* pl*1 rl*1
159
loopstop* var nn* nd* dn* p*
161
plogabs *zexptsimp? *scflag*
162
*sin-cos-recur* *rad-poly-recur* *dintlog-recur*
163
*dintexp-recur* defintdebug *defint-assumptions*
164
*current-assumptions*
165
*global-defint-assumptions*)
163
(ARRAY* (NOTYPE *I* 1 *J* 1))
167
(special $intanalysis $abconvtest $noprincipal $nointegrate)
169
(special $solveradcan $solvetrigwarn *roots *failures
170
$logabs $tlimswitch $maxposex $maxnegex
171
$trigsign $savefactors $radexpand $breakup $%emode
172
$float $exptsubst dosimp context rp-polylogp
173
%P%I HALF%PI %PI2 HALF%PI3 VARLIST genvar
174
$domain $m1pbranch errorsw errrjfflag raterr
176
;;LIMITP T Causes $ASKSIGN to do special things
177
;;For DEFINT like eliminate epsilon look for prin-inf
178
;;take realpart and imagpart.
180
;;If LIMITP is non-null ask-integer conses
181
;;its assumptions onto this list.
183
;If this switch is () then RPART returns ATAN's
167
(array* (notype *i* 1 *j* 1))
171
(special $intanalysis $abconvtest $noprincipal $nointegrate)
173
(special $solveradcan $solvetrigwarn *roots *failures
174
$logabs $tlimswitch $maxposex $maxnegex
175
$trigsign $savefactors $radexpand $breakup $%emode
176
$float $exptsubst dosimp context rp-polylogp
177
%p%i half%pi %pi2 half%pi3 varlist genvar
178
$domain $m1pbranch errorsw errrjfflag raterr
180
;;LIMITP T Causes $ASKSIGN to do special things
181
;;For DEFINT like eliminate epsilon look for prin-inf
182
;;take realpart and imagpart.
184
;;If LIMITP is non-null ask-integer conses
185
;;its assumptions onto this list.
187
;If this switch is () then RPART returns ATAN's
186
190
(declare-top(special infinities real-infinities infinitesimals))
187
(cond ;These are really defined in LIMIT but DEFINT uses them also.
191
(cond ;These are really defined in LIMIT but DEFINT uses them also.
188
192
((not (boundp 'infinities))
189
(setq INFINITIES '($INF $MINF $INFINITY))
190
(setq REAL-INFINITIES '($INF $MINF))
191
(setq INFINITESIMALS '($ZEROA $ZEROB))))
193
(setq infinities '($inf $minf $infinity))
194
(setq real-infinities '($inf $minf))
195
(setq infinitesimals '($zeroa $zerob))))
193
197
(defmvar defintdebug () "If true Defint prints out debugging information")
195
(DEFMVAR INTEGERL NIL
196
"An integer-list for non-atoms found out to be INTEGERs")
198
(DEFMVAR NONINTEGERL NIL
199
"A non-integer-list for non-atoms found out to be NONINTEGERs")
201
(DEFUN $DEFINT (EXP VAR LL UL)
202
(let ((global-defint-assumptions ())
203
(integer-info ()) (integerl integerl) (nonintegerl nonintegerl))
204
(with-new-context (context)
206
(let ((defint-assumptions ()) (*def2* ()) (rad-poly-recur ())
207
(sin-cos-recur ()) (dintexp-recur ()) (dintlog-recur 0.)
208
(ans nil) (orig-exp exp) (orig-var var)
209
(orig-ll ll) (orig-ul ul)
210
(pcprntd nil) (*NODIVerg nil) ($logabs t) (limitp t)
212
($domain '$real) ($m1pbranch ())) ;Try this out.
214
(FIND-FUNCTION '$LIMIT)
215
(make-global-assumptions) ;sets global-defint-assumptions
216
(FIND-FUNCTION '$RESIDUE)
217
(SETQ EXP (RATDISREP EXP))
218
(SETQ VAR (RATDISREP VAR))
219
(SETQ LL (RATDISREP LL))
220
(SETQ UL (RATDISREP UL))
221
(COND (($CONSTANTP VAR)
222
(merror "Variable of integration not a variable: ~M"
224
(($SUBVARP VAR) (SETQ VAR (STRIPDOLLAR (CAAR VAR)))
225
(SETQ EXP ($SUBSTITUTE VAR orig-VAR EXP))))
226
(COND ((NOT (ATOM VAR))
227
(merror "Improper variable of integration: ~M" VAR))
230
(SETQ VAR (STRIPDOLLAR VAR))
231
(SETQ EXP ($SUBSTITUTE VAR orig-VAR EXP))))
232
(cond ((not (equal ($ratsimp ($imagpart ll)) 0))
233
(merror "Defint: Lower limit of integration must be real."))
234
((not (equal ($ratsimp ($imagpart ul)) 0))
236
"Defint: Upper limit of integration must be real.")))
238
(COND ((SETQ ANS (DEFINT EXP VAR LL UL))
239
(SETQ ANS (SUBST orig-VAR VAR ANS))
240
(COND ((atom ans) ans)
241
((and (free ans '%limit)
242
(free ans '%integrate)
243
(OR (not (free ans '$INF))
244
(not (free ans '$MINF))
245
(not (free ans '$INFINITY))))
247
((not (free ans '$und))
248
`((%integrate) ,orig-exp ,orig-var ,orig-ll ,orig-ul))
250
(t `((%integrate) ,orig-exp ,orig-var ,orig-ll ,orig-ul))))
251
(forget-global-assumptions)))))
253
(DEFUN EEZZ (EXP LL UL)
254
(COND ((OR (polyinx EXP VAR nil)
255
(CATCH 'PIN%EX (PIN%EX EXP)))
256
(SETQ EXP (ANTIDERIV EXP))
257
;;;If antideriv can't do it, returns nil
258
;;;use limit to evaluate every answer returned by antideriv.
259
(COND ((NULL EXP) NIL)
260
(t (INTSUBS EXP LL UL))))))
199
(defmvar integerl nil
200
"An integer-list for non-atoms found out to be `integer's")
202
(defmvar nonintegerl nil
203
"A non-integer-list for non-atoms found out to be `noninteger's")
205
(defun $defint (exp var ll ul)
206
(let ((*global-defint-assumptions* ())
207
(integer-info ()) (integerl integerl) (nonintegerl nonintegerl))
208
(with-new-context (context)
210
(let ((*defint-assumptions* ()) (*def2* ()) (*rad-poly-recur* ())
211
(*sin-cos-recur* ()) (*dintexp-recur* ()) (*dintlog-recur* 0.)
212
(ans nil) (orig-exp exp) (orig-var var)
213
(orig-ll ll) (orig-ul ul)
214
(pcprntd nil) (*nodiverg nil) ($logabs t) (limitp t)
216
($domain '$real) ($m1pbranch ())) ;Try this out.
218
(find-function '$limit)
219
(make-global-assumptions) ;sets *global-defint-assumptions*
220
(find-function '$residue)
221
(setq exp (ratdisrep exp))
222
(setq var (ratdisrep var))
223
(setq ll (ratdisrep ll))
224
(setq ul (ratdisrep ul))
225
(cond (($constantp var)
226
(merror "Variable of integration not a variable: ~M"
228
(($subvarp var) (setq var (stripdollar (caar var)))
229
(setq exp ($substitute var orig-var exp))))
230
(cond ((not (atom var))
231
(merror "Improper variable of integration: ~M" var))
234
(setq var (stripdollar var))
235
(setq exp ($substitute var orig-var exp))))
236
(cond ((not (equal ($ratsimp ($imagpart ll)) 0))
237
(merror "Defint: Lower limit of integration must be real."))
238
((not (equal ($ratsimp ($imagpart ul)) 0))
240
"Defint: Upper limit of integration must be real.")))
242
(cond ((setq ans (defint exp var ll ul))
243
(setq ans (subst orig-var var ans))
244
(cond ((atom ans) ans)
245
((and (free ans '%limit)
246
(free ans '%integrate)
247
(or (not (free ans '$inf))
248
(not (free ans '$minf))
249
(not (free ans '$infinity))))
251
((not (free ans '$und))
252
`((%integrate) ,orig-exp ,orig-var ,orig-ll ,orig-ul))
254
(t `((%integrate) ,orig-exp ,orig-var ,orig-ll ,orig-ul))))
255
(forget-global-assumptions)))))
257
(defun eezz (exp ll ul)
258
(cond ((or (polyinx exp var nil)
259
(catch 'pin%ex (pin%ex exp)))
260
(setq exp (antideriv exp))
261
;; If antideriv can't do it, returns nil
262
;; use limit to evaluate every answer returned by antideriv.
263
(cond ((null exp) nil)
264
(t (intsubs exp ll ul))))))
261
265
;;;Hack the expression up for exponentials.
263
267
(defun sinintp (expr var)
264
;;; Is this expr a candidate for SININT ?
265
(let ((expr (factor expr))
268
(setq numer ($num expr))
269
(setq denom ($denom expr))
270
(cond ((polyinx numer var nil)
271
(cond ((and (polyinx denom var nil)
272
(deg-lessp denom var 2))
274
;;;ERF type things go here.
275
((let ((exponent (%einvolve numer)))
276
(and (polyinx exponent var nil)
277
(deg-lessp exponent var 2)))
278
(cond ((free denom var)
268
;; Is this expr a candidate for SININT ?
269
(let ((expr (factor expr))
272
(setq numer ($num expr))
273
(setq denom ($denom expr))
274
(cond ((polyinx numer var nil)
275
(cond ((and (polyinx denom var nil)
276
(deg-lessp denom var 2))
278
;;ERF type things go here.
279
((let ((exponent (%einvolve numer)))
280
(and (polyinx exponent var nil)
281
(deg-lessp exponent var 2)))
282
(cond ((free denom var)
281
285
(defun deg-lessp (expr var power)
282
(cond ((or (atom expr)
286
(do ((ops (cdr expr) (cdr ops)))
288
(cond ((not (deg-lessp (car ops) var power))
291
(and (or (not (alike1 (cadr expr) var))
292
(and (numberp (caddr expr))
293
(not (eq (asksign (m+ power (m- (caddr expr))))
295
(deg-lessp (cadr expr) var power)))))
286
(cond ((or (atom expr)
290
(do ((ops (cdr expr) (cdr ops)))
292
(cond ((not (deg-lessp (car ops) var power))
295
(and (or (not (alike1 (cadr expr) var))
296
(and (numberp (caddr expr))
297
(not (eq (asksign (m+ power (m- (caddr expr))))
299
(deg-lessp (cadr expr) var power)))
301
(member 'array (car expr))
302
(not (eq var (caar expr))))
303
;; We have some subscripted variable that's not our variable
304
;; (I think), so it's deg-lessp.
306
;; FIXME: Is this the best way to handle this? Are there
307
;; other cases we're mising here?
301
(setq ans (SININT A VAR))
302
(COND ((AMONG '%INTEGRATE Ans) NIL)
303
(T (SIMPLIFY Ans)))))
314
(setq ans (sinint a var))
315
(cond ((among '%integrate ans) nil)
316
(t (simplify ans)))))
305
318
;;;This routine tries to take a limit a couple of ways.
306
319
(defmfun get-limit nargs
310
323
(t (cond ((and (or (equal nargs 3) (equal nargs 4))
311
324
(memq (setq val (arg 3)) '($inf $minf)))
312
325
(setq var (arg 2))
313
(setq exp (MAXIMA-SUBSTITUTE (m^t var -1) var (arg 1)))
326
(setq exp (maxima-substitute (m^t var -1) var (arg 1)))
314
327
(cond ((eq val '$inf) (setq dir '$plus))
315
328
(t (setq dir '$minus)))
316
;(setq ans (apply 'limit-no-err `(,exp ,var 0 ,dir)))
329
;(setq ans (apply 'limit-no-err `(,exp ,var 0 ,dir)))
317
330
(setq ans (limit-no-err exp var 0 dir))
318
331
(cond ((not (among '%limit ans)) ans)
321
;(defun limit-no-err nargs
322
; (let ((errorsw t) (ans ()))
323
; (setq ans (catch 'errorsw (apply '$limit (listify nargs))))
324
; (cond ((not (eq ans t)) ans)
334
;;(defun limit-no-err nargs
335
;; (let ((errorsw t) (ans ()))
336
;; (setq ans (catch 'errorsw (apply '$limit (listify nargs))))
337
;; (cond ((not (eq ans t)) ans)
326
339
(defun limit-no-err (&rest argvec)
327
340
(declare (special errorsw))
328
341
(let ((errorsw t) (ans nil))
329
342
(setq ans (catch 'errorsw (apply '$limit argvec)))
330
343
(if (eq ans t) nil ans)))
332
(DEFUN INTCV (NV IND FLAG)
333
(let ((D (bx**n+a nv))
334
(*ROOTS ()) (*FAILURES ()) ($BREAKUP ()))
335
(cond ((AND (EQ UL '$INF)
337
(EQUAL (CADR D) 1)) ())
338
(t (SOLVE (m+t 'YX (m*t -1. NV)) VAR 1.)
339
(COND (*ROOTS (SETQ D (SUBST VAR 'YX (CADDAR *ROOTS)))
340
(COND (FLAG (INTCV2 D IND NV))
341
(T (INTCV1 D IND NV))))
346
(defun intcv (nv ind flag)
347
(let ((d (bx**n+a nv))
348
(*roots ()) (*failures ()) ($breakup ()))
349
(cond ((and (eq ul '$inf)
351
(equal (cadr d) 1)) ())
352
(t (solve (m+t 'yx (m*t -1. nv)) var 1.)
353
(format t "*roots = ~A~%" *roots)
354
(format t "subst ~A~%" (caddar *roots))
356
(setq d (subst var 'yx (caddar *roots)))
357
(format t "d = ~A~%" d)
358
(cond (flag (intcv2 d ind nv))
359
(t (intcv1 d ind nv))))
344
(DEFUN INTCV1 (D IND NV)
345
(COND ((AND (INTCV2 D IND NV)
346
(NOT (alike1 LL1 UL1)))
348
(DEFINT EXP1 VAR LL1 UL1)))))
350
(DEFUN INTCV2 (D IND NV)
352
(AND (COND ((AND (ZEROP1 (m+ LL UL))
354
(SETQ EXP1 (m* 2 EXP1)
355
LL1 (LIMCP NV VAR 0 '$PLUS)))
356
(T (SETQ LL1 (LIMCP NV VAR LL '$PLUS))))
357
(SETQ UL1 (LIMCP NV VAR UL '$MINUS))))
359
(DEFUN LIMCP (A B C D)
360
(let ((Ans ($LIMIT A B C D)))
361
(COND ((NOT (OR (null ans)
367
(DEFUN INTCV3 (D IND NV)
368
(SETQ NN* ($RATSIMP (SDIFF D VAR)))
369
(SETQ EXP1 (SUBST 'YX NV EXP))
370
(SETQ EXP1 (m* NN* (COND (IND EXP)
371
(T (SUBST D VAR EXP1)))))
372
(SETQ EXP1 (sRATSIMP (SUBST VAR 'YX EXP1))))
374
(DEFUN DEFINT (EXP VAR LL UL)
375
(let ((old-assumptions defint-assumptions) (current-assumptions ()))
378
(setq current-assumptions (make-defint-assumptions 'noask))
379
(let ((exp (resimplify exp))
380
(var (resimplify var))
383
;; D (not used? -- cwh)
384
ANS NN* DN* ND* $NOPRINCIPAL)
385
(cond ((setq ans (defint-list exp var ll ul))
390
((NOT (AMONG VAR EXP))
391
(COND ((OR (MEMQ UL '($INF $MINF))
392
(MEMQ LL '($INF $MINF)))
394
(T (SETQ ANS (m* EXP (m+ UL (m- LL))))
396
(let* ((EXP (RMCONST1 EXP))
398
(EXP (%i-out-of-denom (CDR EXP))))
399
(COND ((AND (NOT $NOINTEGRATE)
401
(or (among 'mqapply exp)
402
(NOT (MEMQ (CAAR EXP)
403
'(MEXPT MPLUS MTIMES %SIN %COS
404
%TAN %SINH %COSH %TANH
405
%LOG %ASIN %ACOS %ATAN
409
(COND ((SETQ ANS (ANTIDERIV EXP))
410
(SETQ ANS (INTSUBS ANS LL UL))
413
(SETQ EXP (TANSC EXP))
414
(cond ((setq ans (initial-analysis exp var ll ul))
415
(return (m* c ans))))
417
(restore-defint-assumptions old-assumptions current-assumptions))))
362
(defun intcv (nv ind flag)
363
(let ((d (bx**n+a nv))
364
(*roots ()) (*failures ()) ($breakup ()))
365
(cond ((and (eq ul '$inf)
367
(equal (cadr d) 1)) ())
369
;; This is a hack! If nv is of the form b*x^n+a, we can
370
;; solve the equation manually instead of using solve.
371
;; Why? Because solve asks us for the sign of yx and
374
;; Solve yx = b*x^n+a, for x. Any root will do. So we
375
;; have x = ((yx-a)/b)^(1/n).
376
(destructuring-bind (a n b)
378
(let ((root (power* (div (sub 'yx a) b) (inv n))))
380
(setq d (subst var 'yx root))
381
(cond (flag (intcv2 d ind nv))
382
(t (intcv1 d ind nv))))
385
(solve (m+t 'yx (m*t -1. nv)) var 1.)
387
(setq d (subst var 'yx (caddar *roots)))
388
(cond (flag (intcv2 d ind nv))
389
(t (intcv1 d ind nv))))
392
(defun intcv1 (d ind nv)
393
(cond ((and (intcv2 d ind nv)
394
(not (alike1 *ll1* *ul1*)))
396
(defint exp1 var *ll1* *ul1*)))))
398
(defun intcv2 (d ind nv)
400
(and (cond ((and (zerop1 (m+ ll ul))
402
(setq exp1 (m* 2 exp1)
403
*ll1* (limcp nv var 0 '$plus)))
404
(t (setq *ll1* (limcp nv var ll '$plus))))
405
(setq *ul1* (limcp nv var ul '$minus))))
407
(defun limcp (a b c d)
408
(let ((ans ($limit a b c d)))
409
(cond ((not (or (null ans)
415
(defun intcv3 (d ind nv)
416
(setq nn* ($ratsimp (sdiff d var)))
417
(setq exp1 (subst 'yx nv exp))
418
(setq exp1 (m* nn* (cond (ind exp)
419
(t (subst d var exp1)))))
420
(setq exp1 (sratsimp (subst var 'yx exp1))))
422
(defun defint (exp var ll ul)
423
(let ((old-assumptions *defint-assumptions*) (*current-assumptions* ()))
426
(setq *current-assumptions* (make-defint-assumptions 'noask))
427
(let ((exp (resimplify exp))
428
(var (resimplify var))
431
;; D (not used? -- cwh)
432
ans nn* dn* nd* $noprincipal)
433
(cond ((setq ans (defint-list exp var ll ul))
438
((not (among var exp))
439
(cond ((or (memq ul '($inf $minf))
440
(memq ll '($inf $minf)))
442
(t (setq ans (m* exp (m+ ul (m- ll))))
444
(let* ((exp (rmconst1 exp))
446
(exp (%i-out-of-denom (cdr exp))))
447
(cond ((and (not $nointegrate)
449
(or (among 'mqapply exp)
450
(not (memq (caar exp)
451
'(mexpt mplus mtimes %sin %cos
452
%tan %sinh %cosh %tanh
453
%log %asin %acos %atan
457
(cond ((setq ans (antideriv exp))
458
(setq ans (intsubs ans ll ul))
461
(setq exp (tansc exp))
462
(cond ((setq ans (initial-analysis exp var ll ul))
463
(return (m* c ans))))
465
(restore-defint-assumptions old-assumptions *current-assumptions*))))
419
467
(defun defint-list (exp var ll ul)
420
(COND ((AND (NOT (ATOM EXP))
422
'(MEQUAL MLIST $MATRIX)))
423
(let ((ANS (CONS (CAR EXP)
426
(DEFINT SUB-EXP VAR LL UL))
428
(COND (ANS (simplify ANS))
468
(cond ((and (not (atom exp))
470
'(mequal mlist $matrix)))
471
(let ((ans (cons (car exp)
474
(defint sub-exp var ll ul))
476
(cond (ans (simplify ans))
432
480
(defun initial-analysis (exp var ll ul)
433
(let ((pole (cond ((not $intanalysis)
434
'$no) ;don't do any checking.
435
(t (poles-in-interval exp var ll ul)))))
436
(cond ((eq pole '$no)
437
(cond ((and (oddfn exp var)
438
(or (and (eq ll '$minf)
440
(eq ($sign (m+ ll ul))
442
(t (parse-integrand exp var ll ul))))
443
((eq pole '$unknown) ())
444
(t (principal-value-integral exp var ll ul pole)))))
481
(let ((pole (cond ((not $intanalysis)
482
'$no) ;don't do any checking.
483
(t (poles-in-interval exp var ll ul)))))
484
(cond ((eq pole '$no)
485
(cond ((and (oddfn exp var)
486
(or (and (eq ll '$minf)
488
(eq ($sign (m+ ll ul))
490
(t (parse-integrand exp var ll ul))))
491
((eq pole '$unknown) ())
492
(t (principal-value-integral exp var ll ul pole)))))
446
494
(defun parse-integrand (exp var ll ul)
448
(COND ((SETQ ANS (EEZZ EXP LL UL)) ans)
450
(setq ans (method-by-limits EXP VAR LL UL))) ans)
452
(setq ans (INTBYTERM EXP T))) ans)
453
((setq ans (method-by-limits EXP VAR LL UL)) ans)
496
(cond ((setq ans (eezz exp ll ul)) ans)
498
(setq ans (method-by-limits exp var ll ul))) ans)
500
(setq ans (intbyterm exp t))) ans)
501
((setq ans (method-by-limits exp var ll ul)) ans)
462
510
(defun method-by-limits (exp var ll ul)
463
(let ((old-assumptions defint-assumptions))
464
(setq current-assumptions (make-defint-assumptions 'noask))
465
;;Should be a PROG inside of unwind-protect, but Multics has a compiler
466
;;bug wrt. and I want to test this code now.
468
(COND ((and (and (EQ UL '$INF)
469
(eq ll '$minf)) (mtoinf exp var)))
470
((and (and (eq ul '$inf)
471
(equal ll 0.)) (ztoinf exp var)))
511
(let ((old-assumptions *defint-assumptions*))
512
(setq *current-assumptions* (make-defint-assumptions 'noask))
513
;;Should be a PROG inside of unwind-protect, but Multics has a compiler
514
;;bug wrt. and I want to test this code now.
516
(cond ((and (and (eq ul '$inf)
519
((and (and (eq ul '$inf)
472
522
;;;This seems((and (and (eq ul '$inf)
473
523
;;;fairly losing (setq exp (subin (m+ ll var) exp))
474
524
;;; (setq ll 0.))
475
525
;;; (ztoinf exp var)))
477
(freeof var ul) (eq ($asksign ul) '$pos) (zto1 exp)))
478
; ((and (and (equal ul 1.)
479
; (equal ll 0.)) (zto1 exp)))
480
(t (dintegrate exp var ll ul)))
481
(restore-defint-assumptions old-assumptions defint-assumptions))))
528
(eq ($asksign ul) '$pos)
530
;; ((and (and (equal ul 1.)
531
;; (equal ll 0.)) (zto1 exp)))
532
(t (dintegrate exp var ll ul)))
533
(restore-defint-assumptions old-assumptions *defint-assumptions*))))
484
(DEFUN DINTEGRATE (EXP VAR LL UL)
485
(let ((ans nil) (arg nil) (scflag nil)
486
(*dflag nil) ($%emode t))
536
(defun dintegrate (exp var ll ul)
537
(let ((ans nil) (arg nil) (*scflag* nil)
538
(*dflag nil) ($%emode t))
487
539
;;;NOT COMPLETE for sin's and cos's.
488
(cond ((and (not sin-cos-recur)
492
((and (not rad-poly-recur)
493
(notinvolve exp '(%log))
494
(not (%einvolve exp))
495
(method-radical-poly exp var ll ul)))
496
((and (not (equal dintlog-recur 2.))
497
(SETQ arg (INVOLVE EXP '(%log)))
499
((and (not dintexp-recur)
500
(SETQ arg (%EINVOLVE exp))
502
((and (NOT (RATP EXP VAR))
503
(SETQ ANS ($EXPAND EXP))
504
(NOT (ALIKE1 ANS EXP))
506
((setq ans (antideriv exp))
540
(cond ((and (not *sin-cos-recur*)
544
((and (not *rad-poly-recur*)
545
(notinvolve exp '(%log))
546
(not (%einvolve exp))
547
(method-radical-poly exp var ll ul)))
548
((and (not (equal *dintlog-recur* 2.))
549
(setq arg (involve exp '(%log)))
551
((and (not *dintexp-recur*)
552
(setq arg (%einvolve exp))
554
((and (not (ratp exp var))
555
(setq ans ($expand exp))
556
(not (alike1 ans exp))
558
((setq ans (antideriv exp))
510
562
(defun method-radical-poly (exp var ll ul)
511
563
;;;Recursion stopper
512
(let ((rad-poly-recur t) ;recursion stopper
564
(let ((*rad-poly-recur* t) ;recursion stopper
514
(cond ((and (sinintp EXP VAR)
566
(cond ((and (sinintp exp var)
515
567
(setq result (antideriv exp))
516
568
(intsubs result ll ul)))
518
(setq result (RATFNT EXP))))
570
(setq result (ratfnt exp))))
519
571
((and (setq result (antideriv exp))
520
572
(intsubs result ll ul)))
525
(setq result (CV EXP))))
577
(setq result (cv exp))))
528
580
;;; LIMIT loss can't set logabs to true for these cases.
651
732
((and (eq ul '$minf)
652
733
(not (eq ll '$minf))) '$neg)
653
734
(t ($sign ($limit (m+ ul (m- ll))))))))
654
(cond ((eq sign-ul-ll '$pos)
655
(setq defint-assumptions
656
`(,(assume `((mgreaterp) ,var ,ll))
657
,(assume `((mgreaterp) ,ul ,var)))))
658
((eq sign-ul-ll '$neg)
659
(setq defint-assumptions
660
`(,(assume `((mgreaterp) ,var ,ul))
661
,(assume `((mgreaterp) ,ll ,var))))))
662
(cond ((and (eq sign-ll '$pos)
664
(setq defint-assumptions
665
`(,(assume `((mgreaterp) ,var 0))
666
,@defint-assumptions)))
667
((and (eq sign-ll '$neg)
669
(setq defint-assumptions
670
`(,(assume `((mgreaterp) 0 ,var))
671
,@defint-assumptions)))
672
(t defint-assumptions))))))
735
(cond ((eq sign-ul-ll '$pos)
736
(setq *defint-assumptions*
737
`(,(assume `((mgreaterp) ,var ,ll))
738
,(assume `((mgreaterp) ,ul ,var)))))
739
((eq sign-ul-ll '$neg)
740
(setq *defint-assumptions*
741
`(,(assume `((mgreaterp) ,var ,ul))
742
,(assume `((mgreaterp) ,ll ,var))))))
743
(cond ((and (eq sign-ll '$pos)
745
(setq *defint-assumptions*
746
`(,(assume `((mgreaterp) ,var 0))
747
,@*defint-assumptions*)))
748
((and (eq sign-ll '$neg)
750
(setq *defint-assumptions*
751
`(,(assume `((mgreaterp) 0 ,var))
752
,@*defint-assumptions*)))
753
(t *defint-assumptions*))))))
674
755
(defun restore-defint-assumptions (old-assumptions assumptions)
675
756
(do ((llist assumptions (cdr llist)))
677
(forget (car llist)))
758
(forget (car llist)))
678
759
(do ((llist old-assumptions (cdr llist)))
680
(assume (car llist))))
761
(assume (car llist))))
682
763
(defun make-global-assumptions ()
683
(setq global-defint-assumptions
684
(cons (assume '((mgreaterp) *z* 0.))
685
global-defint-assumptions))
686
;;; *Z* is a "zero parameter" for this package.
687
;;; Its also used to transform.
688
;; limit(exp,var,val,dir) -- limit(exp,tvar,0,dir)
689
(setq global-defint-assumptions
690
(cons (assume '((mgreaterp) epsilon 0.))
691
global-defint-assumptions))
692
(setq global-defint-assumptions
693
(cons (assume '((mlessp) epsilon 1.0e-8))
694
global-defint-assumptions))
695
;;; EPSILON is used in principal vaule code to denote the familiar
696
;;; mathematical entity.
697
(setq global-defint-assumptions
698
(cons (assume '((mgreaterp) prin-inf 1.0e+8))
699
global-defint-assumptions)))
764
(setq *global-defint-assumptions*
765
(cons (assume '((mgreaterp) *z* 0.))
766
*global-defint-assumptions*))
767
;; *Z* is a "zero parameter" for this package.
768
;; Its also used to transform.
769
;; limit(exp,var,val,dir) -- limit(exp,tvar,0,dir)
770
(setq *global-defint-assumptions*
771
(cons (assume '((mgreaterp) epsilon 0.))
772
*global-defint-assumptions*))
773
(setq *global-defint-assumptions*
774
(cons (assume '((mlessp) epsilon 1.0e-8))
775
*global-defint-assumptions*))
776
;; EPSILON is used in principal vaule code to denote the familiar
777
;; mathematical entity.
778
(setq *global-defint-assumptions*
779
(cons (assume '((mgreaterp) prin-inf 1.0e+8))
780
*global-defint-assumptions*)))
700
782
;;; PRIN-INF Is a special symbol in the principal value code used to
701
783
;;; denote an end-point which is proceeding to infinity.
703
785
(defun forget-global-assumptions ()
704
(do ((llist global-defint-assumptions (cdr llist)))
706
(forget (car llist)))
707
(cond ((not (null integer-info))
708
(do ((llist integer-info (cdr llist)))
710
(I-$remove `(,(cadar llist) ,(caddar llist)))))))
712
(DEFUN order-limits (ask-or-not)
713
(cond ((or (not (equal ($imagpart ll) 0))
714
(not (equal ($imagpart ul) 0))) ())
715
(t (COND ((ALIKE1 LL (m*t -1 '$INF))
717
(COND ((ALIKE1 UL (m*t -1 '$INF))
719
(cond ((alike1 ll (m*t -1 '$minf))
721
(cond ((alike1 ul (m*t -1 '$minf))
723
(COND ((EQ UL '$INF) NIL)
725
(SETQ EXP (SUBIN (m- VAR) EXP))
731
(setq exp (m- exp))))
732
;;;Fix limits so that ll < ul.
733
(let ((D (COMPLM ask-or-not)))
741
(DEFUN COMPLM (ask-or-not)
786
(do ((llist *global-defint-assumptions* (cdr llist)))
788
(forget (car llist)))
789
(cond ((not (null integer-info))
790
(do ((llist integer-info (cdr llist)))
792
(i-$remove `(,(cadar llist) ,(caddar llist)))))))
794
(defun order-limits (ask-or-not)
795
(cond ((or (not (equal ($imagpart ll) 0))
796
(not (equal ($imagpart ul) 0))) ())
797
(t (cond ((alike1 ll (m*t -1 '$inf))
799
(cond ((alike1 ul (m*t -1 '$inf))
801
(cond ((alike1 ll (m*t -1 '$minf))
803
(cond ((alike1 ul (m*t -1 '$minf))
805
(cond ((eq ul '$inf) nil)
807
(setq exp (subin (m- var) exp))
813
(setq exp (m- exp))))
814
;;Fix limits so that ll < ul.
815
(let ((d (complm ask-or-not)))
823
(defun complm (ask-or-not)
742
824
(let ((askflag (cond ((eq ask-or-not 'ask) t)
745
(COND ((ALIKE1 UL LL) 0.)
746
((EQ (SETQ A (cond (askflag ($asksign ($limit (m+t UL (m- LL)))))
747
(t ($sign ($limit (m+t UL (m- LL)))))))
753
(DEFUN INTSUBS (E A B)
827
(cond ((alike1 ul ll) 0.)
828
((eq (setq a (cond (askflag ($asksign ($limit (m+t ul (m- ll)))))
829
(t ($sign ($limit (m+t ul (m- ll)))))))
836
(defun intsubs (e a b)
754
837
(cond ((easy-subs e a b))
755
(t (setq current-assumptions
838
(t (setq *current-assumptions*
756
839
(make-defint-assumptions 'ask)) ;get forceful!
757
840
(let ((generate-atan2 ()) ($algebraic t)
758
841
(rpart ()) (ipart ()))
759
(desetq (rpart . ipart) (cond ((not (free e '$%i))
842
(desetq (rpart . ipart)
843
(cond ((not (free e '$%i))
762
846
(cond ((not (equal (sratsimp ipart) 0))
763
847
(let ((rans (cond ((limit-subs rpart a b))
812
896
(defun atan-pole1 (exp ll ul &aux ipart)
815
((matanp exp) ;neglect multiplicity and '$unknowns for now.
816
(desetq (exp . ipart) (trisplit exp))
818
((not (equal (sratsimp ipart) 0)) ())
819
(t (let ((pole (poles-in-interval (let (($algebraic t))
820
(sratsimp (cadr exp)))
822
(cond ((and pole (not (or (eq pole '$unknown)
824
(do ((l pole (cdr l)) (llist ()))
827
((eq (caar l) ll) t) ;Skip this one by definition.
828
(t (let ((low-lim ($limit (cadr exp) var (caar l) '$minus))
829
(up-lim ($limit (cadr exp) var (caar l) '$plus)))
830
(cond ((and (not (eq low-lim up-lim))
831
(real-infinityp low-lim)
832
(real-infinityp up-lim))
833
(let ((change (if (eq low-lim '$minf)
836
(setq llist (cons `((mequal simp) ,exp ,(m+ exp change))
838
(t (do ((l (cdr exp) (cdr l))
841
(setq llist (append llist (atan-pole1 (car l) ll ul)))))))
843
(DEFUN DIFAPPLY (N D S FN1)
844
(PROG (K M R $NOPRINCIPAL)
845
(COND ((eq ($asksign (m+ (DEG D) (m- S) (m- 2.))) '$neg)
847
(SETQ $NOPRINCIPAL T)
848
(COND ((OR (NOT (mexptp D))
849
(NOT (NUMBERP (SETQ R (CADDR D)))))
853
(EQUAL 1. (DEG (CADR D))))
855
(SETQ M (DEG (SETQ D (CADR D))))
856
(SETQ K (m// (m+ S 2.) M))
857
(COND ((eq (ask-integer (m// (m+ S 2.) M) '$any) '$yes)
859
(T (SETQ K (M+ 1 K))))
860
(COND ((eq ($sign (m+ r (m- k))) '$pos)
861
(RETURN (DIFFHK FN1 N D K (m+ R (m- K))))))))
863
(DEFUN DIFFHK (FN1 N D R M)
866
(SETQ D1 (FUNCALL FN1 N
869
(COND (D1 (RETURN (DIFAP1 D1 R '*Z* M 0.))))))
872
(COND ($NOPRINCIPAL (Diverg))
874
(PRINC "Principal Value")
879
(COND ((OR (MNUMP E) (CONSTANT E))
880
(SETQ BPTU (CONS E BPTU)))
881
(t (SETQ E (RMCONST1 E))
885
(SETQ E (CATCH 'PTIMES%E (PTIMES%E NN* ND*)))
888
(COND (*UPDN (SETQ BPTU (CONS E BPTU)))
889
(T (SETQ BPTD (CONS E BPTD))))))))))
891
(DEFUN PTIMES%E (TERM N)
892
(COND ((POLYINX TERM VAR NIL) TERM)
894
(EQ (CADR TERM) '$%E)
895
(POLYINX (CADDR TERM) VAR NIL)
896
(eq ($sign (m+ (DEG ($realpart (CADDR TERM))) -1.))
898
(eq ($sign (m+ (DEG (SETQ NN* ($imagpart (CADDR TERM))))
901
(COND ((EQ ($askSIGN (RATCOEF NN* VAR)) '$pos)
903
(T (SETQ *UPDN NIL)))
906
(SETQ NN* (POLFACTORS TERM))
908
(eq ($sign (m+ n (m- (deg (car nn*)))))
910
(PTIMES%E (CADR NN*) N)
912
(T (THROW 'PTIMES%E NIL))))
914
(DEFUN CSEMIDOWN (N D VAR)
899
((matanp exp) ;neglect multiplicity and '$unknowns for now.
900
(desetq (exp . ipart) (trisplit exp))
902
((not (equal (sratsimp ipart) 0)) ())
903
(t (let ((pole (poles-in-interval (let (($algebraic t))
904
(sratsimp (cadr exp)))
906
(cond ((and pole (not (or (eq pole '$unknown)
908
(do ((l pole (cdr l)) (llist ()))
911
((eq (caar l) ll) t) ;Skip this one by definition.
912
(t (let ((low-lim ($limit (cadr exp) var (caar l) '$minus))
913
(up-lim ($limit (cadr exp) var (caar l) '$plus)))
914
(cond ((and (not (eq low-lim up-lim))
915
(real-infinityp low-lim)
916
(real-infinityp up-lim))
917
(let ((change (if (eq low-lim '$minf)
920
(setq llist (cons `((mequal simp) ,exp ,(m+ exp change))
922
(t (do ((l (cdr exp) (cdr l))
925
(setq llist (append llist (atan-pole1 (car l) ll ul)))))))
927
(defun difapply (n d s fn1)
928
(prog (k m r $noprincipal)
929
(cond ((eq ($asksign (m+ (deg d) (m- s) (m- 2.))) '$neg)
931
(setq $noprincipal t)
932
(cond ((or (not (mexptp d))
933
(not (numberp (setq r (caddr d)))))
937
(equal 1. (deg (cadr d))))
939
(setq m (deg (setq d (cadr d))))
940
(setq k (m// (m+ s 2.) m))
941
(cond ((eq (ask-integer (m// (m+ s 2.) m) '$any) '$yes)
943
(t (setq k (m+ 1 k))))
944
(cond ((eq ($sign (m+ r (m- k))) '$pos)
945
(return (diffhk fn1 n d k (m+ r (m- k))))))))
947
(defun diffhk (fn1 n d r m)
950
(setq d1 (funcall fn1 n
953
(cond (d1 (return (difap1 d1 r '*z* m 0.))))))
956
(cond ($noprincipal (diverg))
958
(princ "Principal Value")
963
(cond ((or (mnump e) (constant e))
964
(setq bptu (cons e bptu)))
965
(t (setq e (rmconst1 e))
969
(setq e (catch 'ptimes%e (ptimes%e nn* nd*)))
972
(cond (*updn (setq bptu (cons e bptu)))
973
(t (setq bptd (cons e bptd))))))))))
975
(defun ptimes%e (term n)
976
(cond ((polyinx term var nil) term)
978
(eq (cadr term) '$%e)
979
(polyinx (caddr term) var nil)
980
(eq ($sign (m+ (deg ($realpart (caddr term))) -1.))
982
(eq ($sign (m+ (deg (setq nn* ($imagpart (caddr term))))
985
(cond ((eq ($asksign (ratcoef nn* var)) '$pos)
987
(t (setq *updn nil)))
990
(setq nn* (polfactors term))
992
(eq ($sign (m+ n (m- (deg (car nn*)))))
994
(ptimes%e (cadr nn*) n)
996
(t (throw 'ptimes%e nil))))
998
(defun csemidown (n d var)
915
999
(let ((pcprntd t)) ;Not sure what to do about PRINCIPAL values here.
916
(PRINCIP (RES N D #'LOWERHALF #'(lambda (X)
1000
(princip (res n d #'lowerhalf #'(lambda (x)
917
1001
(cond ((equal ($imagpart x) 0) t)
920
(DEFUN LOWERHALF (J) (EQ ($askSIGN ($imagpart J)) '$neg))
922
(DEFUN UPPERHALF (J) (EQ ($askSIGN ($imagpart J)) '$pos))
925
(DEFUN CSEMIUP (N D VAR)
1004
(defun lowerhalf (j)
1005
(eq ($asksign ($imagpart j)) '$neg))
1007
(defun upperhalf (j)
1008
(eq ($asksign ($imagpart j)) '$pos))
1011
(defun csemiup (n d var)
926
1012
(let ((pcprntd t)) ;I'm not sure what to do about PRINCIPAL values here.
927
(PRINCIP (RES N D #'UPPERHALF #'(lambda (X)
1013
(princip (res n d #'upperhalf #'(lambda (x)
928
1014
(cond ((equal ($imagpart x) 0) t)
933
(t (m*t '$%I ($RECTFORM (m+ (COND ((CAR N)
944
((POLYINX E VAR NIL) E)
948
(m+t (m^t '$%E (m*t '$%I (CADR E)))
949
(m- (m^t '$%E (m*t (m- '$%I) (CADR E)))))))
952
(m+t (m^t '$%E (m*t '$%I (CADR E)))
953
(m^t '$%E (m*t (m- '$%I) (CADR E))))))
955
(CONS (LIST (CAAR E)) (MAPCAR #'SCONVERT (CDR E)))))))
957
(DEFUN POLFACTORS (EXP)
1018
(cond ((null n) nil)
1019
(t (m*t '$%i ($rectform (m+ (cond ((car n)
1030
((polyinx e var nil) e)
1031
((eq (caar e) '%sin)
1034
(m+t (m^t '$%e (m*t '$%i (cadr e)))
1035
(m- (m^t '$%e (m*t (m- '$%i) (cadr e)))))))
1036
((eq (caar e) '%cos)
1037
(mul* '((rat) 1. 2.)
1038
(m+t (m^t '$%e (m*t '$%i (cadr e)))
1039
(m^t '$%e (m*t (m- '$%i) (cadr e))))))
1041
(cons (list (caar e)) (mapcar #'sconvert (cdr e)))))))
1043
(defun polfactors (exp)
958
1044
(let (poly rest)
959
(COND ((mplusp exp) nil)
960
(t (cond ((mtimesp EXP)
961
(SETQ EXP (REVERSE (CDR EXP))))
962
(T (SETQ EXP (LIST EXP))))
963
(mapc #'(lambda (term)
964
(cond ((polyinx term var nil)
966
(t (push term rest))))
968
(list (m*l poly) (m*l rest))))))
972
(COND ((ATOM E) (RETURN E))
973
((NOT (AMONG '$%E E)) (RETURN E))
976
(SETQ D ($imagpart (CADDR E)))
977
(RETURN (m* (m^t '$%E ($realpART (CADDR E)))
979
(m*t '$%I `((%SIN) ,D))))))
980
(T (RETURN (simplify (CONS (LIST (CAAR E))
981
(MAPCAR #'ESAP (CDR E)))))))))
1045
(cond ((mplusp exp) nil)
1046
(t (cond ((mtimesp exp)
1047
(setq exp (reverse (cdr exp))))
1048
(t (setq exp (list exp))))
1049
(mapc #'(lambda (term)
1050
(cond ((polyinx term var nil)
1052
(t (push term rest))))
1054
(list (m*l poly) (m*l rest))))))
1058
(cond ((atom e) (return e))
1059
((not (among '$%e e)) (return e))
1062
(setq d ($imagpart (caddr e)))
1063
(return (m* (m^t '$%e ($realpart (caddr e)))
1065
(m*t '$%i `((%sin) ,d))))))
1066
(t (return (simplify (cons (list (caar e))
1067
(mapcar #'esap (cdr e)))))))))
1069
(defun mtosc (grand)
987
1073
plf bptu bptd s upans downans)
988
(COND ((not (or (POLYINX D VAR NIL)
989
(AND (SETQ GRAND (%EINVOLVE D))
991
(POLYINX (SETQ D ($RATSIMP (M// D (m^t '$%E GRAND))))
994
(SETQ N (M// N (m^t '$%E GRAND)))))) nil)
995
((EQUAL (SETQ S (DEG D)) 0) NIL)
1074
(cond ((not (or (polyinx d var nil)
1075
(and (setq grand (%einvolve d))
1077
(polyinx (setq d ($ratsimp (m// d (m^t '$%e grand))))
1080
(setq n (m// n (m^t '$%e grand)))))) nil)
1081
((equal (setq s (deg d)) 0) nil)
996
1082
;;;Above tests for applicability of this method.
997
((and (or (SETQ PLF (POLFACTORS N)) t)
998
(SETQ N ($EXPAND (COND ((CAR PLF)
999
(m*t 'X* (SCONVERT (CADR PLF))))
1001
(COND ((mplusp N) (SETQ N (CDR N)))
1002
(T (SETQ N (LIST N))))
1083
((and (or (setq plf (polfactors n)) t)
1084
(setq n ($expand (cond ((car plf)
1085
(m*t 'x* (sconvert (cadr plf))))
1087
(cond ((mplusp n) (setq n (cdr n)))
1088
(t (setq n (list n))))
1003
1089
(do ((do-var n (cdr do-var)))
1004
1090
((null do-var) t)
1005
1091
(cond ((rib (car do-var) s))
1006
1092
(t (return nil))))
1007
1093
;;;Function RIB sets up the values of BPTU and BPTD
1009
(SETQ BPTU (SUBST (CAR PLF) 'X* BPTU))
1010
(SETQ BPTD (SUBST (CAR PLF) 'X* BPTD))
1011
t) ;CROCK, CROCK. This is TERRIBLE code.
1095
(setq bptu (subst (car plf) 'x* bptu))
1096
(setq bptd (subst (car plf) 'x* bptd))
1097
t) ;CROCK, CROCK. This is TERRIBLE code.
1013
1099
;;;If there is BPTU then CSEMIUP must succeed.
1014
1100
;;;Likewise for BPTD.
1021
1107
(sratsimp (m* '$%pi (m+ upans (m- downans))))))))
1024
(DEFUN EVENFN (E var)
1025
(let ((temp (m+ (m- E) (cond ((atom var)
1026
($substitute (m- var) var e))
1027
(t ($ratsubst (m- var) var E))))))
1028
(cond ((zerop1 temp)
1030
((zerop1 ($ratsimp temp))
1110
(defun evenfn (e var)
1111
(let ((temp (m+ (m- e)
1113
($substitute (m- var) var e))
1114
(t ($ratsubst (m- var) var e))))))
1115
(cond ((zerop1 temp)
1117
((zerop1 ($ratsimp temp))
1034
(DEFUN ODDFN (E VAR)
1035
(let ((temp (m+ e (cond ((atom var)
1036
($SUBSTITUTE (m- VAR) var E))
1037
(t ($ratsubst (m- var) var e))))))
1038
(cond ((zerop1 temp)
1040
((zerop1 ($ratsimp temp))
1044
(DEFUN ZTOINF (GRAND VAR)
1045
(PROG (N D SN* SD* VARLIST
1047
ANS R $SAVEFACTORS CHECKFACTORS temp test-var)
1048
(SETQ $SAVEFACTORS T SN* (SETQ SD* (LIST 1.)))
1049
(COND ((eq ($sign (m+ LOOPSTOP* -1))
1052
((SETQ temp (OR (SCAXN GRAND)
1055
((INVOLVE GRAND '(%SIN %COS %TAN))
1056
(SETQ GRAND (SCONVERT GRAND))
1059
(COND ((POLYINX GRAND VAR NIL)
1061
((AND (RATP GRAND VAR)
1063
(ANDMAPCAR #'SNUMDEN (CDR GRAND)))
1121
(defun oddfn (e var)
1122
(let ((temp (m+ e (cond ((atom var)
1123
($substitute (m- var) var e))
1124
(t ($ratsubst (m- var) var e))))))
1125
(cond ((zerop1 temp)
1127
((zerop1 ($ratsimp temp))
1131
(defun ztoinf (grand var)
1132
(prog (n d sn* sd* varlist
1134
ans r $savefactors checkfactors temp test-var)
1135
(setq $savefactors t sn* (setq sd* (list 1.)))
1136
(cond ((eq ($sign (m+ loopstop* -1))
1139
((setq temp (or (scaxn grand)
1142
((involve grand '(%sin %cos %tan))
1143
(setq grand (sconvert grand))
1146
(cond ((polyinx grand var nil)
1148
((and (ratp grand var)
1150
(andmapcar #'snumden (cdr grand)))
1070
1157
;;;New section.
1071
(SETQ N (RMCONST1 NN*))
1072
(SETQ D (RMCONST1 DN*))
1077
(COND ((POLYINX D VAR NIL)
1080
(COND ((AND (SETQ R (FINDP N))
1081
(eq (ask-integer R '$integer) '$yes)
1082
(SETQ test-var (BXM D S))
1083
(SETQ ans (APPLY 'FAN (CONS (m+ 1. R) test-var))))
1084
(RETURN (m* (M// NC DC) ($RATSIMP ans))))
1085
((and (RATP GRAND VAR)
1086
(SETQ ANS (ZMTORAT N (COND ((mtimesp d) d)
1089
(RETURN (m* (M// NC DC) ANS)))
1090
((AND (EVENFN D VAR)
1091
(SETQ NN* (P*LOGNXP N S)))
1092
(SETQ ANS (LOG*RAT (CAR NN*) D (CADR NN*)))
1093
(RETURN (m* (M// NC DC) ANS)))
1094
((INVOLVE GRAND '(%LOG))
1095
(COND ((SETQ ANS (LOGQUAD0 GRAND))
1096
(RETURN (m* (M// NC DC) ANS)))
1099
(COND ((SETQ temp (BATAPP GRAND))
1103
(COND ((let ((MTOINF* nil))
1104
(SETQ temp (GGR GRAND T)))
1107
(COND ((let ((*NODIVERG t))
1108
(SETQ ANS (CATCH 'Divergent
1109
(ANDMAPCAR #'(LAMBDA (G)
1112
(COND ((EQ ANS 'Divergent) NIL)
1113
(T (RETURN (sRATSIMP (m+l ANS)))))))))
1158
(setq n (rmconst1 nn*))
1159
(setq d (rmconst1 dn*))
1164
(cond ((polyinx d var nil)
1167
(cond ((and (setq r (findp n))
1168
(eq (ask-integer r '$integer) '$yes)
1169
(setq test-var (bxm d s))
1170
(setq ans (apply 'fan (cons (m+ 1. r) test-var))))
1171
(return (m* (m// nc dc) ($ratsimp ans))))
1172
((and (ratp grand var)
1173
(setq ans (zmtorat n (cond ((mtimesp d) d)
1176
(return (m* (m// nc dc) ans)))
1177
((and (evenfn d var)
1178
(setq nn* (p*lognxp n s)))
1179
(setq ans (log*rat (car nn*) d (cadr nn*)))
1180
(return (m* (m// nc dc) ans)))
1181
((involve grand '(%log))
1182
(cond ((setq ans (logquad0 grand))
1183
(return (m* (m// nc dc) ans)))
1186
(cond ((setq temp (batapp grand))
1190
(cond ((let ((*mtoinf* nil))
1191
(setq temp (ggr grand t)))
1194
(cond ((let ((*nodiverg t))
1195
(setq ans (catch 'divergent
1196
(andmapcar #'(lambda (g)
1199
(cond ((eq ans 'divergent) nil)
1200
(t (return (sratsimp (m+l ans)))))))))
1115
(COND ((AND (EVENFN GRAND VAR)
1116
(SETQ LOOPSTOP* (M+ 1 LOOPSTOP*))
1117
(SETQ ANS (method-by-limits grand var '$minf '$inf)))
1118
(return (m*t '((RAT) 1. 2.) ANS)))
1202
(cond ((and (evenfn grand var)
1203
(setq loopstop* (m+ 1 loopstop*))
1204
(setq ans (method-by-limits grand var '$minf '$inf)))
1205
(return (m*t '((rat) 1. 2.) ans)))
1121
(DEFUN ZTORAT (N D S)
1122
(COND ((AND (NULL *DFLAG)
1123
(SETQ S (DIFAPPLY N D NN* #'ZTORAT)))
1125
((SETQ N (let ((PLOGABS ()))
1126
(KEYHOLE (m* `((%PLOG) ,(m- VAR)) N) D VAR)))
1128
(T (MERROR "Keyhole failed"))))
1132
(DEFUN LOGQUAD0 (EXP)
1208
(defun ztorat (n d s)
1209
(cond ((and (null *dflag)
1210
(setq s (difapply n d nn* #'ztorat)))
1212
((setq n (let ((plogabs ()))
1213
(keyhole (m* `((%plog) ,(m- var)) n) d var)))
1215
(t (merror "Keyhole failed"))))
1219
(defun logquad0 (exp)
1133
1220
(let ((a ()) (b ()) (c ()))
1134
(COND ((SETQ EXP (LOGQUAD EXP))
1135
(SETQ A (CAR EXP) B (CADR EXP) C (CADDR EXP))
1136
($asksign b) ;let the data base know about the sign of B.
1137
(COND ((EQ ($askSIGN C) '$pos)
1138
(SETQ C (m^ (M// C A) '((RAT) 1. 2.)))
1221
(cond ((setq exp (logquad exp))
1222
(setq a (car exp) b (cadr exp) c (caddr exp))
1223
($asksign b) ;let the data base know about the sign of B.
1224
(cond ((eq ($asksign c) '$pos)
1225
(setq c (m^ (m// c a) '((rat) 1. 2.)))
1139
1226
(setq b (simplify
1140
`((%ACOS) ,(add* 'epsilon (M// B (mul* 2. A C))))))
1141
(setq a (M// (m* B `((%LOG) ,C))
1142
(mul* A (Simplify `((%SIN) ,B)) C)))
1227
`((%acos) ,(add* 'epsilon (m// b (mul* 2. a c))))))
1228
(setq a (m// (m* b `((%log) ,c))
1229
(mul* a (simplify `((%sin) ,b)) c)))
1143
1230
(get-limit a 'epsilon 0 '$plus))))
1146
(DEFUN LOGQUAD (EXP)
1147
(let ((VARLIST (list var)))
1149
(SETQ EXP (CDR (RATREP* EXP)))
1150
(COND ((AND (ALIKE1 (PDIS (CAR EXP))
1152
(NOT (ATOM (CDR EXP)))
1153
(EQUAL (CADR (CDR EXP)) 2.)
1233
(defun logquad (exp)
1234
(let ((varlist (list var)))
1236
(setq exp (cdr (ratrep* exp)))
1237
(cond ((and (alike1 (pdis (car exp))
1239
(not (atom (cdr exp)))
1240
(equal (cadr (cdr exp)) 2.)
1154
1241
(not (equal (pterm (cddr exp) 0.) 0.)))
1155
(SETQ EXP (MAPCAR 'PDIS (CDR (ODDELM (CDR EXP)))))))))
1157
(DEFUN MTOINF (GRAND VAR)
1158
(PROG (ANS SD* SN* P* PE* N D S NC DC $SAVEFACTORS CHECKFACTORS temp)
1159
(SETQ $SAVEFACTORS T)
1160
(SETQ SN* (SETQ SD* (LIST 1.)))
1161
(COND ((eq ($sign (m+ LOOPSTOP* -1)) '$pos)
1163
((INVOLVE GRAND '(%SIN %COS))
1164
(COND ((AND (EVENFN GRAND VAR)
1165
(or (SETQ temp (SCAXN GRAND))
1166
(setq temp (ssp grand))))
1167
(RETURN (m*t 2. temp)))
1168
((SETQ temp (MTOSC GRAND))
1171
((among '$%i (%EINVOLVE GRAND))
1172
(COND ((SETQ temp (MTOSC GRAND))
1175
(COND ((POLYINX GRAND VAR NIL)
1177
((AND (RATP GRAND VAR)
1179
(ANDMAPCAR #'SNUMDEN (CDR GRAND)))
1180
(SETQ NN* (M*L SN*) SN* NIL)
1181
(SETQ DN* (M*L SD*) SD* NIL))
1183
(SETQ N (RMCONST1 NN*))
1184
(SETQ D (RMCONST1 DN*))
1189
(COND ((POLYINX D VAR NIL)
1191
(COND ((AND (NOT (%einvolve grand))
1192
(NOTINVOLVE EXP '(%SINH %COSH %TANH))
1194
(eq (ask-integer P* '$integer) '$yes)
1195
(SETQ PE* (BXM D S)))
1196
(COND ((AND (eq (ask-integer (CADDR PE*) '$even) '$yes)
1197
(eq (ask-integer P* '$even) '$yes))
1198
(COND ((SETQ ANS (APPLY 'FAN (CONS (m+ 1. P*) PE*)))
1199
(SETQ ANS (m*t 2. ANS))
1200
(RETURN (m* (M// NC DC) ANS)))))
1201
((EQUAL (CAR PE*) 1.)
1202
(COND ((AND (SETQ ANS (APPLY 'FAN (CONS (m+ 1. P*) PE*)))
1203
(SETQ NN* (FAN (m+ 1. P*)
1208
(SETQ ANS (m+ ANS (m*t (m^ -1. P*) NN*)))
1209
(RETURN (m* (M// NC DC) ANS))))))))
1210
(COND ((RATP GRAND VAR)
1211
(SETQ ANS (m*t '$%PI (ZMTORAT N (COND ((mtimesp d) d)
1215
(RETURN (m* (M// NC DC) ANS)))
1216
((AND (OR (%einvolve grand)
1217
(INVOLVE GRAND '(%SINH %COSH %TANH)))
1218
(P*PIN%EX N) ;setq's P* and PE*...Barf again.
1219
(SETQ ANS (CATCH 'PIN%EX (PIN%EX D))))
1221
(RETURN (DINTEXP GRAND VAR)))
1222
((AND (ZEROP1 (get-LIMIT GRAND VAR '$INF))
1223
(ZEROP1 (get-LIMIT GRAND VAR '$MINF))
1224
(SETQ ANS (RECTZTO%PI2 (M*L P*) (M*L PE*) D)))
1225
(RETURN (m* (M// NC DC) ANS)))
1228
(COND ((SETQ ANS (GGRM GRAND))
1230
((AND (EVENFN GRAND VAR)
1231
(SETQ LOOPSTOP* (M+ 1 LOOPSTOP*))
1232
(SETQ ANS (method-by-limits grand var 0 '$inf)))
1233
(RETURN (m*t 2. ANS)))
1236
(DEFUN LINPOWER0 (EXP VAR)
1237
(COND ((AND (SETQ EXP (LINPOWER EXP VAR))
1238
(eq (ask-integer (CADDR EXP) '$even)
1240
(RATGREATERP 0. (CAR EXP)))
1242
(setq exp (mapcar 'pdis (cdr (oddelm (cdr exp)))))))))
1244
(defun mtoinf (grand var)
1245
(prog (ans sd* sn* p* pe* n d s nc dc $savefactors checkfactors temp)
1246
(setq $savefactors t)
1247
(setq sn* (setq sd* (list 1.)))
1248
(cond ((eq ($sign (m+ loopstop* -1)) '$pos)
1250
((involve grand '(%sin %cos))
1251
(cond ((and (evenfn grand var)
1252
(or (setq temp (scaxn grand))
1253
(setq temp (ssp grand))))
1254
(return (m*t 2. temp)))
1255
((setq temp (mtosc grand))
1258
((among '$%i (%einvolve grand))
1259
(cond ((setq temp (mtosc grand))
1262
(cond ((polyinx grand var nil)
1264
((and (ratp grand var)
1266
(andmapcar #'snumden (cdr grand)))
1267
(setq nn* (m*l sn*) sn* nil)
1268
(setq dn* (m*l sd*) sd* nil))
1270
(setq n (rmconst1 nn*))
1271
(setq d (rmconst1 dn*))
1276
(cond ((polyinx d var nil)
1278
(cond ((and (not (%einvolve grand))
1279
(notinvolve exp '(%sinh %cosh %tanh))
1281
(eq (ask-integer p* '$integer) '$yes)
1282
(setq pe* (bxm d s)))
1283
(cond ((and (eq (ask-integer (caddr pe*) '$even) '$yes)
1284
(eq (ask-integer p* '$even) '$yes))
1285
(cond ((setq ans (apply 'fan (cons (m+ 1. p*) pe*)))
1286
(setq ans (m*t 2. ans))
1287
(return (m* (m// nc dc) ans)))))
1288
((equal (car pe*) 1.)
1289
(cond ((and (setq ans (apply 'fan (cons (m+ 1. p*) pe*)))
1290
(setq nn* (fan (m+ 1. p*)
1295
(setq ans (m+ ans (m*t (m^ -1. p*) nn*)))
1296
(return (m* (m// nc dc) ans))))))))
1297
(cond ((ratp grand var)
1298
(setq ans (m*t '$%pi (zmtorat n (cond ((mtimesp d) d)
1302
(return (m* (m// nc dc) ans)))
1303
((and (or (%einvolve grand)
1304
(involve grand '(%sinh %cosh %tanh)))
1305
(p*pin%ex n) ;setq's P* and PE*...Barf again.
1306
(setq ans (catch 'pin%ex (pin%ex d))))
1307
;; We have an integral of the form p(x)*F(exp(x)), where
1308
;; p(x) is a polynomial.
1311
(return (dintexp grand var)))
1312
((not (and (zerop1 (get-limit grand var '$inf))
1313
(zerop1 (get-limit grand var '$minf))))
1314
;; These limits must exist for the integral to converge.
1316
((setq ans (rectzto%pi2 (m*l p*) (m*l pe*) d))
1317
;; This only handles the case when the F(z) is a
1318
;; rational function.
1319
(return (m* (m// nc dc) ans)))
1320
((setq ans (log-transform (m*l p*) (m*l pe*) d))
1321
;; If we get here, F(z) is not a rational function.
1322
;; We transform it using the substitution x=log(y)
1323
;; which gives us an integral of the form
1324
;; p(log(y))*F(y)/y, which maxima should be able to
1326
(return (m* (m// nc dc) ans)))
1328
;; Give up. We don't know how to handle this.
1331
(cond ((setq ans (ggrm grand))
1333
((and (evenfn grand var)
1334
(setq loopstop* (m+ 1 loopstop*))
1335
(setq ans (method-by-limits grand var 0 '$inf)))
1336
(return (m*t 2. ans)))
1339
(defun linpower0 (exp var)
1340
(cond ((and (setq exp (linpower exp var))
1341
(eq (ask-integer (caddr exp) '$even)
1343
(ratgreaterp 0. (car exp)))
1243
1346
;;; given (b*x+a)^n+c returns (a b n c)
1244
(DEFUN LINPOWER (EXP VAR)
1245
(let (LINPART DEG LC C VARLIST)
1246
(COND ((NOT (POLYP EXP)) NIL)
1247
(t (let ((varlist (list var)))
1249
(SETQ LINPART (CADR (RATREP* EXP)))
1250
(COND ((ATOM LINPART)
1252
(t (SETQ DEG (CADR LINPART))
1347
(defun linpower (exp var)
1348
(let (linpart deg lc c varlist)
1349
(cond ((not (polyp exp)) nil)
1350
(t (let ((varlist (list var)))
1352
(setq linpart (cadr (ratrep* exp)))
1353
(cond ((atom linpart)
1355
(t (setq deg (cadr linpart))
1253
1356
;;;get high degree of poly
1254
(SETQ LINPART ($DIFF EXP VAR (m+ deg -1)))
1357
(setq linpart ($diff exp var (m+ deg -1)))
1255
1358
;;;diff down to linear.
1256
(SETQ LC (SDIFF LINPART VAR))
1359
(setq lc (sdiff linpart var))
1257
1360
;;;all the way to constant.
1258
(SETQ LINPART ($RATSIMP (M// LINPART lc)))
1259
(SETQ LC ($RATSIMP (M// LC `((MFACTORIAL) ,DEG))))
1361
(setq linpart ($ratsimp (m// linpart lc)))
1362
(setq lc ($ratsimp (m// lc `((mfactorial) ,deg))))
1260
1363
;;;get rid of factorial from differentiation.
1261
(SETQ C ($RATSIMP (m+ EXP (m* (m- LC)
1262
(m^ LINPART DEG)))))))
1364
(setq c ($ratsimp (m+ exp (m* (m- lc)
1365
(m^ linpart deg)))))))
1263
1366
;;;Sees if can be expressed as (a*x+b)^n + part freeof x.
1264
(COND ((NOT (AMONG VAR c))
1265
`(,LC ,LINPART ,DEG ,C))
1268
(DEFUN MTORAT (N D S)
1270
(COND ((AND (NULL *DFLAG)
1271
(SETQ S (DIFAPPLY N D s #'MTORAT)))
1273
(T (CSEMIUP N D VAR)))))
1275
(DEFUN ZMTORAT (N D S FN1)
1277
(COND ((eq ($sign (m+ S (M+ 1 (SETQ NN* (DEG N)))))
1280
((eq ($sign (m+ s -4))
1283
(SETQ D ($FACTOR D))
1284
(SETQ C (RMCONST1 D))
1290
(SETQ N (PARTNUM N D))
1292
(SETQ N ($XTHRU (M+L
1293
(MAPCAR #'(LAMBDA (A B)
1294
(M// (FUNCALL FN1 (CAR A) B (DEG B))
1298
(RETURN (COND (C (M// N C))
1302
(SETQ N (FUNCALL FN1 N D S))
1303
(RETURN (sRATSIMP (COND (C (M// N C))
1306
(DEFUN PFRNUM (F G N N2 VAR)
1307
(let ((VARLIST (list var)) GENVAR)
1308
(SETQ F (POLYFORM F)
1312
(SETQ VAR (CAADR (RATREP* VAR)))
1313
(SETQ F (RESPROG0 F G N N2))
1314
(LIST (LIST (PDIS (CADR F)) (PDIS (CDDR F)))
1315
(LIST (PDIS (CAAR F)) (PDIS (CDAR F))))))
1320
(SETQ F (RATREP* E))
1321
(AND (EQUAL (CDDR F) 1)
1323
(AND (EQUAL (LENGTH (SETQ D (CDDR F))) 3)
1326
(RETURN (LIST (CAR D)
1328
(PTIMES (CADR F) (CADDR D)))))
1329
(MERROR "Bug from PFRNUM in RESIDU")))
1331
(DEFUN PARTNUM (N DL)
1332
(let ((n2 1) ANS NL)
1333
(do ((dl dl (cdr dl)))
1335
(nconc ans (ncons (list n n2))))
1336
(SETQ NL (PFRNUM (CAR DL) (M*L (CDR DL)) N N2 VAR))
1337
(SETQ ANS (NCONC ANS (NCONS (CAR NL))))
1338
(SETQ N2 (CADADR NL) N (CAADR NL) NL NIL))))
1341
(PROG (POLY EXPO MTOINF* MB VARLIST GENVAR L C GVAR)
1342
(SETQ VARLIST (LIST VAR))
1344
(COND ((AND (SETQ EXPO (%EINVOLVE E))
1345
(POLYP (SETQ POLY ($RATSIMP (M// E (m^t '$%E EXPO)))))
1346
(SETQ L (CATCH 'GGRM (GGR (m^t '$%E EXPO) NIL))))
1348
(SETQ MB (m- (SUBIN 0. (CADR L))))
1349
(SETQ POLY (m+ (SUBIN (m+t MB VAR) POLY)
1350
(SUBIN (m+t MB (m*t -1. VAR)) POLY))))
1352
(SETQ EXPO (CADDR L)
1357
(SETQ POLY (CDR (RATREP* POLY)))
1358
(SETQ MB (m^ (PDIS (CDR POLY)) -1.)
1360
(SETQ GVAR (CAADR (RATREP* VAR)))
1361
(COND ((OR (ATOM POLY)
1362
(POINTERGP GVAR (CAR POLY)))
1363
(SETQ POLY (LIST 0. POLY)))
1364
(T (SETQ POLY (CDR POLY))))
1365
(return (do ((poly poly (cddr poly)))
1367
(mul* (m^t '$%E C) (m^t EXPO -1.) MB (M+L E)))
1368
(SETQ E (CONS (GGRM1 (CAR POLY) (PDIS (CADR POLY)) L EXPO)
1371
(DEFUN GGRM1 (D K A B)
1372
(SETQ B (M// (m+t 1. D) B))
1373
(m* K `((%GAMMA) ,B) (m^ A (m- B))))
1376
;;;If rd* is t the m^ts must just be free of var.
1377
;;;If rd* is () the m^ts must be mnump's.
1381
(DEFUN KEYHOLE (N D VAR)
1382
(let ((SEMIRAT* ()))
1383
(SETQ N (RES N D #'(LAMBDA (J)
1384
(OR (NOT (equal ($imagpart j) 0))
1385
(EQ ($askSIGN J) '$neg)))
1387
(COND ((EQ ($askSIGN J) '$pos)
1391
($rectform ($multthru (m+ (COND ((CAR N) (CAR N))
1393
(COND ((CADR N) (CADR N))
1398
(COND ((ATOM E) (RETURN NIL)))
1399
(SETQ E (PARTITION E VAR 1))
1402
(COND ((SETQ R (SINRX E)) (RETURN (LIST M R 1)))
1404
(eq (ask-integer (SETQ K (CADDR E)) '$integer) '$yes)
1405
(SETQ R (SINRX (CADR E))))
1406
(RETURN (LIST M R K))))))
1409
(COND ((and (consp e) (EQ (CAAR E) '%SIN))
1410
(COND ((EQ (CADR E) VAR) 1.)
1411
((AND (SETQ E (PARTITION (CADR E) VAR 1)) (EQ (CDR E) VAR))
1414
(DECLARE-TOP(SPECIAL N))
1419
(setq exp ($substitute (m^t `((%sin) ,var) 2.)
1420
(m+t 1. (m- (m^t `((%cos) ,var) 2.)))
1424
(COND ((AND (SETQ N (FINDP DN*))
1425
(eq (ask-integer N '$integer) '$yes))
1426
(COND ((SETQ C (SKR U))
1427
(RETURN (SCMP C N)))
1429
(SETQ C (ANDMAPCAR #'SKR (CDR U))))
1430
(RETURN (m+l (MAPCAR #'(LAMBDA (J) (SCMP J N))
1433
(DECLARE-TOP(UNSPECIAL N))
1436
(m* (CAR C) (m^ (CADR C) (m+ N -1)) `((%SIGNUM) ,(CADR C))
1437
(SINSP (CADDR C) N)))
1440
(m* HALF%PI ($MAKEGAMMA `((%BINOMIAL) ,(m+t (m+ N -1) '((RAT) -1. 2.))
1447
(T (BYGAMMA (m+ N -1) 0.))))
1452
(COND ((eq ($sign (m+ L (m- (m+ K -1))))
1367
(cond ((not (among var c))
1368
`(,lc ,linpart ,deg ,c))
1371
(defun mtorat (n d s)
1373
(cond ((and (null *dflag)
1374
(setq s (difapply n d s #'mtorat)))
1376
(t (csemiup n d var)))))
1378
(defun zmtorat (n d s fn1)
1380
(cond ((eq ($sign (m+ s (m+ 1 (setq nn* (deg n)))))
1383
((eq ($sign (m+ s -4))
1386
(setq d ($factor d))
1387
(setq c (rmconst1 d))
1393
(setq n (partnum n d))
1395
(setq n ($xthru (m+l
1396
(mapcar #'(lambda (a b)
1397
(m// (funcall fn1 (car a) b (deg b))
1401
(return (cond (c (m// n c))
1405
(setq n (funcall fn1 n d s))
1406
(return (sratsimp (cond (c (m// n c))
1409
(defun pfrnum (f g n n2 var)
1410
(let ((varlist (list var)) genvar)
1411
(setq f (polyform f)
1415
(setq var (caadr (ratrep* var)))
1416
(setq f (resprog0 f g n n2))
1417
(list (list (pdis (cadr f)) (pdis (cddr f)))
1418
(list (pdis (caar f)) (pdis (cdar f))))))
1423
(setq f (ratrep* e))
1424
(and (equal (cddr f) 1)
1426
(and (equal (length (setq d (cddr f))) 3)
1429
(return (list (car d)
1431
(ptimes (cadr f) (caddr d)))))
1432
(merror "Bug from `pfrnum' in `residu'")))
1434
(defun partnum (n dl)
1435
(let ((n2 1) ans nl)
1436
(do ((dl dl (cdr dl)))
1438
(nconc ans (ncons (list n n2))))
1439
(setq nl (pfrnum (car dl) (m*l (cdr dl)) n n2 var))
1440
(setq ans (nconc ans (ncons (car nl))))
1441
(setq n2 (cadadr nl) n (caadr nl) nl nil))))
1444
(prog (poly expo *mtoinf* mb varlist genvar l c gvar)
1445
(setq varlist (list var))
1447
(cond ((and (setq expo (%einvolve e))
1448
(polyp (setq poly ($ratsimp (m// e (m^t '$%e expo)))))
1449
(setq l (catch 'ggrm (ggr (m^t '$%e expo) nil))))
1451
(setq mb (m- (subin 0. (cadr l))))
1452
(setq poly (m+ (subin (m+t mb var) poly)
1453
(subin (m+t mb (m*t -1. var)) poly))))
1455
(setq expo (caddr l)
1460
(setq poly (cdr (ratrep* poly)))
1461
(setq mb (m^ (pdis (cdr poly)) -1.)
1463
(setq gvar (caadr (ratrep* var)))
1464
(cond ((or (atom poly)
1465
(pointergp gvar (car poly)))
1466
(setq poly (list 0. poly)))
1467
(t (setq poly (cdr poly))))
1468
(return (do ((poly poly (cddr poly)))
1470
(mul* (m^t '$%e c) (m^t expo -1.) mb (m+l e)))
1471
(setq e (cons (ggrm1 (car poly) (pdis (cadr poly)) l expo)
1474
(defun ggrm1 (d k a b)
1475
(setq b (m// (m+t 1. d) b))
1476
(m* k `((%gamma) ,b) (m^ a (m- b))))
1479
;;If rd* is t the m^ts must just be free of var.
1480
;;If rd* is () the m^ts must be mnump's.
1484
(defun keyhole (n d var)
1485
(let ((semirat* ()))
1486
(setq n (res n d #'(lambda (j)
1487
(or (not (equal ($imagpart j) 0))
1488
(eq ($asksign j) '$neg)))
1490
(cond ((eq ($asksign j) '$pos)
1494
($rectform ($multthru (m+ (cond ((car n) (car n))
1496
(cond ((cadr n) (cadr n))
1499
;; Look at an expression e of the form sin(r*x)^k, where k is an
1500
;; integer. Return the list (1 r k). (Not sure if the first element
1501
;; of the list is always 1 because I'm not sure what partition is
1502
;; trying to do here.)
1505
(cond ((atom e) (return nil)))
1506
(setq e (partition e var 1))
1509
(cond ((setq r (sinrx e))
1510
(return (list m r 1)))
1512
(eq (ask-integer (setq k (caddr e)) '$integer) '$yes)
1513
(setq r (sinrx (cadr e))))
1514
(return (list m r k))))))
1516
;; Look at an expression e of the form sin(r*x) and return r.
1518
(cond ((and (consp e) (eq (caar e) '%sin))
1519
(cond ((eq (cadr e) var)
1521
((and (setq e (partition (cadr e) var 1))
1527
;; integrate(a*sc(r*x)^k/x^n,x,0,inf).
1530
;; I don't think this needs to be special.
1532
(declare (special n))
1533
;; Replace (1-cos(x)^2) with sin(x)^2.
1534
(setq exp ($substitute (m^t `((%sin) ,var) 2.)
1535
(m+t 1. (m- (m^t `((%cos) ,var) 2.)))
1539
(cond ((and (setq n (findp dn*))
1540
(eq (ask-integer n '$integer) '$yes))
1541
;; n is the power of the denominator.
1542
(cond ((setq c (skr u))
1544
(return (scmp c n)))
1546
(setq c (andmapcar #'skr (cdr u))))
1547
;; Do this for a sum of such terms.
1548
(return (m+l (mapcar #'(lambda (j) (scmp j n))
1551
;; We have an integral of the form sin(r*x)^k/x^n. C is the list (1 r k).
1553
;; The substitution y=r*x converts this integral to
1555
;; r^(n-1)*integral(sin(y)^k/y^n,y,0,inf)
1557
;; (If r is negative, we need to negate the result.)
1559
;; The recursion Wang gives on p. 87 has a typo. The second term
1560
;; should be subtracted from the first. This formula is given in G&R,
1561
;; 3.82, formula 12.
1563
;; integrate(sin(x)^r/x^s,x) =
1564
;; r*(r-1)/(s-1)/(s-2)*integrate(sin(x)^(r-2)/x^(s-2),x)
1565
;; - r^2/(s-1)/(s-2)*integrate(sin(x)^r/x^(s-2),x)
1567
;; (Limits are assumed to be 0 to inf.)
1569
;; This recursion ends up with integrals with s = 1 or 2 and
1571
;; integrate(sin(x)^p/x,x,0,inf) = integrate(sin(x)^(p-1),x,0,%pi/2)
1573
;; with p > 0 and odd. This latter integral is known to maxima, and
1574
;; it's value is beta(p/2,1/2)/2.
1576
;; integrate(sin(x)^2/x^2,x,0,inf) = %pi/2*binomial(q-3/2,q-1)
1581
;; Compute sign(r)*r^(n-1)*integrate(sin(y)^k/y^n,y,0,inf)
1583
(m^ (cadr c) (m+ n -1))
1584
`((%signum) ,(cadr c))
1585
(sinsp (caddr c) n)))
1587
;; integrate(sin(x)^n/x^2,x,0,inf) = pi/2*binomial(n-3/2,n-1).
1588
;; Express in terms of Gamma functions, though.
1590
(m* half%pi ($makegamma `((%binomial) ,(m+t (m+ n -1) '((rat) -1. 2.))
1594
;; integrate(sin(x)^n/x,x,0,inf) = beta((n+1)/2,1/2)/2, for n odd and
1599
(t (bygamma (m+ n -1) 0.))))
1601
;; This implements the recursion for computing
1602
;; integrate(sin(y)^l/y^k,y,0,inf). (Note the change in notation from
1607
(cond ((eq ($sign (m+ l (m- (m+ k -1))))
1455
((NOT (EVEN1 (m+ L K)))
1461
((eq ($sign (m+ K -2.))
1609
;; Integral diverges if l-(k-1) < 0.
1611
((not (even1 (m+ l k)))
1612
;; If l + k is not even, return NIL. (Is this the right
1616
;; We have integrate(sin(y)^l/y^2). Use sevn to evaluate.
1619
;; We have integrate(sin(y)^l/y,y)
1621
((eq ($sign (m+ k -2.))
1463
(SETQ I (m* (m+ K -1) (SETQ J (m+ K -2.))))
1464
(m+ (m* L (m+ L -1) (m^t I -1.) (SINSP (m+ L -2.) J))
1465
(m* (m- (m^ L 2)) (m^t I -1.)
1472
(LIST (CAR A) (REMAINDER (CADR A) (CADDR A)) (CADDR A)))
1473
((AND (ATOM A) (ABLESS1 A)) A)
1476
(ABLESS1 (CADDR A)))
1480
(COND ((POLYINX E VAR NIL) 0.)
1486
(m+l (MAPCAR #'THRAD E)))))
1623
(setq i (m* (m+ k -1)
1624
(setq j (m+ k -2.))))
1625
;; j = k-2, i = (k-1)*(k-2)
1628
;; The main recursion:
1631
;; = l*(l-1)/(k-1)/(k-2)*i(sin(y)^(l-2)/y^k)
1632
;; - l^2/(k-1)/(k-1)*i(sin(y)^l/y^(k-2))
1635
(sinsp (m+ l -2.) j))
1644
(list (car a) (remainder (cadr a) (caddr a)) (caddr a)))
1645
((and (atom a) (abless1 a)) a)
1648
(abless1 (caddr a)))
1652
(cond ((polyinx e var nil) 0.)
1658
(m+l (mapcar #'thrad e)))))
1489
1661
;;; THE FOLLOWING FUNCTION IS FOR TRIG FUNCTIONS OF THE FOLLOWING TYPE:
1490
1662
;;; LOWER LIMIT=0 B A MULTIPLE OF %PI SCA FUNCTION OF SIN (X) COS (X)
1493
(DEFUN PERIOD (P E VAR)
1494
(and (ALIKE1 (no-err-sub var E) (setq e (NO-ERR-SUB (m+ P VAR) E)))
1665
(defun period (p e var)
1666
(and (alike1 (no-err-sub var e) (setq e (no-err-sub (m+ p var) e)))
1495
1667
;; means there was no error
1496
1668
(not (eq e t))))
1499
1671
(let ((var '$%i)
1500
1672
(r (subin 0. a))
1502
(SETQ C (SUBIN 1. (m+ A (m*t -1. R))))
1503
(SETQ A (IGPRT (m* '((RAT) 1. 2.) C)))
1504
(CONS A (m+ R (m*t (M+ C (m* -2. A)) '$%PI)))))
1674
(setq c (subin 1. (m+ a (m*t -1. r))))
1675
(setq a (igprt (m* '((rat) 1. 2.) c)))
1676
(cons a (m+ r (m*t (m+ c (m* -2. a)) '$%pi)))))
1507
(M+ R (m* -1. (FPART R))))
1679
(m+ r (m* -1. (fpart r))))
1510
1682
;;;Try making exp(%i*var) --> yy, if result is rational then do integral
1511
1683
;;;around unit circle. Make corrections for limits of integration if possible.
1513
(let* ((exp-form (sconvert sc)) ;Exponentialize
1514
(rat-form (MAXIMA-SUBSTITUTE 'yy (m^t '$%e (m*t '$%i var))
1515
exp-form))) ;Try to make Rational fun.
1516
(COND ((AND (RATP rat-form 'YY)
1517
(NOT (AMONG VAR rat-form)))
1518
(COND ((ALIKE1 B %PI2)
1519
(let ((ans (ZTO%PI2 rat-form 'YY)))
1523
(EVENFN exp-form VAR))
1524
(let ((ans (ZTO%PI2 rat-form 'YY)))
1525
(cond (ans (m*t '((RAT) 1. 2.) ans))
1527
((AND (ALIKE1 B HALF%PI)
1528
(EVENFN exp-form VAR)
1530
(NO-ERR-SUB (m+t '$%PI (m*t -1. VAR))
1532
(let ((ans (ZTO%PI2 rat-form 'yy)))
1533
(cond (ans (m*t '((RAT) 1. 4.) ans))
1685
(let* ((exp-form (sconvert sc)) ;Exponentialize
1686
(rat-form (maxima-substitute 'yy (m^t '$%e (m*t '$%i var))
1687
exp-form))) ;Try to make Rational fun.
1688
(cond ((and (ratp rat-form 'yy)
1689
(not (among var rat-form)))
1690
(cond ((alike1 b %pi2)
1691
(let ((ans (zto%pi2 rat-form 'yy)))
1695
(evenfn exp-form var))
1696
(let ((ans (zto%pi2 rat-form 'yy)))
1697
(cond (ans (m*t '((rat) 1. 2.) ans))
1699
((and (alike1 b half%pi)
1700
(evenfn exp-form var)
1702
(no-err-sub (m+t '$%pi (m*t -1. var))
1704
(let ((ans (zto%pi2 rat-form 'yy)))
1705
(cond (ans (m*t '((rat) 1. 4.) ans))
1536
1708
;;; Do integrals of sin and cos. this routine makes sure lower limit
1538
(defun INTSC1 (A B E)
1710
(defun intsc1 (a b e)
1539
1711
(let ((limit-diff (m+ b (m* -1 a)))
1542
(sin-cos-recur t)) ;recursion stopper
1543
(PROG (ANS D NZP2 L)
1544
(COND ((OR (NOT (MNUMP (M// limit-diff '$%PI)))
1545
(NOT (PERIOD %PI2 E VAR)))
1548
(setq e (MAXIMA-SUBSTITUTE (m+ a var) var e))
1550
(setq b limit-diff)))
1714
(*sin-cos-recur* t)) ;recursion stopper
1715
(prog (ans d nzp2 l)
1716
(cond ((or (not (mnump (m// limit-diff '$%pi)))
1717
(not (period %pi2 e var)))
1720
(setq e (maxima-substitute (m+ a var) var e))
1722
(setq b limit-diff)))
1551
1723
;;;Multiples of 2*%pi in limits.
1552
(COND ((eq (ask-integer (SETQ D (let (($float nil))
1553
(m// limit-diff %pi2))) '$integer)
1555
(SETQ ANS (m* D (COND ((SETQ ANS (INTSC E %PI2 VAR))
1557
(T (RETURN NIL)))))))
1558
(COND ((RATGREATERP %PI2 B)
1559
(RETURN (INTSC E B VAR)))
1565
(SETQ limit-diff 0.)
1569
(m*t -1. (COND ((SETQ ANS (INTSC E (CDR L) VAR))
1572
(SETQ NZP2 (m+ (CAR B) (m- (CAR L))))
1573
OUT (SETQ ANS (add* (COND ((ZEROP1 NZP2) 0.)
1574
((SETQ ANS (INTSC E %PI2 VAR))
1577
(COND ((ZEROP1 (CDR B)) 0.)
1578
((SETQ ANS (INTSC E (CDR B) VAR))
1584
(DEFUN INTSC (SC B VAR)
1585
(COND ((EQ ($sign B) '$neg)
1586
(SETQ B (m*t -1. B))
1587
(SETQ SC (m* -1. (SUBIN (m*t -1. VAR) SC)))))
1588
(SETQ SC (PARTITION SC VAR 1))
1589
(COND ((SETQ B (INTSC0 (CDR SC) B VAR))
1590
(m* (resimplify (CAR SC)) B))))
1592
(DEFUN INTSC0 (SC B VAR)
1724
(cond ((eq (ask-integer (setq d (let (($float nil))
1725
(m// limit-diff %pi2))) '$integer)
1727
(setq ans (m* d (cond ((setq ans (intsc e %pi2 var))
1729
(t (return nil)))))))
1730
(cond ((ratgreaterp %pi2 b)
1731
(return (intsc e b var)))
1737
(setq limit-diff 0.)
1741
(m*t -1. (cond ((setq ans (intsc e (cdr l) var))
1744
(setq nzp2 (m+ (car b) (m- (car l))))
1745
out (setq ans (add* (cond ((zerop1 nzp2) 0.)
1746
((setq ans (intsc e %pi2 var))
1749
(cond ((zerop1 (cdr b)) 0.)
1750
((setq ans (intsc e (cdr b) var))
1756
(defun intsc (sc b var)
1757
(cond ((eq ($sign b) '$neg)
1758
(setq b (m*t -1. b))
1759
(setq sc (m* -1. (subin (m*t -1. var) sc)))))
1760
(setq sc (partition sc var 1))
1761
(cond ((setq b (intsc0 (cdr sc) b var))
1762
(m* (resimplify (car sc)) b))))
1764
(defun intsc0 (sc b var)
1593
1765
(let ((nn* (scprod sc))
1595
(COND (NN* (COND ((ALIKE1 B HALF%PI)
1596
(BYGAMMA (CAR NN*) (CADR NN*)))
1598
(cond ((eq (real-branch (cadr nn*) -1.) '$yes)
1599
(m* (m+ 1. (m^ -1 (CADR NN*)))
1600
(BYGAMMA (car nn*) (cadr nn*))))))
1602
(COND ((or (AND (eq (ask-integer (CAR NN*) '$even)
1604
(eq (ask-integer (CADR NN*) '$even)
1606
(and (ratnump (car nn*))
1607
(eq (real-branch (car nn*) -1.)
1609
(ratnump (cadr nn*))
1610
(eq (real-branch (cadr nn*) -1.)
1612
(m* 4. (BYGAMMA (car nn*) (cadr nn*))))
1613
((or (eq (ask-integer (car nn*) '$odd) '$yes)
1614
(eq (ask-integer (cadr nn*) '$odd) '$yes))
1617
((ALIKE1 B HALF%PI3)
1618
(m* (m+ 1. (m^ -1 (CADR NN*)) (m^ -1 (m+l NN*)))
1619
(bygamma (car nn*) (cadr nn*))))))
1620
(t (cond ((AND (OR (EQ B '$%PI)
1623
(SETQ DN* (SCRAT SC B)))
1625
((SETQ NN* (ANTIDERIV SC))
1626
(sin-cos-intsubs nn* var 0. b))
1767
(cond (nn* (cond ((alike1 b half%pi)
1768
(bygamma (car nn*) (cadr nn*)))
1770
(cond ((eq (real-branch (cadr nn*) -1.) '$yes)
1771
(m* (m+ 1. (m^ -1 (cadr nn*)))
1772
(bygamma (car nn*) (cadr nn*))))))
1774
(cond ((or (and (eq (ask-integer (car nn*) '$even)
1776
(eq (ask-integer (cadr nn*) '$even)
1778
(and (ratnump (car nn*))
1779
(eq (real-branch (car nn*) -1.)
1781
(ratnump (cadr nn*))
1782
(eq (real-branch (cadr nn*) -1.)
1784
(m* 4. (bygamma (car nn*) (cadr nn*))))
1785
((or (eq (ask-integer (car nn*) '$odd) '$yes)
1786
(eq (ask-integer (cadr nn*) '$odd) '$yes))
1789
((alike1 b half%pi3)
1790
(m* (m+ 1. (m^ -1 (cadr nn*)) (m^ -1 (m+l nn*)))
1791
(bygamma (car nn*) (cadr nn*))))))
1792
(t (cond ((and (or (eq b '$%pi)
1795
(setq dn* (scrat sc b)))
1797
((setq nn* (antideriv sc))
1798
(sin-cos-intsubs nn* var 0. b))
1629
1801
;;;Is careful about substitution of limits where the denominator may be zero
1630
1802
;;;because of various assumptions made.
1639
1811
(denom (pdis (cddr rat-exp))))
1640
1812
(cond ((not (equal (intsubs num ll ul) 0.))
1641
1813
(intsubs exp ll ul))
1814
;; Why do we want to return zero when the denom is not zero?
1815
;; That doesn't seem to make sense to me (rtoy). Checking
1816
;; for a zero denominator makes sense, but what we should
1817
;; return in that case? 0 seems like a bad choice. $inf or
1818
;; $undefined seem like better choices. Or maybe just
1819
;; signaling an error?
1642
1821
((not (equal ($asksign denom) '$zero))
1823
((equal ($asksign denom) '$zero)
1644
1825
(t (let (($%piargs ()))
1645
1826
(intsubs exp ll ul))))))
1648
(let ((great-minus-1 #'(lambda (temp)
1649
(ratgreaterp temp -1)))
1652
((SETQ M (Powerofx E `((%SIN) ,VAR) great-minus-1 VAR))
1654
((SETQ N (Powerofx E `((%COS) ,VAR) great-minus-1 VAR))
1658
(OR (SETQ M (Powerofx (CADR E) `((%SIN) ,VAR) great-minus-1 VAR))
1659
(SETQ N (Powerofx (CADR E) `((%COS) ,VAR) great-minus-1 VAR)))
1662
(SETQ M (Powerofx (CADDR E) `((%SIN) ,VAR) great-minus-1 VAR)))
1663
(T (SETQ N (Powerofx (CADDR E) `((%COS) ,VAR) great-minus-1 VAR))))
1829
(let ((great-minus-1 #'(lambda (temp)
1830
(ratgreaterp temp -1)))
1833
((setq m (powerofx e `((%sin) ,var) great-minus-1 var))
1835
((setq n (powerofx e `((%cos) ,var) great-minus-1 var))
1839
(or (setq m (powerofx (cadr e) `((%sin) ,var) great-minus-1 var))
1840
(setq n (powerofx (cadr e) `((%cos) ,var) great-minus-1 var)))
1843
(setq m (powerofx (caddr e) `((%sin) ,var) great-minus-1 var)))
1844
(t (setq n (powerofx (caddr e) `((%cos) ,var) great-minus-1 var))))
1668
1849
(defun real-branch (exponent value)
1669
;;;Says wether (m^t value exponent) has at least one real branch.
1670
;;;Only works for values of 1 and -1 now.
1671
;;;Returns $yes $no $unknown.
1672
(cond ((equal value 1.)
1674
((eq (ask-integer exponent '$integer) '$yes)
1677
(cond ((eq ($oddp (caddr exponent)) t)
1682
(DEFUN BYGAMMA (M N)
1683
(let ((one-half (m//t 1. 2.)))
1684
(m* one-half `(($BETA) ,(m* one-half (m+t 1. M))
1685
,(m* one-half (m+t 1. N))))))
1687
;Seems like Guys who call this don't agree on what it should return.
1688
(DEFUN Powerofx (E X P VAR)
1689
(SETQ E (COND ((NOT (AMONG VAR E)) NIL)
1694
(NOT (AMONG VAR (CADDR E))))
1696
(COND ((NULL E) NIL)
1699
(DECLARE-TOP(SPECIAL L C K))
1701
(COMMENT (THE FOLLOWING FUNC IS NOT COMPLETE))
1705
; (COND ((NOTINVOLVE E '(%SIN %COS %TAN %LOG))
1706
; (cond ((SETQ ANS (BATAP E))
1708
; (cond ((AND (NOTINVOLVE E '(%SIN %COS %TAN))
1710
; (COND ((SETQ ANS (BATAP (M// E `((%LOG) ,VAR))))
1711
; (SETQ K NN* L DN*)
1713
; (m+ (subfunmake '$PSI '(0) (list K))
1714
; (m* -1. (subfunmake '$PSI
1718
; (return ans)))))))
1723
(COND ((ATOM E) NIL)
1726
(OR (AND (SETQ K (FINDP (CADR E)))
1727
(SETQ C (BXM (CADDR E) (POLYINX (CADDR E) VAR NIL))))
1728
(AND (SETQ K (FINDP (CADDR E)))
1729
(SETQ C (BXM (CADR E) (POLYINX (CADR E) VAR NIL))))))
1731
((SETQ C (BXM E (POLYINX E VAR NIL)))
1736
; (COND ((NOT (BATA0 E)) (RETURN NIL))
1737
; ((AND (EQUAL -1. (CADDDR C))
1738
; (EQ ($askSIGN (SETQ K (m+ 1. K)))
1740
; (EQ ($askSIGN (SETQ L (m+ 1. (CAR C))))
1743
; (m^ UL (CADDR C)))
1745
; (EQ ($askSIGN (SETQ C (CADDR C))) '$pos))
1746
; (RETURN (M// (m* (m^ UL (m+t K (m* C (m+t -1. L))))
1747
; `(($BETA) ,(SETQ NN* (M// K C))
1752
;;; Integrals of the form i(log(x)^m*x^k*(a+b*x^n)^l,x,0,ul) with
1753
;;; ul>0, b<0, a=-b*ul^n, k>-1, l>-1, n>0, m a nonnegative integer.
1754
;;; These integrals are essentially partial derivatives of the
1755
;;; Beta function (i.e. the Eulerian integral of the first kind).
1756
;;; Note that, currently, with the default setting intanalysis:true,
1757
;;; this function might not even be called for some of these integrals.
1758
;;; However, this can be palliated by setting intanalysis:false.
1850
;; Says wether (m^t value exponent) has at least one real branch.
1851
;; Only works for values of 1 and -1 now. Returns $yes $no
1853
(cond ((equal value 1.)
1855
((eq (ask-integer exponent '$integer) '$yes)
1858
(cond ((eq ($oddp (caddr exponent)) t)
1863
;; Compute beta((m+1)/2,(n+1)/2)/2.
1864
(defun bygamma (m n)
1865
(let ((one-half (m//t 1. 2.)))
1866
(m* one-half `(($beta) ,(m* one-half (m+t 1. m))
1867
,(m* one-half (m+t 1. n))))))
1869
;;Seems like Guys who call this don't agree on what it should return.
1870
(defun powerofx (e x p var)
1871
(setq e (cond ((not (among var e)) nil)
1876
(not (among var (caddr e))))
1878
(cond ((null e) nil)
1882
(comment (the following func is not complete))
1886
;; (COND ((NOTINVOLVE E '(%SIN %COS %TAN %LOG))
1887
;; (cond ((SETQ ANS (BATAP E))
1889
;; (cond ((AND (NOTINVOLVE E '(%SIN %COS %TAN))
1891
;; (COND ((SETQ ANS (BATAP (M// E `((%LOG) ,VAR))))
1892
;; (SETQ K NN* L DN*)
1893
;; (SETQ ANS (m* ANS
1894
;; (m+ (subfunmake '$PSI '(0) (list K))
1895
;; (m* -1. (subfunmake '$PSI
1899
;; (return ans)))))))
1903
;; Check e for an expression of the form x^kk*(b*x^n+a)^l. If it
1904
;; matches, Return the two values kk and (list l a n b).
1907
(cond ((atom e) nil)
1910
(or (and (setq k (findp (cadr e)))
1911
(setq c (bxm (caddr e) (polyinx (caddr e) var nil))))
1912
(and (setq k (findp (caddr e)))
1913
(setq c (bxm (cadr e) (polyinx (cadr e) var nil))))))
1915
((setq c (bxm e (polyinx e var nil)))
1922
;; (COND ((NOT (BATA0 E)) (RETURN NIL))
1923
;; ((AND (EQUAL -1. (CADDDR C))
1924
;; (EQ ($askSIGN (SETQ K (m+ 1. K)))
1926
;; (EQ ($askSIGN (SETQ L (m+ 1. (CAR C))))
1929
;; (m^ UL (CADDR C)))
1930
;; (SETQ E (CADR C))
1931
;; (EQ ($askSIGN (SETQ C (CADDR C))) '$pos))
1932
;; (RETURN (M// (m* (m^ UL (m+t K (m* C (m+t -1. L))))
1933
;; `(($BETA) ,(SETQ NN* (M// K C))
1938
;; Integrals of the form i(log(x)^m*x^k*(a+b*x^n)^l,x,0,ul). There
1939
;; are two cases to consider: One case has ul>0, b<0, a=-b*ul^n, k>-1,
1940
;; l>-1, n>0, m a nonnegative integer. The second case has ul=inf, l < 0.
1942
;; These integrals are essentially partial derivatives of the Beta
1943
;; function (i.e. the Eulerian integral of the first kind). Note
1944
;; that, currently, with the default setting intanalysis:true, this
1945
;; function might not even be called for some of these integrals.
1946
;; However, this can be palliated by setting intanalysis:false.
1760
1948
(defun zto1 (e)
1761
1949
(when (or (mtimesp e) (mexptp e))
1762
(let ((m 0) (log (list '(%log) var)))
1763
(flet ((set-m (p) (setq m p)))
1764
(find-if #'(lambda (fac) (powerofx fac log #'set-m var)) (cdr e)))
1951
(log (list '(%log) var)))
1954
(find-if #'(lambda (fac)
1955
(powerofx fac log #'set-m var))
1765
1957
(when (and (freeof var m)
1766
1958
(eq (ask-integer m '$integer) '$yes)
1767
1959
(not (eq ($asksign m) '$neg)))
1768
1960
(setq e (m//t e (list '(mexpt) log m)))
1769
(multiple-value-bind
1770
(k/n l n b) (batap-new e)
1772
(let ((beta (simplify (list '($beta) k/n l)))
1773
(m (if (eq ($asksign m) '$zero) 0 m)))
1774
;; The result looks like B(k/n,l) ( ... ).
1775
;; Perhaps, we should just $factor, instead of
1776
;; pulling out beta like this.
1782
(m^t (m-t b) (m1-t l))
1783
(m^t ul (m*t n (m1-t l)))
1784
(m^t n (m-t (m1+t m)))
1785
($at ($diff (m*t (m^t ul (m*t n var)) (list '($beta) var l))
1786
var m) (list '(mequal) var k/n)))
1963
(multiple-value-bind (kk s d r cc)
1965
;; We have i(x^kk/(d+cc*x^r)^s,x,0,inf) =
1966
;; beta(aa,bb)/(cc^aa*d^bb*r). Compute this, and then
1967
;; differentiate it m times to get the log term
1970
(let* ((aa (div (add 1 var) r))
1972
(m (if (eq ($asksign m) '$zero)
1975
(let ((res (div `(($beta) ,aa ,bb)
1979
($at ($diff res var m)
1980
(list '(mequal) var kk)))))))
1982
(multiple-value-bind
1983
(k/n l n b) (batap-new e)
1985
(let ((beta (simplify (list '($beta) k/n l)))
1986
(m (if (eq ($asksign m) '$zero) 0 m)))
1987
;; The result looks like B(k/n,l) ( ... ).
1988
;; Perhaps, we should just $factor, instead of
1989
;; pulling out beta like this.
1995
(m^t (m-t b) (m1-t l))
1996
(m^t ul (m*t n (m1-t l)))
1997
(m^t n (m-t (m1+t m)))
1998
($at ($diff (m*t (m^t ul (m*t n var))
1999
(list '($beta) var l))
2001
(list '(mequal) var k/n)))
1790
2005
;;; If e is of the form given below, make the obvious change
1796
2011
;;; of course would give wrong results.
1798
2013
(defun batap-new (e)
1800
(declare (special k c))
1803
(multiple-value-bind
1804
;; e=x^k*(a+b*x^n)^l
1805
(l a n b) (values-list c)
1806
(when (and (freeof var k) (freeof var n) (freeof var l)
1807
(alike1 a (m-t (m*t b (m^t ul n))))
1808
(eq ($asksign b) '$neg)
1809
(eq ($asksign (setq k (m1+t k))) '$pos)
1810
(eq ($asksign (setq l (m1+t l))) '$pos)
1811
(eq ($asksign n) '$pos))
1812
(values (m//t k n) l n b))))))
1818
(COND ((NOT (OR (EQUAL LL 0) (EQ LL '$MINF)))
1819
(SETQ E (SUBIN (m+ LL VAR) E))))
1820
(COND ((NOT (BATA0 E)) (RETURN NIL))
1821
((AND (RATGREATERP (SETQ AL (CADDR C)) 0.)
1822
(EQ ($askSIGN (SETQ K (M// (m+ 1. K)
1825
(RATGREATERP (SETQ L (m* -1. (CAR C)))
1827
(EQ ($askSIGN (m* (SETQ D (CADR C))
1828
(SETQ C (CADDDR C))))
1830
(SETQ L (m+ L (m*t -1. K)))
1831
(RETURN (M// `(($BETA) ,K ,L)
1832
(mul* AL (m^ C K) (m^ D L))))))))
1834
(DECLARE-TOP(UNSPECIAL L C K))
1836
(DEFUN GAMMA1 (C A B D)
1838
(m^ (m* B (m^ A (SETQ C (M// (m+t C 1.) B))))
2015
(multiple-value-bind (k c)
2018
;; e=x^k*(a+b*x^n)^l
2019
(destructuring-bind (l a n b)
2021
(when (and (freeof var k)
2024
(alike1 a (m-t (m*t b (m^t ul n))))
2025
(eq ($asksign b) '$neg)
2026
(eq ($asksign (setq k (m1+t k))) '$pos)
2027
(eq ($asksign (setq l (m1+t l))) '$pos)
2028
(eq ($asksign n) '$pos))
2029
(values (m//t k n) l n b))))))
2032
;; Wang p. 71 gives the following formula for a beta function:
2034
;; integrate(x^(k-1)/(c*x^r+d)^s,x,0,inf)
2035
;; = beta(a,b)/(c^a*d^b*r)
2037
;; where a = k/r > 0, b = s - a > 0, s > k > 0, r > 0, c*d > 0.
2039
;; This function matches this and returns k-1, d, r, c, a, b. And
2040
;; also checks that all the conditions hold. If not, NIL is returned.
2042
(defun batap-inf (e)
2043
(multiple-value-bind (k c)
2046
(destructuring-bind (l d r cc)
2048
(let* ((s (mul -1 l))
2052
(when (and (freeof var k)
2055
(eq ($asksign kk) '$pos)
2056
(eq ($asksign a) '$pos)
2057
(eq ($asksign b) '$pos)
2058
(eq ($asksign (sub s k)) '$pos)
2059
(eq ($asksign r) '$pos)
2060
(eq ($asksign (mul cc d)) '$pos))
2061
(values k s d r cc)))))))
2064
;; Handles beta integrals.
2066
(cond ((not (or (equal ll 0)
2068
(setq e (subin (m+ ll var) e))))
2069
(multiple-value-bind (k c)
2074
(destructuring-bind (l d al c)
2076
;; e = x^k*(d+c*x^al)^l.
2077
(let ((new-k (m// (m+ 1 k) al)))
2078
(when (and (ratgreaterp al 0.)
2079
(eq ($asksign new-k) '$pos)
2080
(ratgreaterp (setq l (m* -1. l))
2082
(eq ($asksign (m* d c))
2084
(setq l (m+ l (m*t -1. new-k)))
2085
(m// `(($beta) ,new-k ,l)
2086
(mul* al (m^ c new-k) (m^ d l))))))))))
2089
;; Compute exp(d)*gamma((c+1)/b)/b/a^((c+1)/b). In essence, this is
2090
;; the value of integrate(x^c*exp(d-a*x^b),x,0,inf).
2091
(defun gamma1 (c a b d)
2093
(m^ (m* b (m^ a (setq c (m// (m+t c 1.) b))))
1842
(DEFUN ZTO%PI2 (GRAND VAR)
1843
(let ((result (UNITCIR ($RATSIMP (M// GRAND VAR)) VAR)))
1844
(cond (result (sratsimp (m* (m- '$%I) result)))
2097
(defun zto%pi2 (grand var)
2098
(let ((result (unitcir ($ratsimp (m// grand var)) var)))
2099
(cond (result (sratsimp (m* (m- '$%i) result)))
1847
(DEFUN UNITCIR (GRAND VAR)
1849
(let ((result (PRINCIP (RES NN* DN* #'(LAMBDA (PT)
1850
(RATGREATERP 1 (CABS PT)))
1852
(ALIKE1 1 (CABS PT)))))))
1853
(cond (result (m* '$%PI result))
2102
(defun unitcir (grand var)
2104
(let ((result (princip (res nn* dn* #'(lambda (pt)
2105
(ratgreaterp 1 (cabs pt)))
2107
(alike1 1 (cabs pt)))))))
2108
(cond (result (m* '$%pi result))
1857
(DEFUN LOGX1 (EXP LL UL)
2112
(defun logx1 (exp ll ul)
1858
2113
(let ((arg nil))
1860
((AND (NOTINVOLVE EXP '(%SIN %COS %TAN %ATAN %ASIN %ACOS))
1861
(SETQ ARG (INVOLVE EXP '(%LOG))))
1863
(COND ((RATGREATERP 1. LL)
1864
(COND ((NOT (EQ UL '$INF))
1865
(INTCV1 (m^t '$%E (m- VAR)) () (m- `((%LOG) ,VAR))))
1866
(T (INTCV1 (m^t '$%E VAR) () `((%LOG) ,VAR)))))))
1867
(t (INTCV ARG NIL nil)))))))
2115
((and (notinvolve exp '(%sin %cos %tan %atan %asin %acos))
2116
(setq arg (involve exp '(%log))))
2118
(cond ((ratgreaterp 1. ll)
2119
(cond ((not (eq ul '$inf))
2120
(intcv1 (m^t '$%e (m- var)) () (m- `((%log) ,var))))
2121
(t (intcv1 (m^t '$%e var) () `((%log) ,var)))))))
2122
(t (intcv arg nil nil)))))))
1872
(COND ((ATOM E) NIL)
1873
((AND (OR (EQ (CAAR E) '%SIN)
1874
(EQ (CAAR E) '%COS))
1876
(SETQ E (BX**N (CADR E))))
1877
(COND ((EQUAL (CAR E) 1.) '$IND)
1878
((ZEROP (SETQ S (let ((sign ($askSIGN (CADR E))))
1879
(cond ((eq sign '$pos) 1)
1880
((eq sign '$neg) -1)
1881
((eq sign '$zero) 0)))))
1883
((not (EQ ($askSIGN (m+ -1 (CAR E))) '$pos))
1885
(t (SETQ G (GAMMA1 0. (m* S (CADR E)) (CAR E) 0.))
1886
(SETQ E (m* G `((,IND) ,(M// HALF%PI (CAR E)))))
1887
(m* (COND ((AND (EQ IND '%SIN)
2125
;; Wang 81-83. Unfortunately, the pdf version has page 82 as all
2126
;; black, so here is, as best as I can tell, what Wang is doing.
2127
;; Fortunately, p. 81 has the necessary hints.
2129
;; First consider integrate(exp(%i*k*x^n),x) around the closed contour
2130
;; consisting of the real axis from 0 to R, the arc from the angle 0
2131
;; to %pi/(2*n) and the ray from the arc back to the origin.
2133
;; There are no poles in this region, so the integral must be zero.
2134
;; But consider the integral on the three parts. The real axis is the
2135
;; integral we want. The return ray is
2137
;; exp(%i*%pi/2/n) * integrate(exp(%i*k*(t*exp(%i*%pi/2/n))^n),t,R,0)
2138
;; = exp(%i*%pi/2/n) * integrate(exp(%i*k*t^n*exp(%i*%pi/2)),t,R,0)
2139
;; = -exp(%i*%pi/2/n) * integrate(exp(-k*t^n),t,0,R)
2141
;; As R -> infinity, this last integral is gamma(1/n)/k^(1/n)/n.
2143
;; We assume the integral on the circular arc approaches 0 as R ->
2144
;; infinity. (Need to prove this.)
2148
;; integrate(exp(%i*k*t^n),t,0,inf)
2149
;; = exp(%i*%pi/2/n) * gamma(1/n)/k^(1/n)/n.
2151
;; Equating real and imaginary parts gives us the desired results:
2153
;; integrate(cos(k*t^n),t,0,inf) = G * cos(%pi/2/n)
2154
;; integrate(sin(k*t^n),t,0,inf) = G * sin(%pi/2/n)
2156
;; where G = gamma(1/n)/k^(1/n)/n.
2160
(cond ((atom e) nil)
2161
((and (or (eq (caar e) '%sin)
2162
(eq (caar e) '%cos))
2164
(setq e (bx**n (cadr e))))
2165
;; Ok, we have cos(b*x^n) or sin(b*x^n), and we set e = (a
2166
;; b n) where the arg of the trig function is b*x^n+a.
2167
(cond ((equal (car e) 1.)
2170
((zerop (setq s (let ((sign ($asksign (cadr e))))
2171
(cond ((eq sign '$pos) 1)
2172
((eq sign '$neg) -1)
2173
((eq sign '$zero) 0)))))
2174
;; s is the sign of b. Give up if it's zero.
2176
((not (eq ($asksign (m+ -1 (car e))) '$pos))
2177
;; Give up if a-1 <= 0
2180
;; We can apply our formula now. g = gamma(1/n)/n/b^(1/n)
2181
(setq g (gamma1 0. (m* s (cadr e)) (car e) 0.))
2182
(setq e (m* g `((,ind) ,(m// half%pi (car e)))))
2183
(m* (cond ((and (eq ind '%sin)
1894
(COMMENT THIS IS THE SECOND PART OF THE DEFINITE INTEGRAL PACKAGE)
1896
(DECLARE-TOP(SPECIAL VAR PLM* PL* RL* PL*1 RL*1))
1898
(DEFUN P*LOGNXP (A S)
1900
(COND ((NOT (AMONG '%LOG A))
2190
(comment this is the second part of the definite integral package)
2192
(declare-top(special var plm* pl* rl* pl*1 rl*1))
2194
(defun p*lognxp (a s)
2196
(cond ((not (among '%log a))
1902
((AND (POLYINX (SETQ B (MAXIMA-SUBSTITUTE 1. `((%LOG) ,VAR) A))
1904
(eq ($sign (m+ S (M+ 1 (DEG B))))
2198
((and (polyinx (setq b (maxima-substitute 1. `((%log) ,var) a))
2200
(eq ($sign (m+ s (m+ 1 (deg b))))
1907
(SETQ A (LOGNXP ($RATSIMP (M// A B)))))
1911
(COND ((ATOM A) NIL)
1912
((AND (EQ (CAAR A) '%LOG)
1913
(EQ (CADR A) VAR)) 1.)
1919
(COMMENT CHECK THE FOLLOWING FUNCTION FOR UNUSED PROG VAR A)
1921
(DEFUN LOGCPI0 (N D)
1923
(SETQ PL (POLELIST D #'UPPERHALF #'(LAMBDA (J)
1924
(COND ((ZEROP1 J) NIL)
1925
((equal ($imagpart j) 0)
1929
(SETQ FACTORS (CAR PL)
1931
(COND ((OR (CADR PL)
1933
(SETQ DP (SDIFF D VAR))))
1934
(COND ((SETQ PLM* (CAR PL))
1935
(SETQ RLM* (RESIDUE N (COND (LEADCOEF FACTORS)
1938
(COND ((SETQ PL* (CADR PL))
1939
(SETQ RL* (RES1 N DP PL*))))
1940
(COND ((SETQ PL*1 (CADDR PL))
1941
(SETQ RL*1 (RES1 N DP PL*1))))
1942
(RETURN (m*t (m//t 1. 2.)
1945
(LIST (COND ((SETQ NN* (APPEND RL* RLM*))
1947
(COND (RL*1 (M+L RL*1))))))))))
2203
(setq a (lognxp ($ratsimp (m// a b)))))
2207
(cond ((atom a) nil)
2208
((and (eq (caar a) '%log)
2209
(eq (cadr a) var)) 1.)
2215
(comment check the following function for unused prog var a)
2217
(defun logcpi0 (n d)
2219
(setq pl (polelist d #'upperhalf #'(lambda (j)
2220
(cond ((zerop1 j) nil)
2221
((equal ($imagpart j) 0)
2225
(setq factors (car pl)
2227
(cond ((or (cadr pl)
2229
(setq dp (sdiff d var))))
2230
(cond ((setq plm* (car pl))
2231
(setq rlm* (residue n (cond (leadcoef factors)
2234
(cond ((setq pl* (cadr pl))
2235
(setq rl* (res1 n dp pl*))))
2236
(cond ((setq pl*1 (caddr pl))
2237
(setq rl*1 (res1 n dp pl*1))))
2238
(return (m*t (m//t 1. 2.)
2241
(list (cond ((setq nn* (append rl* rlm*))
2243
(cond (rl*1 (m+l rl*1))))))))))
1949
(DEFUN LOGNX2 (NN DN PL RL)
2245
(defun lognx2 (nn dn pl rl)
1950
2246
(do ((pl pl (cdr pl))
1951
2247
(rl rl (cdr rl))
1954
2250
(null rl)) ans)
1955
(SETQ ANS (CONS (m* DN (car rl) (m^ `((%PLOG) ,(car pl)) NN))
2251
(setq ans (cons (m* dn (car rl) (m^ `((%plog) ,(car pl)) nn))
1958
(DEFUN LOGCPJ (N D I)
1960
(COND (PLM* (LIST (mul* (m*t '$%i %pi2)
1962
(RESIDUE (m* (m^ `((%PLOG) ,VAR) I)
1966
(LOGNX2 I (m*t '$%i %pi2) PL* RL*)
1967
(LOGNX2 I %P%I PL*1 RL*1)))
1969
(T (SIMPLIFY (M+L N)))))
2254
(defun logcpj (n d i)
2256
(cond (plm* (list (mul* (m*t '$%i %pi2)
2258
(residue (m* (m^ `((%plog) ,var) i)
2262
(lognx2 i (m*t '$%i %pi2) pl* rl*)
2263
(lognx2 i %p%i pl*1 rl*1)))
2265
(t (simplify (m+l n)))))
1972
2268
;;should replace these references to *i* and *j* to symbol-value arrays.
1973
2269
;;here and in SUMI, and LOGCPI. These are the only references in this file.
1974
2270
;;I did change I to *I*
1976
#-cl ;in case other lisps don't understand internal declares.
2272
#-cl ;in case other lisps don't understand internal declares.
1977
2273
(declare-top(special *i* *j*))
1979
(DEFUN LOG*RAT (N D M)
1980
(PROG (LEADCOEF FACTORS C PLM* PL* RL* PL*1 RL*1 RLM*)
1981
(declare (special *i* *j*))
1982
; (ARRAY *I* T (M+ 1 M))
1983
; (ARRAY *J* T (M+ 1 M))
1984
(setq *i* (*ARRAY nil T (M+ 1 M)))
1985
(setq *j* (*ARRAY nil T (M+ 1 M)))
1987
(STORE (aref *J* C) 0.)
1988
(do ((c 0. (m+ 1 c)))
1990
(return (logcpi n d m)))
1991
(STORE (aref *I* C) (LOGCPI N D C))
1992
(STORE (aref *J* C) (LOGCPJ N FACTORS C)))))
2275
;; Handle integral(n(x)/d(x)*log(x)^m,x,0,inf). n and d are
2277
(defun log*rat (n d m)
2278
(prog (leadcoef factors c plm* pl* rl* pl*1 rl*1 rlm*)
2279
(declare (special *i* *j*))
2280
;; (ARRAY *I* T (M+ 1 M))
2281
;; (ARRAY *J* T (M+ 1 M))
2282
(setq *i* (*array nil t (m+ 1 m)))
2283
(setq *j* (*array nil t (m+ 1 m)))
2285
(store (aref *j* c) 0.)
2286
(do ((c 0. (m+ 1 c)))
2288
(return (logcpi n d m)))
2289
(store (aref *i* c) (logcpi n d c))
2290
(store (aref *j* c) (logcpj n factors c)))))
1994
(DEFUN LOGCPI (N D C)
2292
(defun logcpi (n d c)
1995
2293
(declare (special *j*))
1998
(T (m* '((RAT) 1. 2.)
1999
(m+ (aref *J* C) (m* -1. (SUMI C)))))))
2296
(t (m* '((rat) 1. 2.)
2297
(m+ (aref *j* c) (m* -1. (sumi c)))))))
2002
2300
(declare (special *i*))
2003
2301
(do ((k 1. (m+ 1 k))
2007
(SETQ ANS (CONS (mul* ($MAKEGAMMA `((%BINOMIAL) ,C ,K))
2010
(aref *I* (m+ C (m- K))))
2305
(setq ans (cons (mul* ($makegamma `((%binomial) ,c ,k))
2308
(aref *i* (m+ c (m- k))))
2013
(DEFUN FAN (P M A N B)
2311
(defun fan (p m a n b)
2014
2312
(let ((povern (m// p n))
2015
2313
(ab (m// a b)))
2017
((OR (eq (ask-integer POVERN '$integer) '$yes)
2018
(NOT (equal ($imagpart ab) 0))) ())
2019
(t (let ((IND ($askSIGN AB)))
2020
(COND ((eq IND '$zero) NIL)
2021
((eq ind '$neg) nil)
2022
((NOT (RATGREATERP M POVERN)) NIL)
2024
($MAKEGAMMA `((%BINOMIAL) ,(m+ -1. M (m- POVERN))
2026
`((mabs) ,(m^ A (m+ POVERN (m- m)))))
2029
`((%SIN) ,(m*t '$%PI POVERN)))))))))))
2315
((or (eq (ask-integer povern '$integer) '$yes)
2316
(not (equal ($imagpart ab) 0))) ())
2317
(t (let ((ind ($asksign ab)))
2318
(cond ((eq ind '$zero) nil)
2319
((eq ind '$neg) nil)
2320
((not (ratgreaterp m povern)) nil)
2322
($makegamma `((%binomial) ,(m+ -1. m (m- povern))
2324
`((mabs) ,(m^ a (m+ povern (m- m)))))
2327
`((%sin) ,(m*t '$%pi povern)))))))))))
2032
2330
;;Makes a new poly such that np(x)-np(x+2*%i*%pi)=p(x).
2035
2333
;;of the varibale of integration.
2036
2334
;;Can probably be made simpler now.
2039
2337
(let ((n (deg p)) (ans ()) (varlist ()) (gp ()) (cl ()) (zz ()))
2040
(SETQ ANS (GENPOLY (M+ 1 N))) ;Make poly with gensyms of 1 higher deg.
2041
(setq cl (cdr ans)) ;Coefficient list
2042
(SETQ VARLIST (APPEND CL (LIST VAR))) ;Make VAR most important.
2043
(SETQ GP (CAR ANS)) ;This is the poly with gensym coeffs.
2338
(setq ans (genpoly (m+ 1 n))) ;Make poly with gensyms of 1 higher deg.
2339
(setq cl (cdr ans)) ;Coefficient list
2340
(setq varlist (append cl (list var))) ;Make VAR most important.
2341
(setq gp (car ans)) ;This is the poly with gensym coeffs.
2044
2342
;;;Now, poly(x)-poly(x+2*%i*%pi)=p(x), P is the original poly.
2045
(SETQ ANS (m+ GP (SUBIN (m+t (m*t '$%i %pi2) VAR) (m- GP)) (m- P)))
2047
(SETQ ANS (RATREP* ANS)) ;Rational rep with VAR leading.
2048
(SETQ ZZ (COEFSOLVE N CL (COND ((NOT (EQ (CAADR ANS) ;What is Lead Var.
2049
(GENFIND (CAR ANS) VAR)))
2050
(LIST 0 (CADR ANS)));No VAR in ans.
2051
((CDADR ANS)))));The real Poly.
2343
(setq ans (m+ gp (subin (m+t (m*t '$%i %pi2) var) (m- gp)) (m- p)))
2345
(setq ans (ratrep* ans)) ;Rational rep with VAR leading.
2346
(setq zz (coefsolve n cl (cond ((not (eq (caadr ans) ;What is Lead Var.
2347
(genfind (car ans) var)))
2348
(list 0 (cadr ans))) ;No VAR in ans.
2349
((cdadr ans))))) ;The real Poly.
2052
2350
(if (or (null zz) (null gp))
2054
($SUBSTITUTE ZZ GP))));Substitute Values for gensyms.
2056
(DEFUN COEFSOLVE (N CL E)
2058
(EQL (NCONS (PDIS (PTERM E N))) (CONS (PDIS (PTERM E M)) EQL))
2059
(M (m+ N -1) (m+ M -1)))
2060
((SIGNP L M) (SOLVEX EQL CL NIL NIL))))
2062
(DEFUN RECTZTO%PI2 (P PE D)
2063
(PROG (DP N PL A B C denom-exponential)
2064
(if (not (and (SETQ denom-exponential (CATCH 'PIN%EX (PIN%EX D)))
2065
(%e-integer-coeff pe)
2066
(%e-integer-coeff d)))
2068
(SETQ N (m* (COND ((NULL P) -1.)
2069
(T ($EXPAND (m*t '$%i %pi2 (MAKPOLY P)))))
2073
(SETQ PL (CDR (POLELIST DENOM-EXPONENTIAL #'(LAMBDA (J)
2074
(OR (NOT (equal ($imagpart J) 0))
2075
(EQ ($askSIGN ($realpart J)) '$neg)))
2077
(NOT (EQ ($askSIGN ($realpart J)) '$zero)))))))
2078
(COND ((NULL PL) (RETURN NIL))
2080
(CADDR PL)) (SETQ DP (SDIFF D VAR))))
2081
(COND ((CADR PL) (SETQ B (MAPCAR #'log-imag-0-2%pi (CADR PL)))
2082
(SETQ B (RES1 N DP B))
2086
(let ((temp (MAPCAR #'log-imag-0-2%pi (CADDR PL))))
2087
(setq c (append temp (mapcar #'(lambda (j)
2088
(m+ (m*t '$%i %pi2) j))
2090
(SETQ C (RES1 N DP C))
2094
(let ((poles (mapcar #'log-imag-0-2%pi (caar pl)))
2095
(exp (m// n (subst (m^t '$%e var) 'z* denom-exponential))))
2096
(SETQ A (mapcar #'(lambda (j)
2097
($RESIDUE exp var j))
2101
(RETURN (sRATSIMP (m+ a b (m* '((rat) 1. 2.) c))))))
2352
($substitute zz gp)))) ;Substitute Values for gensyms.
2354
(defun coefsolve (n cl e)
2356
(eql (ncons (pdis (pterm e n))) (cons (pdis (pterm e m)) eql))
2357
(m (m+ n -1) (m+ m -1)))
2358
((signp l m) (solvex eql cl nil nil))))
2360
;; Integrate(p(x)*f(exp(x))/g(exp(x)),x,minf,inf) by applying the
2361
;; transformation y = exp(x) to get
2362
;; integrate(p(log(y))*f(y)/g(y)/y,y,0,inf). This should be handled
2364
(defun log-transform (p pe d)
2365
(let ((new-p (subst (list '(%log) var) var p))
2366
(new-pe (subst var 'z* (catch 'pin%ex (pin%ex pe))))
2367
(new-d (subst var 'z* (catch 'pin%ex (pin%ex d)))))
2368
(defint (div (div (mul new-p new-pe) new-d) var) var 0 ul)))
2370
;; This implements Wang's algorithm in Chapter 5.2, pp. 98-100.
2372
;; This is a very brief description of the algorithm. Basically, we
2373
;; have integrate(R(exp(x))*p(x),x,minf,inf), where R(x) is a rational
2374
;; function and p(x) is a polynomial.
2376
;; We find a polynomial q(x) such that q(x) - q(x+2*%i*%pi) = p(x).
2377
;; Then consider a contour integral of R(exp(z))*q(z) over a
2378
;; rectangular contour. Opposite corners of the rectangle are (-R,
2379
;; 2*%i*%pi) and (R, 0).
2381
;; Wang shows that this contour integral, in the limit, is the
2382
;; integral of R(exp(x))*q(x)-R(exp(x))*q(x+2*%i*%pi), which is
2383
;; exactly the integral we're looking for.
2385
;; Thus, to find the value of the contour integral, we just need the
2386
;; residues of R(exp(z))*q(z). The only tricky part is that we want
2387
;; the log function to have an imaginary part between 0 and 2*%pi
2388
;; instead of -%pi to %pi.
2389
(defun rectzto%pi2 (p pe d)
2390
(prog (dp n pl a b c denom-exponential)
2391
(if (not (and (setq denom-exponential (catch 'pin%ex (pin%ex d)))
2392
(%e-integer-coeff pe)
2393
(%e-integer-coeff d)))
2395
(setq n (m* (cond ((null p) -1.)
2396
(t ($expand (m*t '$%i %pi2 (makpoly p)))))
2400
;; Find the poles of the denominator. denom-exponential is the
2401
;; denominator of R(x).
2403
;; It seems as if polelist returns a list of several items.
2404
;; The first element is a list consisting of the pole and (z -
2405
;; pole). We don't care about this, so we take the rest of the
2406
;; result. I think the second element of the list is an alist
2407
;; consisting of the pole and it's multiplicity. I don't know
2408
;; what the rest of the list is.
2409
(setq pl (cdr (polelist denom-exponential
2411
;; The imaginary part is nonzero,
2412
;; or the realpart is negative.
2413
(or (not (equal ($imagpart j) 0))
2414
(eq ($asksign ($realpart j)) '$neg)))
2416
;; The realpart is not zero.
2417
(not (eq ($asksign ($realpart j)) '$zero)))))))
2418
;; Not sure what this does.
2423
(setq dp (sdiff d var))))
2424
;; Not sure what this does.
2426
(setq b (mapcar #'log-imag-0-2%pi (cadr pl)))
2427
(setq b (res1 n dp b))
2430
;; Not sure what this does either.
2432
(let ((temp (mapcar #'log-imag-0-2%pi (caddr pl))))
2433
(setq c (append temp (mapcar #'(lambda (j)
2434
(m+ (m*t '$%i %pi2) j))
2436
(setq c (res1 n dp c))
2440
;; We have the poles of deonom-exponential, so we need to
2441
;; convert them to the actual pole values for R(exp(x)),
2442
;; by taking the log of the value of poles.
2443
(let ((poles (mapcar #'(lambda (p)
2444
(log-imag-0-2%pi (car p)))
2446
(exp (m// n (subst (m^t '$%e var) 'z* denom-exponential))))
2447
;; Compute the residues at all of these poles and sum
2449
(setq a (mapcar #'(lambda (j)
2450
($residue exp var j))
2454
(return (sratsimp (m+ a b (m* '((rat) 1. 2.) c))))))
2104
2457
(do ((i i (m+ i -1))
2105
2458
(c (gensym) (gensym))
2109
2462
(cons (m+l ans) cl))
2110
(SETQ ANS (CONS (m* C (m^t VAR I)) ANS))
2111
(SETQ CL (CONS C CL))))
2113
(DECLARE-TOP(SPECIAL *FAILFLAG *LHFLAG LHV *INDICATOR CNT *DISCONFLAG))
2463
(setq ans (cons (m* c (m^t var i)) ans))
2464
(setq cl (cons c cl))))
2466
;;(declare-top(special *failflag *lhflag lhv *indicator cnt *disconflag))
2115
2469
(defun %e-integer-coeff (exp)
2116
2470
(cond ((mapatom exp) t)
2117
2471
((and (mexptp exp)
2121
2475
(t (andmapc '%e-integer-coeff (cdr exp)))))
2123
(DEFUN WLINEARPOLY (E VAR)
2124
(COND ((AND (SETQ E (POLYINX E VAR T))
2128
(DECLARE-TOP(SPECIAL E $EXPONENTIALIZE))
2131
(PIN%EX0 (COND ((NOTINVOLVE EXP '(%SINH %COSH %TANH)) EXP)
2132
(T (SETQ EXP (let (($EXPONENTIALIZE t))
2136
(COND ((NOT (AMONG VAR E)) E)
2137
((ATOM E) (THROW 'PIN%EX NIL))
2140
(COND ((EQ (CADDR E) VAR) 'Z*)
2477
;; Check to see if each term in exp that is of the form exp(k*x) has
2478
;; an integer value for k.
2479
(defun %e-integer-coeff (exp)
2480
(cond ((mapatom exp) t)
2482
(eq (cadr exp) '$%e))
2483
(eq (ask-integer ($coeff (caddr exp) var) '$integer)
2485
(t (andmapc '%e-integer-coeff (cdr exp)))))
2487
(defun wlinearpoly (e var)
2488
(cond ((and (setq e (polyinx e var t))
2492
;; Test to see if exp is of the form f(exp(x)), and if so, replace
2493
;; exp(x) with 'z*. That is, basically return f(z*).
2495
(declare (special $exponentialize))
2496
(pin%ex0 (cond ((notinvolve exp '(%sinh %cosh %tanh))
2499
(let (($exponentialize t))
2500
(setq exp ($expand exp)))))))
2503
;; Does e really need to be special here? Seems to be ok without
2504
;; it; testsuite works.
2506
(declare (special e))
2507
(cond ((not (among var e))
2510
(throw 'pin%ex nil))
2513
(cond ((eq (caddr e) var)
2141
2515
((let ((linterm (wlinearpoly (caddr e) var)))
2143
2517
(m* (subin 0 e) (m^t 'z* linterm)))))
2144
(T (THROW 'PIN%EX NIL))))
2145
((mtimesp E) (m*l (MAPCAR #'PIN%EX0 (CDR E))))
2146
((mplusp E) (m+l (MAPCAR #'PIN%EX0 (CDR E))))
2147
(T (THROW 'PIN%EX NIL))))
2149
(DECLARE-TOP (UNSPECIAL E))
2151
(DEFUN P*PIN%EX (ND*)
2152
(SETQ ND* ($FACTOR ND*))
2153
(COND ((POLYINX ND* VAR NIL) (SETQ P* (CONS ND* P*)) T)
2154
((CATCH 'PIN%EX (PIN%EX ND*)) (SETQ PE* (CONS ND* PE*)) T)
2156
(andMAPCAR #'P*PIN%EX (CDR ND*)))))
2159
(COND ((FINDP P) NIL)
2160
((SETQ ND* (BX**N P))
2161
(m^t VAR (CAR ND*)))
2162
((SETQ P (BX**N+A P))
2163
(m* (CADDR P) (m^t VAR (CADR P))))))
2165
(DEFUN FUNCLOGOR%E (E)
2166
(PROG (ANS ARG NVAR R)
2167
(COND ((OR (RATP E VAR)
2168
(INVOLVE E '(%SIN %COS %TAN))
2169
(NOT (SETQ ARG (XOR (AND (SETQ ARG (INVOLVE E '(%LOG)))
2173
AG (SETQ NVAR (COND ((EQ R '%LOG) `((%LOG) ,ARG))
2174
(T (m^t '$%E ARG))))
2175
(SETQ ANS (MAXIMA-SUBSTITUTE (m^t 'YX -1.) (m^t NVAR -1.) (MAXIMA-SUBSTITUTE 'YX nvar E)))
2176
(COND ((NOT (AMONG VAR ANS)) (RETURN (LIST (SUBST VAR 'YX ANS) NVAR)))
2178
(SETQ ARG (FINDSUB ARG)))
2181
(DEFUN DINTBYPART (U V A B)
2519
(throw 'pin%ex nil))))
2521
(m*l (mapcar #'pin%ex0 (cdr e))))
2523
(m+l (mapcar #'pin%ex0 (cdr e))))
2525
(throw 'pin%ex nil))))
2527
;; Test to see if exp is of the form p(x)*f(exp(x)). If so, set p* to
2528
;; be p(x) and set pe* to f(exp(x)).
2529
(defun p*pin%ex (nd*)
2530
(setq nd* ($factor nd*))
2531
(cond ((polyinx nd* var nil)
2532
(setq p* (cons nd* p*)) t)
2533
((catch 'pin%ex (pin%ex nd*))
2534
(setq pe* (cons nd* pe*)) t)
2536
(andmapcar #'p*pin%ex (cdr nd*)))))
2539
(cond ((findp p) nil)
2540
((setq nd* (bx**n p))
2541
(m^t var (car nd*)))
2542
((setq p (bx**n+a p))
2543
(m* (caddr p) (m^t var (cadr p))))))
2545
;; I think this is looking at f(exp(x)) and tries to find some
2546
;; rational function R and some number k such that f(exp(x)) =
2548
(defun funclogor%e (e)
2549
(prog (ans arg nvar r)
2550
(cond ((or (ratp e var)
2551
(involve e '(%sin %cos %tan))
2552
(not (setq arg (xor (and (setq arg (involve e '(%log)))
2556
ag (setq nvar (cond ((eq r '%log) `((%log) ,arg))
2557
(t (m^t '$%e arg))))
2558
(setq ans (maxima-substitute (m^t 'yx -1.) (m^t nvar -1.) (maxima-substitute 'yx nvar e)))
2559
(cond ((not (among var ans)) (return (list (subst var 'yx ans) nvar)))
2561
(setq arg (findsub arg)))
2564
;; Integration by parts.
2566
;; integrate(u(x)*diff(v(x),x),x,a,b)
2568
;; = u(x)*v(x)| - integrate(v(x)*diff(u(x),x))
2571
(defun dintbypart (u v a b)
2182
2572
;;;SINCE ONLY CALLED FROM DINTLOG TO get RID OF LOGS - IF LOG REMAINS, QUIT
2183
(let ((ad (antideriv v)))
2184
(COND ((or (NULL AD)
2185
(INVOLVE AD '(%LOG)))
2187
(t (let ((P1 (m* U AD))
2188
(P2 (m* AD (SDIFF U VAR))))
2189
(let ((P1-part1 (get-LIMit P1 VAR B '$MINUS))
2190
(p1-part2 (get-LIMIT P1 VAR A '$PLUS)))
2191
(cond ((or (null p1-part1)
2194
(t (let ((P2 (Let ((*DEF2* t))
2195
(DEFINT P2 VAR A B))))
2196
(COND (P2 (add* p1-part1
2573
(let ((ad (antideriv v)))
2574
(cond ((or (null ad)
2575
(involve ad '(%log)))
2577
(t (let ((p1 (m* u ad))
2578
(p2 (m* ad (sdiff u var))))
2579
(let ((p1-part1 (get-limit p1 var b '$minus))
2580
(p1-part2 (get-limit p1 var a '$plus)))
2581
(cond ((or (null p1-part1)
2584
(t (let ((p2 (let ((*def2* t))
2585
(defint p2 var a b))))
2586
(cond (p2 (add* p1-part1
2201
(DEFUN DINTEXP (EXP arg &aux ans)
2202
(let ((dintexp-recur t)) ;recursion stopper
2203
(COND ((and (sinintp exp var) ;To be moved higher in the code.
2591
;; integrate(f(exp(k*x)),x,a,b), where f(z) is rational.
2593
;; See Wang p. 96-97.
2595
;; If the limits are minf to inf, we use the substitution y=exp(k*x)
2596
;; to get integrate(f(y)/y,y,0,inf)/k. If the limits are 0 to inf,
2597
;; use the substitution s+1=exp(k*x) to get
2598
;; integrate(f(s+1)/(s+1),s,0,inf).
2599
(defun dintexp (exp arg &aux ans)
2600
(let ((*dintexp-recur* t)) ;recursion stopper
2601
(cond ((and (sinintp exp var) ;To be moved higher in the code.
2204
2602
(setq ans (antideriv exp))
2205
(setq ans (intsubs ans ll ul))))
2206
((setq ANS (FUNCLOGOR%E EXP))
2207
(COND ((AND (EQUAL LL 0.) (EQ UL '$INF))
2208
(SETQ EXP (SUBIN (m+t 1. arg) (CAR ANS)))
2209
(SETQ ANS (m+t -1. (CADR ANS))))
2210
(T (SETQ EXP (CAR ANS))
2211
(SETQ ANS (CADR ANS))))
2212
(INTCV ANS T NIL)))))
2603
(setq ans (intsubs ans ll ul)))
2604
;; If we can integrate it directly, do so and take the
2605
;; appropriate limits.
2607
((setq ans (funclogor%e exp))
2608
;; ans is the list (f(x) exp(k*x)).
2609
(cond ((and (equal ll 0.)
2611
;; Use the substitution s + 1 = exp(k*x). The
2612
;; integral becomes integrate(f(s+1)/(s+1),s,0,inf)
2613
(setq exp (subin (m+t 1. arg) (car ans)))
2614
(setq ans (m+t -1. (cadr ans))))
2616
;; Use the substitution y=exp(k*x) because the
2617
;; limits are minf to inf.
2618
(setq exp (car ans))
2619
(setq ans (cadr ans))))
2620
;; Apply the substitution and integrate it.
2621
(intcv ans t nil)))))
2214
(DEFUN DINTLOG (EXP ARG)
2215
(let ((dintlog-recur (f1+ dintlog-recur))) ;recursion stopper
2217
(COND ((AND (EQ UL '$INF)
2220
(EQUAL 1. ($RATSIMP (M// EXP (m* (m- (SUBIN (m^t VAR -1.)
2623
;; integrate(log(g(x))*f(x),x,0,inf)
2624
(defun dintlog (exp arg)
2625
(let ((*dintlog-recur* (f1+ *dintlog-recur*))) ;recursion stopper
2627
(cond ((and (eq ul '$inf)
2630
(equal 1. ($ratsimp (m// exp (m* (m- (subin (m^t var -1.)
2633
;; Make the substitution y=1/x. If the integrand has
2634
;; exactly the same form, the answer has to be 0.
2224
2636
((setq ans (antideriv exp))
2637
;; It's easy if we have the antiderivative.
2225
2638
(return (intsubs ans ll ul)))
2226
((SETQ ANS (LOGX1 EXP LL UL))
2228
(SETQ ANS (M// EXP `((%LOG) ,ARG)))
2229
(COND ((INVOLVE ANS '(%LOG)) (RETURN NIL))
2231
(EQUAL 0. (NO-ERR-SUB 0. ANS))
2232
(SETQ D (let ((*DEF2* t))
2233
(DEFINT (m* ANS (m^t VAR '*Z*))
2235
(RETURN (DERIVAT '*Z* 1. D 0.)))
2236
((SETQ ANS (DINTBYPART `((%LOG) ,ARG) ANS LL UL))
2239
(DEFUN DERIVAT (VAR N E PT) (SUBIN PT (APPLY '$DIFF (LIST E VAR N))))
2241
;;;MAYBPC RETURNS (COEF EXPO CONST)
2242
(DEFUN MAYBPC (E VAR)
2243
(COND (MTOINF* (THROW 'GGRM (LINPOWER0 E VAR)))
2245
(NULL (SETQ E (BX**N+A E)))) ;bx**n+a --> (a b n) or nil.
2246
NIL) ;with var being x.
2247
((AND (AMONG '$%I (CADDR E))
2248
(ZEROP1 ($realpart (CADDR E)))
2249
(SETQ ZN ($imagPART (CADDR E)))
2250
(EQ ($askSIGN (CADR E)) '$pos))
2251
(COND ((EQ ($askSIGN ZN) '$neg)
2255
(SETQ ZD (m^t '$%E (M// (mul* VAR '$%I '$%PI (m+t 1. ND*))
2256
(m*t 2. (CADR E)))))
2257
`(,ZN ,(CADR E) ,(CAR E)))
2258
((AND (OR (EQ (SETQ VAR ($askSIGN ($realPART (CADDR E)))) '$neg)
2260
(equal ($imagpart (CADR E)) 0)
2261
(RATGREATERP (CADR E) 0.))
2262
`(,(m- (CADDR E)) ,(CADR E) ,(CAR E)))))
2265
(PROG (C ZD ZN NN* DN* ND* DOSIMP $%EMODE)
2267
(COND (IND (SETQ E ($EXPAND E))
2268
(COND ((AND (mplusp E)
2269
(let ((*NODIVERG t))
2270
(SETQ E (CATCH 'Divergent
2275
(COND ((EQ E 'Divergent) NIL)
2276
(T (RETURN (sRATSIMP (CONS '(MPLUS) E)))))))))
2277
(SETQ E (RMCONST1 E))
2280
(COND ((SETQ E (GGR1 E VAR))
2281
(SETQ E (APPLY 'GAMMA1 E))
2282
(COND (ZD (SETQ $%EMODE T)
2284
(SETQ E (m* ZD E))))))
2285
(COND (E (RETURN (m* C E))))))
2290
(COND ((ATOM E) NIL)
2293
(COND ((SETQ E (MAYBPC (CADDR E) VAR)) (CONS 0. E))))
2296
(OR (AND (SETQ DN* (XTORTERM (CADR E) VAR))
2297
(RATGREATERP (SETQ ND* ($realPART DN*))
2299
(SETQ NN* (GGR1 (CADDR E) VAR)))
2300
(AND (SETQ DN* (XTORTERM (CADDR E) VAR))
2301
(RATGREATERP (SETQ ND* ($realPART DN*))
2303
(SETQ NN* (GGR1 (CADR E) VAR)))))
2308
;;; returns list of (a b n) or nil.
2313
(t (let ((a (NO-ERR-SUB 0. E)))
2315
(t (SETQ E (m+ E (m*t -1. A)))
2316
(COND ((SETQ E (BX**N E))
2321
;;;returns a list (n e) or nil.
2323
(AND (SETQ N (XEXPONGET E VAR))
2325
(SETQ E (let (($maxposex 1)
2327
($EXPAND (M// E (m^t VAR N)))))))
2330
(DEFUN XEXPONGET (E NN*)
2331
(COND ((ATOM E) (COND ((EQ E VAR) 1.)))
2335
(NOT (AMONG NN* (CADDR E))))
2337
(T (ORMAPC #'(LAMBDA (J)
2639
((setq ans (logx1 exp ll ul))
2641
;; Ok, the easy cases didn't work. We now try integration by
2642
;; parts. Set ANS to f(x).
2643
(setq ans (m// exp `((%log) ,arg)))
2644
(cond ((involve ans '(%log))
2645
;; Bad. f(x) contains a log term, so we give up.
2648
(equal 0. (no-err-sub 0. ans))
2649
(setq d (let ((*def2* t))
2650
(defint (m* ans (m^t var '*z*))
2652
;; The arg of the log function is the same as the
2653
;; integration variable. We can do something a little
2654
;; simpler than integration by parts. We have something
2655
;; like f(x)*log(x). Consider f(x)*x^z. If we
2656
;; differentiate this wrt to z, the integrand becomes
2657
;; f(x)*log(x)*x^z. When we evaluate this at z = 0, we
2658
;; get the desired integrand.
2660
;; So we need f(0) to be 0 at 0. If we can integrate
2661
;; f(x)*x^z, then we differentiate the result and
2662
;; evaluate it at z = 0.
2663
(return (derivat '*z* 1. d 0.)))
2664
((setq ans (dintbypart `((%log) ,arg) ans ll ul))
2665
;; Try integration by parts.
2668
;; Compute diff(e,var,n) at the point pt.
2669
(defun derivat (var n e pt)
2670
(subin pt (apply '$diff (list e var n))))
2674
;; MAYBPC returns (COEF EXPO CONST)
2676
;; This basically picks off b*x^n+a and returns the list
2677
;; (b n a). It may also set the global *zd*.
2678
(defun maybpc (e var)
2679
(declare (special *zd*))
2680
(cond (*mtoinf* (throw 'ggrm (linpower0 e var)))
2681
((and (not *mtoinf*)
2682
(null (setq e (bx**n+a e)))) ;bx**n+a --> (a n b) or nil.
2683
nil) ;with var being x.
2684
;; At this point, e is of the form (a n b)
2685
((and (among '$%i (caddr e))
2686
(zerop1 ($realpart (caddr e)))
2687
(setq zn ($imagpart (caddr e)))
2688
(eq ($asksign (cadr e)) '$pos))
2689
;; If we're here, b is complex, and n > 0. zn = imagpart(b).
2691
;; Set var to the same sign as zn.
2692
(cond ((eq ($asksign zn) '$neg)
2696
;; zd = exp(var*%i*%pi*(1+nd)/(2*n). (ZD is special!)
2697
(setq *zd* (m^t '$%e (m// (mul* var '$%i '$%pi (m+t 1. nd*))
2698
(m*t 2. (cadr e)))))
2700
`(,(caddr e) ,(cadr e) ,(car e)))
2701
((and (or (eq (setq var ($asksign ($realpart (caddr e)))) '$neg)
2703
(equal ($imagpart (cadr e)) 0)
2704
(ratgreaterp (cadr e) 0.))
2705
;; We're here if realpart(b) <= 0, and n >= 0. Then return -b, n, a.
2706
`(,(caddr e) ,(cadr e) ,(car e)))))
2708
;; Integrate x^m*exp(b*x^n+a), with realpart(m) > -1.
2710
;; See Wang, pp. 84-85.
2712
;; I believe the formula Wang gives is incorrect. The derivation is
2713
;; correct except for the last step.
2715
;; Let J = integrate(x^m*exp(%i*k*x^n),x,0,inf), with k < 0.
2717
;; Then J = exp(-%pi*%i*(m+1)/(2*n))*integrate(R^m*exp(k*R^n),R,0,inf)
2719
;; Wang seems to say this last integral is gamma(s/n/(-k)^s) where s =
2720
;; (m+1)/n. But that seems wrong. If we use the substitution x =
2721
;; (y/k)^(1/n), we end up with the result:
2723
;; integrate(y^((m+1)/n-1)*exp(-y),y,0,inf)/k^((m+1)/n)/n.
2725
;; or gamma((m+1)/n)/k^((m+1)/n)/n.
2729
(prog (c *zd* zn nn* dn* nd* dosimp $%emode)
2730
(declare (special *zd*))
2732
(cond (ind (setq e ($expand e))
2733
(cond ((and (mplusp e)
2734
(let ((*nodiverg t))
2735
(setq e (catch 'divergent
2740
(cond ((eq e 'divergent) nil)
2741
(t (return (sratsimp (cons '(mplus) e)))))))))
2742
(setq e (rmconst1 e))
2745
(cond ((setq e (ggr1 e var))
2746
;; e = (m b n a). I think we want to compute
2747
;; gamma((m+1)/n)/k^((m+1)/n)/n.
2749
;; FIXME: If n > m + 1, the integral converges. We need
2750
;; to check for this.
2752
(format t "e = ~A~%" e)
2753
(format t "asksign ~A = ~A~%"
2754
(sub (third e) (add ($realpart (first e)) 1))
2755
($asksign (sub (third e) (add ($realpart (first e)) 1)))))
2757
(setq e (apply #'gamma1 e))
2758
;; NOTE: *zd* (Ick!) is special and might be set by maybpc.
2760
;; FIXME: Why do we set %emode here? Shouldn't we just
2761
;; bind it? And why do we want it bound to T anyway?
2762
;; Shouldn't the user control that? The same goes for
2766
(setq e (m* *zd* e)))))
2767
(cond (e (return (m* c e))))))
2770
(prog (c *zd* zn nn* dn* nd* dosimp $%emode)
2771
(declare (special *zd*))
2773
(cond (ind (setq e ($expand e))
2774
(cond ((and (mplusp e)
2775
(let ((*nodiverg t))
2776
(setq e (catch 'divergent
2781
(cond ((eq e 'divergent) nil)
2782
(t (return (sratsimp (cons '(mplus) e)))))))))
2783
(setq e (rmconst1 e))
2786
(cond ((setq e (ggr1 e var))
2787
;; e = (m b n a). I think we want to compute
2788
;; gamma((m+1)/n)/k^((m+1)/n)/n.
2790
;; FIXME: If n > m + 1, the integral converges. We need
2791
;; to check for this.
2792
(destructuring-bind (m b n a)
2794
;; Check for convergence. If b is complex, we need n -
2795
;; m > 1. If b is real, we need b < 0.
2796
(when (and (zerop1 ($imagpart b))
2797
(not (eq ($asksign b) '$neg)))
2799
(when (and (not (zerop1 ($imagpart b)))
2800
(not (eq ($asksign (sub n (add m 1))) '$pos)))
2802
(setq e (gamma1 m (cond ((zerop1 ($imagpart b))
2803
;; If we're here, b must be negative.
2806
;; Complex b. Take the imaginary part
2807
`((mabs) ,($imagpart b))))
2809
;; NOTE: *zd* (Ick!) is special and might be set by maybpc.
2811
;; FIXME: Why do we set %emode here? Shouldn't we just
2812
;; bind it? And why do we want it bound to T anyway?
2813
;; Shouldn't the user control that? The same goes for
2817
(setq e (m* *zd* e))))))
2818
(cond (e (return (m* c e))))))
2821
;; Match x^m*exp(b*x^n+a). If it does, return (list m b n a).
2823
(cond ((atom e) nil)
2826
;; We're looking at something like exp(f(var)). See if it's
2827
;; of the form b*x^n+a, and return (list 0 b n a). (The 0 is
2828
;; so we can graft something onto it if needed.)
2829
(cond ((setq e (maybpc (caddr e) var))
2832
;; E should be the product of exactly 2 terms
2834
;; Check to see if one of the terms is of the form
2835
;; var^p. If so, make sure the realpart of p > -1. If
2836
;; so, check the other term has the right form via
2837
;; another call to ggr1.
2838
(or (and (setq dn* (xtorterm (cadr e) var))
2839
(ratgreaterp (setq nd* ($realpart dn*))
2841
(setq nn* (ggr1 (caddr e) var)))
2842
(and (setq dn* (xtorterm (caddr e) var))
2843
(ratgreaterp (setq nd* ($realpart dn*))
2845
(setq nn* (ggr1 (cadr e) var)))))
2846
;; Both terms have the right form and nn* contains the arg of
2847
;; the exponential term. Put dn* as the car of nn*. The
2848
;; result is something like (m b n a) when we have the
2849
;; expression x^m*exp(b*x^n+a).
2853
;; Match b*x^n+a. If a match is found, return the list (a n b).
2854
;; Otherwise, return NIL
2860
(t (let ((a (no-err-sub 0. e)))
2862
(t (setq e (m+ e (m*t -1. a)))
2863
(cond ((setq e (bx**n e))
2867
;; Match b*x^n. Return the list (n b) if found or NIL if not.
2870
(and (setq n (xexponget e var))
2872
(setq e (let (($maxposex 1)
2874
($expand (m// e (m^t var n)))))))
2877
(defun xexponget (e nn*)
2878
(cond ((atom e) (cond ((eq e var) 1.)))
2882
(not (among nn* (caddr e))))
2884
(t (ormapc #'(lambda (j)
2342
2889
;;; given (b*x^n+a)^m returns (m a n b)
2347
(INVOLVE E '(%LOG %SIN %COS %TAN))
2350
((mexptp E) (COND ((AMONG VAR (CADDR E)) NIL)
2351
((SETQ R (BX**N+A (CADR E)))
2352
(CONS (caddr e) R))))
2353
((SETQ R (BX**N+A E)) (CONS 1. R))
2894
(involve e '(%log %sin %cos %tan))
2897
((mexptp e) (cond ((among var (caddr e)) nil)
2898
((setq r (bx**n+a (cadr e)))
2899
(cons (caddr e) r))))
2900
((setq r (bx**n+a e)) (cons 1. r))
2355
2902
;;;Catches Unfactored forms.
2356
(SETQ M (M// (SDIFF E VAR) E))
2361
((AND (SETQ R (BX**N+A ($ratsimp r)))
2362
(NOT (AMONG VAR (SETQ M (M// M (m* (CADR R) (CADDR R)
2365
(SETQ E (M// (SUBIN 0. E) (m^t (CAR R) M))))
2368
(T (SETQ E (m^ E (M// 1. M)))
2369
(LIST M (m* E (CAR R)) (CADR R)
2370
(m* E (CADDR R))))))))
2903
(setq m (m// (sdiff e var) e))
2908
((and (setq r (bx**n+a ($ratsimp r)))
2909
(not (among var (setq m (m// m (m* (cadr r) (caddr r)
2912
(setq e (m// (subin 0. e) (m^t (car r) m))))
2915
(t (setq e (m^ e (m// 1. m)))
2916
(list m (m* e (car r)) (cadr r)
2917
(m* e (caddr r))))))))
2373
2920
;;;Is E = VAR raised to some power? If so return power or 0.
2375
(COND ((NOT (AMONG VAR E)) 0.)
2376
(T (XTORTERM E VAR))))
2922
(cond ((not (among var e)) 0.)
2923
(t (xtorterm e var))))
2378
(DEFUN XTORTERM (E VAR1)
2925
(defun xtorterm (e var1)
2379
2926
;;;Is E = VAR1 raised to some power? If so return power.
2380
(COND ((ALIKE1 E VAR1) 1.)
2383
(ALIKE1 (CADR E) VAR1)
2384
(NOT (AMONG VAR (CADDR E))))
2927
(cond ((alike1 e var1) 1.)
2930
(alike1 (cadr e) var1)
2931
(not (among var (caddr e))))
2388
(m^ (m* (m^ (CADDR L) '((RAT) 1. 2.))
2389
(m+ (CADR L) (m^ (m* (CAR L) (CADDR L))
2935
(m^ (m* (m^ (caddr l) '((rat) 1. 2.))
2936
(m+ (cadr l) (m^ (m* (car l) (caddr l))
2393
2940
(defun radbyterm (d l)
2394
2941
(do ((l l (cdr l))
2398
(let (((const . integrand) (rmconst1 (car l))))
2399
(setq ans (cons (m* const (dintrad0 integrand d))
2402
(DEFUN SQDTC (E IND)
2403
(PROG (A B C VARLIST)
2404
(SETQ VARLIST (LIST VAR))
2406
(SETQ E (CDADR (RATREP* E)))
2407
(SETQ C (PDIS (PTERM E 0.)))
2408
(SETQ B (m*t (m//t 1. 2.) (PDIS (PTERM E 1.))))
2409
(SETQ A (PDIS (PTERM E 2.)))
2410
(COND ((AND (EQ ($askSIGN (m+ B (m^ (m* A C)
2414
(NOT (EQ ($askSIGN A) '$neg))
2415
(EQ ($askSIGN C) '$pos))
2416
(AND (EQ ($askSIGN A) '$pos)
2417
(NOT (EQ ($askSIGN C) '$neg)))))
2418
(RETURN (LIST A B C))))))
2420
(DEFUN DIFAP1 (E PWR VAR M PT)
2421
(M// (mul* (COND ((eq (ask-integer M '$even) '$yes)
2425
(DERIVAT VAR M E PT))
2426
`((%GAMMA) ,(m+ PWR M))))
2428
(DEFUN SQRTINVOLVE (E)
2429
(COND ((ATOM E) NIL)
2432
(AND (MNUMP (CADDR E))
2433
(NOT (NUMBERP (CADDR E)))
2434
(EQUAL (CADDR (CADDR E)) 2.))
2435
(AMONG VAR (CADR E)))
2437
(T (ORMAPC #'SQRTINVOLVE (CDR E)))))
2439
(DEFUN BYDIF (R S D)
2441
(SETQ D (m+ (m*t '*Z* VAR) D))
2442
(COND ((OR (ZEROP1 (SETQ P (m+ S (m*t -1. R))))
2443
(AND (ZEROP1 (m+ 1. P))
2445
(DIFAP1 (DINTRAD0 B (m^ D '((RAT) 3. 2.)))
2446
'((RAT) 3. 2.) '*Z* R 0.))
2447
((EQ ($askSIGN P) '$pos)
2448
(DIFAP1 (DIFAP1 (DINTRAD0 1. (m^ (m+t 'Z** D)
2450
'((RAT) 3. 2.) '*Z* R 0.)
2451
'((RAT) 3. 2.) 'Z** P 0.)))))
2945
(destructuring-let (((const . integrand) (rmconst1 (car l))))
2946
(setq ans (cons (m* const (dintrad0 integrand d))
2949
(defun sqdtc (e ind)
2950
(prog (a b c varlist)
2951
(setq varlist (list var))
2953
(setq e (cdadr (ratrep* e)))
2954
(setq c (pdis (pterm e 0.)))
2955
(setq b (m*t (m//t 1. 2.) (pdis (pterm e 1.))))
2956
(setq a (pdis (pterm e 2.)))
2957
(cond ((and (eq ($asksign (m+ b (m^ (m* a c)
2961
(not (eq ($asksign a) '$neg))
2962
(eq ($asksign c) '$pos))
2963
(and (eq ($asksign a) '$pos)
2964
(not (eq ($asksign c) '$neg)))))
2965
(return (list a b c))))))
2967
(defun difap1 (e pwr var m pt)
2968
(m// (mul* (cond ((eq (ask-integer m '$even) '$yes)
2972
(derivat var m e pt))
2973
`((%gamma) ,(m+ pwr m))))
2975
(defun sqrtinvolve (e)
2976
(cond ((atom e) nil)
2979
(and (mnump (caddr e))
2980
(not (numberp (caddr e)))
2981
(equal (caddr (caddr e)) 2.))
2982
(among var (cadr e)))
2984
(t (ormapc #'sqrtinvolve (cdr e)))))
2986
(defun bydif (r s d)
2988
(setq d (m+ (m*t '*z* var) d))
2989
(cond ((or (zerop1 (setq p (m+ s (m*t -1. r))))
2990
(and (zerop1 (m+ 1. p))
2992
(difap1 (dintrad0 b (m^ d '((rat) 3. 2.)))
2993
'((rat) 3. 2.) '*z* r 0.))
2994
((eq ($asksign p) '$pos)
2995
(difap1 (difap1 (dintrad0 1. (m^ (m+t 'z** d)
2997
'((rat) 3. 2.) '*z* r 0.)
2998
'((rat) 3. 2.) 'z** p 0.)))))
2453
(DEFUN DINTRAD0 (N D)
2455
(COND ((AND (mexptp D)
2456
(EQUAL (DEG (CADR D)) 2.))
2457
(COND ((ALIKE1 (CADDR D) '((RAT) 3. 2.))
2458
(COND ((AND (EQUAL N 1.)
2459
(SETQ L (SQDTC (CADR D) T)))
2462
(SETQ L (SQDTC (CADR D) NIL)))
2463
(TBF (REVERSE L)))))
2464
((AND (SETQ R (FINDP N))
2465
(OR (EQ ($askSIGN (m+ -1. (m- R) (m*t 2.
2469
(SETQ S (m+ '((RAT) -3. 2.) (CADDR D)))
2470
(EQ ($askSIGN S) '$pos)
2471
(eq (ask-integer S '$integer) '$yes))
2472
(BYDIF R S (CADR D)))
2473
((POLYINX N VAR NIL)
2474
(RADBYTERM D (CDR N))))))))
3000
(defun dintrad0 (n d)
3002
(cond ((and (mexptp d)
3003
(equal (deg (cadr d)) 2.))
3004
(cond ((alike1 (caddr d) '((rat) 3. 2.))
3005
(cond ((and (equal n 1.)
3006
(setq l (sqdtc (cadr d) t)))
3009
(setq l (sqdtc (cadr d) nil)))
3010
(tbf (reverse l)))))
3011
((and (setq r (findp n))
3012
(or (eq ($asksign (m+ -1. (m- r) (m*t 2.
3016
(setq s (m+ '((rat) -3. 2.) (caddr d)))
3017
(eq ($asksign s) '$pos)
3018
(eq (ask-integer s '$integer) '$yes))
3019
(bydif r s (cadr d)))
3020
((polyinx n var nil)
3021
(radbyterm d (cdr n))))))))
2477
3024
;;;Looks at the IMAGINARY part of a log and puts it in the interval 0 2*%pi.
added
removed
src/defint.lisp
