On 2012-06-06, David Ronis <David.Ronis at McGill.CA> wrote:
> domain:complex;
> f(z):=sqrt(z+Gamma)/((z+D*k^2)*(zp-z+D*k^2));
> integrate(f(%i*omega)*exp(%i*omega*t), omega, minf, inf);
>
> When I run this I get:
>
> *** - Program stack overflow. RESET
Looks like PTIMES%E attempts to split off part of the integrand and
calls itself to handle the reduced part. But the attempted reduction
just yields the same expression over and over. POLFACTORS is supposed to
compute the reduced integrand, I guess.
Here are the first few calls to PTIMES%E and POLFACTORS. You can see how
the infinite regress gets rolling.
1 Enter ?\?polfactors
[k^4*sqrt(%i*omega+Gamma)*%e^(%i*omega*t)*D^2+k^2*sqrt(%i*omega+Gamma)*%e^(%i*omega*t)*zp*D
-%i*omega*sqrt(%i*omega+Gamma)*%e^(%i*omega*t)*zp
+omega^2*sqrt(%i*omega+Gamma)*%e^(%i*omega*t)]
1 Exit ?\?polfactors false
1 Enter ?\?ptimes%e [omega^2*sqrt(%i*omega+Gamma)*%e^(%i*omega*t),4]
1 Enter ?\?polfactors [omega^2*sqrt(%i*omega+Gamma)*%e^(%i*omega*t)]
1 Exit ?\?polfactors
(?mexpt,?simp,?ratsimp)(sqrt(%i*omega+Gamma)*%e^(%i*omega*t))
2 Enter ?\?ptimes%e [sqrt(%i*omega+Gamma)*%e^(%i*omega*t),4]
1 Enter ?\?polfactors [sqrt(%i*omega+Gamma)*%e^(%i*omega*t)]
1 Exit ?\?polfactors (1,sqrt(%i*omega+Gamma)*%e^(%i*omega*t))
3 Enter ?\?ptimes%e [sqrt(%i*omega+Gamma)*%e^(%i*omega*t),4]
1 Enter ?\?polfactors [sqrt(%i*omega+Gamma)*%e^(%i*omega*t)]
1 Exit ?\?polfactors (1,sqrt(%i*omega+Gamma)*%e^(%i*omega*t))
4 Enter ?\?ptimes%e [sqrt(%i*omega+Gamma)*%e^(%i*omega*t),4]
1 Enter ?\?polfactors [sqrt(%i*omega+Gamma)*%e^(%i*omega*t)]
I'm guessing that the leading 1 in the POLFACTORS return value indicates
it didn't find a nontrivial factor. If so the patch shown below could
fix the bug.
After applying the patch I get
'integrate(sqrt(%i*omega+Gamma)*%e^(%i*omega*t)/((k^2*D+%i*omega)*(k^2*D+zp-%i*omega)),omega,minf,inf)
as the output, which isn't too satisfying but better than an error.
What does anybody think about the patch? Maybe POLFACTORS isn't trying
hard enough or something like that?
best
Robert Dodier
PS. Here is PTIMES%E from src/defint.lisp, in extenso.
(defun ptimes%e (term n)
(cond ((and (mexptp term)
(eq (cadr term) '$%e)
(polyinx (caddr term) var nil)
(eq ($sign (m+ (deg ($realpart (caddr term))) -1))
'$neg)
(eq ($sign (m+ (deg (setq nn* ($imagpart (caddr term))))
-2.))
'$neg))
(cond ((eq ($asksign (ratdisrep (ratcoef nn* var))) '$pos)
(setq *updn t))
(t (setq *updn nil)))
term)
((and (mtimesp term)
(setq nn* (polfactors term))
(or (null (car nn*))
(eq ($sign (m+ n (m- (deg (car nn*)))))
'$pos))
(ptimes%e (cadr nn*) n)
term))
(t (throw 'ptimes%e nil))))
PPS. Here's a patch.
diff --git a/src/defint.lisp b/src/defint.lisp
index 1459caa..3a8f9f7 100644
--- a/src/defint.lisp
+++ b/src/defint.lisp
@@ -1111,6 +1111,7 @@ in the interval of integration.")
term)
((and (mtimesp term)
(setq nn* (polfactors term))
+ (not (and (consp nn*) (equal 1 (car nn*)))) ;; no
nontrivial factor found
(or (null (car nn*))
(eq ($sign (m+ n (m- (deg (car nn*)))))
'$pos))