Neither the Lisp nor the Maxima code given collects terms completely.
Note below that the first three terms should be collected.
ex:1/(x^3-1)$
pfex: partfrac(ex,x) => 1/(3*(x-1))-(x+2)/(3*(x^2+x+1))
pfmult(pfex,pfex+1,x) =>
-(x+2)/(3*(x^2+x+1))
+2*(x+1)/(9*(x^2+x+1))
+1/(9*(x^2+x+1))
+(x+1)/(3*(x^2+x+1)^2)
+1/(9*(x-1))
+1/(9*(x-1)^2)
pfmult_by_distrib(pfex,pfex+1,x) =>
(2*x+2)/(9*(x^2+x+1))
-(x+2)/(3*(x^2+x+1))
+1/(9*(x^2+x+1))
+(x+1)/(3*(x^2+x+1)^2)
+1/(9*(x-1))
+1/(9*(x-1)^2)
partfrac(pfex*(pfex+1),x) =>
-(x+3)/(9*(x^2+x+1))
+(x+1)/(3*(x^2+x+1)^2)
+1/(9*(x-1))
+1/(9*(x-1)^2)
Also, the Maxima version doesn't work for p=q, because distrib doesn't
expand powers.
pfmult_by_distrib(1+1/x,1+1/x,x) => (1/x+1)^2
-s