[Axiom-developer] 20070921.01.tpd.patch

[daly]

The call 

 laplace(log(z),z,w)

causes an infinite loop between lfintegrate and monomialIntegrate.

Waldek modified the code to fail and return unevaluated.

A more reasonable return value, not supported by Axiom, would be:

  -log(w)-gamma
  -------------
       w


======================================================================
diff --git a/changelog b/changelog
index 48493c1..d085c48 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,6 @@
+20070920 tpd src/input/Makefile add bug101.input regression test
+20070920 tpd src/input/bug101.input test laplace(log(z),z,w) bug 101
+20070920 wxh src/algebra/laplace.spad fix laplace(log(z),z,w) bug 101
 20070916 tpd src/input/Makefile add bug103.input regression test
 20070916 tpd src/input/bug103.input test solve(z=z,z) bug fix
 20070916 tpd src/algebra/polycat.spad solve(z=z,z) bug fix
diff --git a/src/algebra/laplace.spad.pamphlet 
b/src/algebra/laplace.spad.pamphlet
index 6f81b2e..76cc837 100644
--- a/src/algebra/laplace.spad.pamphlet
+++ b/src/algebra/laplace.spad.pamphlet
@@ -161,41 +161,68 @@ LaplaceTransform(R, F): Exports == Implementation where
 
     lapkernel(f, t, tt, ss) ==
       (k := retractIfCan(f)@Union(K, "failed")) case "failed" => "failed"
-      empty?(arg := argument(k::K)) or not empty? rest arg => "failed"
+      empty?(arg := argument(k::K)) => "failed"
+      is?(op := operator k, "%diff"::SE) =>
+        not( #arg = 3) => "failed"
+        not(is?(arg.3, t)) => "failed"
+        fint := eval(arg.1, arg.2, tt)
+        s := name operator (kernels(ss).1)
+        ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0)
+      not (empty?(rest arg)) => "failed"
       member?(t, variables(a := first(arg) / tt)) => "failed"
       is?(op := operator k, "Si"::SE) => atan(a / ss) / ss
       is?(op, "Ci"::SE) => log((ss**2 + a**2) / a**2) / (2 * ss)
       is?(op, "Ei"::SE) => log((ss + a) / a) / ss
+      -- digamma (or Gamma) needs SpecialFunctionCategory
+      -- which we do not have here
+      -- is?(op, "log"::SE) => (digamma(1) - log(a) - log(ss)) / ss
       "failed"
 
+    -- Below we try to apply one of the texbook rules for computing
+    -- Laplace transforms, either reducing problem to simpler cases
+    -- or using one of known base cases
     locallaplace(f, t, tt, s, ss) ==
       zero? f => 0
 --      one? f  => inv ss
       (f = 1)  => inv ss
+
+      -- laplace(f(t)/t,t,s) 
+      --              = integrate(laplace(f(t),t,v), v = s..%plusInfinity)
       (x := tdenom(f, tt)) case F =>
         g := locallaplace(x::F, t, tt, vv := new()$SE, vvv := vv::F)
         (x := intlaplace(f, ss, g, vv, vvv)) case F => x::F
         oplap(f, tt, ss)
+
+      -- Use linearity
       (u := mkPlus f) case List(F) =>
         +/[locallaplace(g, t, tt, s, ss) for g in u::List(F)]
       (rec := splitConstant(f, t)).const ^= 1 =>
         rec.const * locallaplace(rec.nconst, t, tt, s, ss)
+
+      -- laplace(t^n*f(t),t,s) = (-1)^n*D(laplace(f(t),t,s), s, n))
       (v := atn(f, t)) case Record(coef:F, deg:PI) =>
         vv := v::Record(coef:F, deg:PI)
         is?(la := locallaplace(vv.coef, t, tt, s, ss), oplap) => oplap(f,tt,ss)
         (-1$Integer)**(vv.deg) * differentiate(la, s, vv.deg)
+
+      -- Complex shift rule
       (w := aexp(f, t)) case Record(coef:F, coef1:F, coef0:F) =>
         ww := w::Record(coef:F, coef1:F, coef0:F)
         exp(ww.coef0) * locallaplace(ww.coef,t,tt,s,ss - ww.coef1)
+
+      -- Try base cases
       (x := lapkernel(f, t, tt, ss)) case F => x::F
-      -- last chance option: try to use the fact that
-      --    laplace(f(t),t,s) = s laplace(g(t),t,s) - g(0)  where dg/dt = f(t)
-      elem?(int := lfintegrate(f, t)) and (rint := retractIfCan int) case F =>
-           fint := rint :: F
-           -- to avoid infinite loops, we don't call laplace recursively
-           -- if the integral has no new logs and f is an algebraic function
-           empty?(logpart int) and algebraic?(f, t) => oplap(fint, tt, ss)
-           ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0)
+
+--    -- The following does not seem to help computing transforms, but
+--    -- quite frequently leads to loops, so I (wh) disabled it for now
+--    -- last chance option: try to use the fact that
+--    --    laplace(f(t),t,s) = s laplace(g(t),t,s) - g(0)  where dg/dt = f(t)
+--    elem?(int := lfintegrate(f, t)) and (rint := retractIfCan int) case F =>
+--        fint := rint :: F
+--        -- to avoid infinite loops, we don't call laplace recursively
+--        -- if the integral has no new logs and f is an algebraic function
+--        empty?(logpart int) and algebraic?(f, t) => oplap(fint, tt, ss)
+--        ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0)
       oplap(f, tt, ss)
 
     setProperty(oplap,SPECIALDIFF,dvlap@((List F,SE)->F) pretend None)
