/* 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") $