# The Poor Man's Integrator, a parallel integration heuristic
# Version 1.1 ---  May 10, 2005  (c) M.Bronstein and INRIA 2004-2005
pmint := proc(f,x)
	local ff, si, li, lin, lout, ld, q, d, l, vars, dx, ls, fint, lc;
	ff := eval(convert(f, tan));		# convert trigs to tan
	si := select(proc(d) diff(d,x) <> 0 end, indets(ff));
	si := select(proc(d) diff(d,x) <> 0 end, indets(map(diff, si, x))) union si;
	li := [op(si)];		# list of terms in integrand and its derivative
	lin := [seq(d=`tools/genglobal`(x), d=li)];	# substitution terms->indets
	lout := [seq(rhs(d)=lhs(d), d=lin)];		# substitution indets->terms
	ld := subs(lin, map(diff, li, x));	# derivatives of all the terms
	q := lcm(seq(denom(d), d=ld));		# denominator of the total derivation
	l := [seq(normal(q * d), d=ld)];	# q * derivatives of all the terms
	vars := map(lhs, lout);
	dx := totalDerivation(vars, l);		# vector field dx = q * d/dx
	ls := [seq(getSpecial(d, lin), d=li)];	# list of known Darboux for dx
	fint := subs(lout, pmIntegrate(subs(lin, ff), dx, q, vars, ls));
	lc := select(proc(d) convert(d,string)[1]="_"; end, indets(fint, name));
	subs({seq(d = 0, d=lc minus si)}, fint);
end;

getSpecial := proc(f, l) local p;		# return known Darboux polys
	p := op(0,f);
	if p = `tan` then [1+subs(l,f)^2, false];
	elif p = `tanh` then [1 + subs(l,f), false], [1 - subs(l,f), false];
	elif p = `LambertW` then [subs(l,f), true];
	else NULL; fi;
end;

totalDerivation := proc(lv, ld)
	proc(f) local fp, i;
		fp := 0; for i to nops(ld) do fp := fp + ld[i] * diff(f, lv[i]); od;
		fp;
	end;
end;

pmIntegrate := proc(f, d, q, vars)
	local ls, splq, s, ff, df, spl, cden, dg, monomials, cand, lunk, sol, i;
	if nargs = 5 then ls := args[5]; else ls := []; fi;
	splq := splitFactor(q, d);
	s := splq[1]; for i to nops(ls) do if ls[i][2] then s := s*ls[i][1]; fi; od;
	ff := normal(f); df := denom(ff); spl := splitFactor(df, d);
	cden := s * spl[1] * deflation(spl[2], d);
	dg := 1 + degree(s) +  max(degree(numer(ff)), degree(denom(ff)));
	monomials := [op(enumerateMonoms(vars, dg))];
	cand := add('_A'[i] * monomials[i], i = 1..nops(monomials)) / cden;
	lunk := { seq('_A'[i], i = 1..nops(monomials)) };
	sol:= tryIntegral(f, d, q, vars, cand, lunk, spl[1], spl[2], splq[1], ls, 0);
	if sol[1] then sol := tryIntegral(f, d, q, vars, cand, lunk,
						spl[1], spl[2], splq[1], ls, I); fi;
	if sol[1] then Int(f); else sol[2]; fi;
end;

tryIntegral := proc(f, d, q, vars, cand, lunk, l1, l2, l3, ls, K)
	local candlog, p, candidate, i, sol;
	candlog := [op({ myfactors(l1, K), myfactors(l2, K), myfactors(l3, K) }
			union { seq(p[1], p=ls) })];
	candidate := cand + add('_B'[i] * log(candlog[i]), i = 1..nops(candlog));
	sol := solve({coeffs(numer(normal(f - d(candidate)/q)), {op(vars)})},
				lunk union { seq('_B'[i], i = 1..nops(candlog)) });
	[evalb(sol=NULL), subs(sol,candidate)];
end;

myfactors := proc(p, K) local l, fact;
	if K = 0 then l := factors(p); else l := factors(p, K); fi;
	seq(fact[1], fact=l[2]);
end;

enumerateMonoms := proc(vars, d) local n, x, i, v, s, w;
	n := nops(vars);
	if n = 0 then {1}; else
		x := vars[n];
		v := [seq(vars[i], i = 1..n-1)];
		s := enumerateMonoms(v, d);
		for i to d do s := s union {seq(x^i*w,w=enumerateMonoms(v,d-i))}; od;
		s;
	fi;
end;

splitFactor := proc(p, d) local si, x, c, q, spl, s, splh;
	si := select(proc(z) d(z) <> 0 end, indets(p,name));
	if si = {} then RETURN([1,p]) fi;
	x := si[1];
	c := content(p, x, 'q');
	spl := splitFactor(c, d);
	s := normal(gcd(q, d(q)) / gcd(q, diff(q, x)));
	if degree(s) = 0 then RETURN([spl[1], q * spl[2]]); fi;
	splh := splitFactor(normal(q / s), d);
	[spl[1] * splh[1] * s, spl[2] * splh[2]];
end;

deflation := proc(p, d) local si, x, c, q;
	si := select(proc(z) d(z) <> 0 end, indets(p,name));
	if si = {} then RETURN(p) fi;
	x := si[1];
	c := content(p, x, 'q');
	deflation(c, d) * gcd(q, diff(q, x));
end;