/* taken out ....... use the lisp version!

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))$   */
    
 matchdeclare(ss, 
   partition_expression("*",lambda([u],freeof(x,u)),[],cons,"[",'ANSss))$
 defrule(ddr1,ss,ANSss)$

(
/* %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)]),
 intable[polygamma]         :lambda([u,v], if freeof(x,u) then
                              [polygamma(u-1,v),diff(v,x)]),


/* 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)]),
/* some extra pieces courtesy of Barton Willis */
intable[abs]       : lambda([u], [u * abs(u) / 2, diff(u, %voi)]),
intable[signum]    : lambda([u], [abs(u), diff(u, %voi)]),
intable[unit_step] : lambda([u], [(u + abs(u))/2, diff(u, %voi)])

)
$
(
 is_sum(x):= is (not(atom(x)) and (inpart(x,0)="+")),
 inargs(z):=substinpart("[",z,0), /* utility like args, but better on a/b */

 matchdeclare(ss, partition_expression("*",lambda([u],freeof(%voi,u)),
                                          [],cons,"[",'ANSss)),
 defrule(ddr1,ss,ANSss), /*this forces the call to partition_expression */

intfudu(exp,%voi):= /*integrate exp=f(u)*du wrt %voi*/
 if is_sum(exp) then map (lambda([r],intfudu(r,%voi)), exp) else 
  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:intable[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 diff(k,%voi) # 0 and freeof(%voi,therest:ratsimp(exp/k/diff(k,%voi)))
               then (result:(therest*k^2/2),
                            return())),
  return(if (result=false) then 'int(exp,%voi) else result)))
 
$
/* ratsimp is optional there.. */
"infudu loaded";

done;