/* this is all the maxima-language code in the partition paper, all in one file */
kill(all);
/* Our partitioning tool. Note that the expression being partitioned
is the last argument. This 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)$
/* inargs is a utility like args() except args(a/b) is [a,1/b]. we want [a,b] */
(matchdeclare (pp, partition_expression("+",constantp,0,"+",g,'ANSpp)),
tellsimp(foo(pp),ANSpp),
declare([a,b,c],constant),
[foo(a+b+c+x+y+3),foo(w),foo(a), foo(a*b+x),foo(a*b)]);
(matchdeclare(qq, partition_expression("*",constantp,1,"*",g,'ANSqq)),
tellsimp(bar(qq),ANSqq),
[bar(a*b*3*x*y),bar(w),bar(a+b)]);
(matchdeclare(ss,
partition_expression("*",constantp,[],cons,"[",'ANSss)),
tellsimp(bar2(ss),ANSss),
tellsimp(bar3(ss), The_Partitions_Are(ANSss)),
[bar2(3*%pi*xx), bar2(w),bar2(%pi+x+3), bar3(x*3)]);
(matchdeclare(odd,oddp, nov, freeof(%v)),
defmatch(OneOddPowerp,nov*%v^odd), /*%v is global */
matchdeclare(tt,
partition_expression("+",OneOddPowerp,{},adjoin,"[",'ANStt)),
tellsimp(separatepowers(tt),ANStt),
[block([%v:x],separatepowers(3*x+4*x^5+7*x^10)),
block([%v:y],separatepowers(3*y+4*y^5+7*x^11)),
block([%v:a],separatepowers(expand((a+x)^5)))]);
/* %voi is ``variable of integration'', a global variable */
/* shorthand for du/dv used in most table entries */
dudv(u):=diff(u,%voi)$
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),dudv(u)])$
intable[sin] : lambda([u], [-cos(u),dudv(u)])$
intable[cos] : lambda([u], [sin(u),dudv(u)])$
intable[tan] : lambda([u], [log(sec(u)),dudv(u)])$
intable[sec] : lambda([u], [log(tan(u)+sec(u)),dudv(u)])$
intable[csc] : lambda([u], [log(tan(u/2)),dudv(u)])$
intable[cot] : lambda([u], [log(sin(u)),dudv(u)])$
intable[atan] : lambda([u], [-(log(u^2+1)-2*u*atan(u))/2,dudv(u)])$
intable[acos] : lambda([u], [-sqrt(1-u^2)+u*acos(u),dudv(u)])$
intable[asin] : lambda([u], [sqrt(1-u^2)+u*asin(u),dudv(u)])$
intable[sinh] : lambda([u], [cosh(u),dudv(u)])$
intable[cosh] : lambda([u], [sinh(u),dudv(u)])$
intable[tanh] : lambda([u], [log(cosh(u)),dudv(u)])$
intable[sech] : lambda([u], [atan(sinh(u)),dudv(u)])$
intable[csch] : lambda([u], [log(tanh(u/2)),dudv(u)])$
intable[coth] : lambda([u], [log(sinh(u)),dudv(u)])$
intable[asinh]: lambda([u], [-sqrt(1-u^2)+u*asinh(u),dudv(u)])$
intable[acosh]: lambda([u], [-sqrt(u^2-1)+u*acosh(u),dudv(u)])$
intable[atanh]: lambda([u], [log(1-u^2)/2+u*atanh(u),dudv(u)])$
intable[acsch]: lambda([u], [u*asinh(u)/abs(u)+u*acsch(u),dudv(u)])$
intable[asech]: lambda([u], [u*asech(u)-atan(sqrt(1/u^2-1)),dudv(u)])$
intable[acoth]: lambda([u], [log(u^2-1)/2+u*acoth(u),dudv(u)])$
intable["^"] : lambda([u,v], if freeof(%voi,u) then [u^v/log(u),dudv(v)]
else if freeof(%voi,v)
then if (v#-1) then [u^(v+1)/(v+1), dudv(u)]
else [log(u),dudv(u)])$
/* Rule for integration of 'diff(f(u),u) with respect to u for "unknown" f */
/* Need a fancier rule for integration of 'diff(f(u^2),u) with respect to u.
Maxima does not have built-in notation for diff with respect to 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),dudv(u)])$
intable[Si]: lambda([u], [u*Si(u)-cos(u),dudv(u)])$
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),dudv(v)])$
intable[nounify(bessel_i)] :lambda([u,v], if(u=1) then [ bessel_i(0,v),dudv(v)])$
intable[nounify(bessel_k)] :lambda([u,v], if(u=1) then [-bessel_k(0,v),dudv(v)])$
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)),dudv(u)])$
/* 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_p(n-1,u))/(2*n+1), dudv(u)])$
/* some extra pieces courtesy of Barton Willis */
intable[abs] : lambda([u], [u * abs(u) / 2, dudv(u)])$
intable[signum] : lambda([u], [abs(u), dudv(u)])$
intable[unit_step] : lambda([u], [(u + abs(u))/2, dudv(u)])$
is_sum(x):= is (not(atom(x)) and inpart(x,0)="+")$
matchdeclare(ss, partition_expression("*",lambda([u],freeof(%voi,u)),
[],cons,"[",'ANSss))$
defrule(ddr1,ss,ANSss)$
/* Finally, here's the main integration program */
intfudu(E,%voi):= /*integrate E=f(u)*du with respect to %voi*/
if is_sum(E) then map (lambda([r],intfudu(r,%voi)), E) else
block([lists,consts,factors, thefuns, therest, thelist, int, df,
result:false],
if freeof(%voi,E) then return (E*%voi),
lists:ddr1(E), /*partition expression into factors*/
if lists=false then lists:[[1],[E]],
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(E/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(E/k/diff(k,%voi)))
then (result:(therest*k^2/2),
return())),
return(if (result=false) then 'int(E,%voi) else result))$
print ("intfudu loaded") $