/* Expression Partitioning tool.
   Included is an application implementing stage 1 of Moses' SIN
   symbolic indefinite integration, intfudu.
  
   Author: Richard Fateman, January, 2009. Suggestions by B. Willis.
   Documentation http://www.cs.berkeley.edu/~fateman/papers/partition.pdf */

  /*  Moving expression arg to end makes matchdeclare simpler */

partition_expression(operator,pred,init,combiner,action,res,E):=
    block([yes:init,no:init],
   if not atom(E) and inpart(E,0)=operator then
       (map(lambda([r], if apply(pred,[r])=false
              then  no:apply(combiner,[r,no])
              else yes:apply(combiner,[r,yes])),
            inargs(E) ),
       res :: apply(action, [yes,no]), /* result stored as requested */
       true))$

 inargs(z):=substinpart("[",z,0)$ /* utility like args, but avoids / or - */

 (
/* %voi = variable of integration, a global variable wrt these entries */

intablefun(op):= if atom(op) then intable[op] else false,

intable[otherwise]:=false, /*backstop for undefined kernels */
 intable[log] : lambda([u],  [-u+u*log(u),diff(u,%voi)]),
 intable[sin] : lambda([u],  [-cos(u),diff(u,%voi)]),
 intable[cos] : lambda([u],  [sin(u),diff(u,%voi)]),
 intable[tan] : lambda([u],  [log(sec(u)),diff(u,%voi)]),
 intable[sec] : lambda([u],  [log(tan(u)+sec(u)),diff(u,%voi)]),
 intable[csc] : lambda([u],  [log(tan(u/2)),diff(u,%voi)]),
 intable[cot] : lambda([u],  [log(sin(u)),diff(u,%voi)]),
 intable[atan]: lambda([u],  [-(log(u^2+1)-2*u*atan(u))/2,diff(u,%voi)]),
 intable[acos]: lambda([u],  [-sqrt(1-u^2)+u*acos(u),diff(u,%voi)]),
 intable[asin]: lambda([u],  [sqrt(1-u^2)+u*asin(u),diff(u,%voi)]),
 intable[sinh]: lambda([u],  [cosh(u),diff(u,%voi)]),
 intable[cosh]: lambda([u],  [sinh(u),diff(u,%voi)]),
 intable[tanh]: lambda([u],  [log(cosh(u)),diff(u,%voi)]),
 intable[sech]: lambda([u],  [atan(sinh(u)),diff(u,%voi)]),
 intable[csch]: lambda([u],  [log(tanh(u/2)),diff(u,%voi)]),
 intable[coth]: lambda([u],  [log(sinh(u)),diff(u,%voi)]),
 intable[asinh]: lambda([u],  [-sqrt(1-u^2)+u*asinh(u),diff(u,%voi)]),
 intable[acosh]: lambda([u],  [-sqrt(u^2-1)+u*acosh(u),diff(u,%voi)]),
 intable[atanh]: lambda([u],  [log(1-u^2)/2+u*atanh(u),diff(u,%voi)]),
 intable[acsch]:lambda([u], [u*asinh(u)/abs(u)+u*acsch(u),diff(u,%voi)]),
 intable[asech]:lambda([u], [u*asech(u)-atan(sqrt(1/u^2-1)),diff(u,%voi)]),
 intable[acoth]:lambda([u], [log(u^2-1)/2+u*acoth(u),diff(u,%voi)]),

 intable["^"] : lambda([u,v], if freeof(%voi,u) then [u^v/log(u),diff(v,%voi)]
                                else if freeof(%voi,v)
                                then if (v#-1) then [u^(v+1)/(v+1), diff(u,%voi)]
                                else [log(u),diff(u,%voi)]),

/* Rule for integration of 'diff(f(u),u) wrt u for "unknown" f */

/* Need a fancier rule for integration of 'diff(f(u^2),u) wrt u.
  Maxima does not have notation for diff wrt first argument,
  though see pdiff share file */

 intable[nounify(diff)]:lambda([[u]], 
  if length(u)=3 and u[2]=%voi then [diff(u[1],%voi,u[3]-1),1]),

/* We can add rules for more operations like this */

 intable[polygamma] :lambda([u,v], if freeof(%voi,u) then [polygamma(u-1,v),diff(v,%voi)]),
 intable[Ci]: lambda([u], [u*Ci(u)-sin(u),diff(u,%voi)]),
 intable[Si]: lambda([u], [u*Si(u)-cos(u),diff(u,%voi)]),
 gradef(Ci(w), cos(w)/w),
 gradef(Si(w), sin(w)/w),
 intable[nounify(bessel_j)] :lambda([u,v], if(u=1) then [-bessel_j(0,v),diff(v,%voi)]),
 intable[nounify(bessel_i)] :lambda([u,v], if(u=1) then [ bessel_i(0,v),diff(v,%voi)]),
 intable[nounify(bessel_k)] :lambda([u,v], if(u=1) then [-bessel_k(0,v),diff(v,%voi)]),

/* Here are Airy functions.   */
intable[nounify(airy_ai)]:lambda([u],
[-(u*(-3*gamma(1/3)*gamma(5/3)*hgfred([1/3], [2/3, 4/3], 
      u^3/9) + 3^(1/3)*u*gamma(2/3)^2*hgfred([2/3], 
      [4/3, 5/3], u^3/9)))/(9*3^(2/3)*gamma(2/3)*gamma(4/3)*gamma(5/3)),diff(u,%voi)]),


/* Here are Legendre polynomials P[n](x).  Presumably if n were an explicit
   integer this symbolism would be removed, so we assume it is of symbolic
   order n. */

intable[legendre_p] :
  lambda([n,u], if freeof(n,u) then 
          [(legendre_p(n+1,u)-legendre_(n-1,u))/(2*n+1),
            diff(u,%voi)])
)$

/*etc */

/* Here's the main integration program */

( matchdeclare(ss, partition_expression("*",lambda([u],freeof(%voi,u)),[],cons,"[",'ANSss)),
 defrule(ddr1,ss,ANSss))$

/*intfudu stands for integrate f(u)*du. This is equivalent to Moses' SIN, stage 1 */

 intfudu(exp,%voi):= /*integrate exp=f(u)*du wrt %voi*/ 
  block([lists,consts,factors, thefuns, therest, thelist, int, df, result:false],
   if freeof(%voi,exp) then return (exp*%voi),
   lists:ddr1(exp), /*partition expression into factors*/
   if lists=false then lists:[[1],[exp]],
   factors:second(lists),
   for k in factors do 
    if not atom(k) and 
      (thefuns:intablefun(inpart(k,0)))#false and
      (thelist:apply(thefuns,inargs(k)))#false
     then( [int,df]:thelist,
           if freeof(%voi,therest:ratsimp(exp/k/df))
               then (result:(therest*int),
                            return())),
   /*if nothing of form f(u)du worked, then try matching u^1*u'-> u^2/2 */
  if result=false then  
      for k in factors do 
          if freeof(%voi,therest:ratsimp(exp/k/diff(k,%voi)))
               then (result:(therest*k^2/2),
                            return()),
  return(result))$