@@ -228,6 +255,7 @@ InverseLaplaceTransform(R, F): Exports == Implementation 
where
       ++ Laplace transform of \spad{f(s)}
       ++ using t as the new variable or "failed" if unable to find
       ++ a closed form.
+      ++ Handles only rational \spad{f(s)}.
 
   Implementation ==> add
 
@@ -246,8 +274,12 @@ InverseLaplaceTransform(R, F): Exports == Implementation 
where
     ilt(expr,var,t) ==
       expr = 0 => 0
       r := univariate(expr,kernel(var))
+
+      -- Check that r is a rational function such that degree of
+      -- the numarator is lower that degree of denominator
       not(numer(r) quo denom(r) = 0) => "failed"
       not( freeOf?(numer r,var) and freeOf?(denom r,var)) => "failed"
+
       ilt1(r,t::F)
 
     hintpac := TranscendentalHermiteIntegration(F, UP)
diff --git a/src/input/Makefile.pamphlet b/src/input/Makefile.pamphlet
index b286aa6..b2d5fe9 100644
--- a/src/input/Makefile.pamphlet
+++ b/src/input/Makefile.pamphlet
@@ -287,7 +287,8 @@ REGRES= algaggr.regress algbrbf.regress  algfacob.regress 
alist.regress  \
     arrows.regress    assign.regress   atansqrt.regress \
     asec.regress      bags.regress      bbtree.regress \
     binary.regress    bop.regress      bstree.regress   bouquet.regress \
-    bug100.regress    bug103.regress   bug10069.regress \
+    bug100.regress    bug101.regress \
+    bug103.regress    bug10069.regress \
     bugs.regress      bug10312.regress bug6357.regress  bug9057.regress \
     calcprob.regress  \
     calculus2.regress calculus.regress cardinal.regress card.regress \
@@ -502,7 +503,8 @@ FILES= ${OUT}/algaggr.input  ${OUT}/algbrbf.input    
${OUT}/algfacob.input \
        ${OUT}/bags.input     ${OUT}/bbtree.input     ${OUT}/bern.input \
        ${OUT}/bernpoly.input ${OUT}/binary.input     ${OUT}/bop.input \
        ${OUT}/bouquet.input  ${OUT}/bstree.input     ${OUT}/bug6357.input \
-       ${OUT}/bug9057.input  ${OUT}/bug100.input     ${OUT}/bug103.input \
+       ${OUT}/bug9057.input  ${OUT}/bug100.input     ${OUT}/bug101.input \
+       ${OUT}/bug103.input \
        ${OUT}/bug10069.input ${OUT}/bug10312.input \
        ${OUT}/calcprob.input ${OUT}/calculus.input \
        ${OUT}/cardinal.input ${OUT}/card.input       ${OUT}/carten.input \
@@ -678,7 +680,8 @@ DOCFILES= \
   ${DOC}/bernpoly.input.dvi    ${DOC}/binary.input.dvi     \
   ${DOC}/bop.input.dvi         ${DOC}/bouquet.input.dvi    \
   ${DOC}/bstree.input.dvi      ${DOC}/bug10069.input.dvi   \
-  ${DOC}/bug100.input.dvi      ${DOC}/bug103.input.dvi \
+  ${DOC}/bug100.input.dvi      ${DOC}/bug101.input.dvi     \
+  ${DOC}/bug103.input.dvi \
   ${DOC}/bug10312.input.dvi    ${DOC}/bug6357.input.dvi    \
   ${DOC}/bug9057.input.dvi     ${DOC}/bugs.input.dvi       \
   ${DOC}/c02aff.input.dvi      ${DOC}/c02agf.input.dvi     \
diff --git a/src/input/bug101.input.pamphlet b/src/input/bug101.input.pamphlet
new file mode 100644
index 0000000..b82caf0
--- /dev/null
+++ b/src/input/bug101.input.pamphlet
@@ -0,0 +1,59 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input bug101.input}
+\author{Timothy Daly}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+@
+The call 
+\begin{verbatim}
+laplace(log(z),z,w)
+\end{verbatim}
+causes an infinite loop between lfintegrate and monomialIntegrate.
+
+Waldek modified the code to fail and return unevaluated.
+
+A more reasonable return value, not supported by Axiom, would be:
+\begin{verbatim}
+  -log(w)-gamma
+  -------------
+       w
+\end{verbatim}
+
+<<*>>=
+)spool bug101.output
+)set message test on
+)set message auto off
+)clear all
+ 
+--S 1 of 2
+laplace(log(z),z,w)
+--R
+--R   (1)  laplace(log(z),z,w)
+--R                                                     Type: Expression 
Integer
+--E 1
+
+--S 2 of 2
+laplace(log(z),w,z)
+--R
+--R        log(z)
+--R   (2)  ------
+--R           z
+--R                                                     Type: Expression 
Integer
+--E 2
+@
+<<*>>=
+)spool 
+)lisp (bye)
+ 
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}