On 6/14/2012 8:33 AM, Dmitry Shkirmanov wrote:
> Thanks, with this improvements maple still wins, but now it is only 3
> times faster than maxima. ( 66 seconds against 196 seconds for N= 6 in
> my loop). Is there possibility to improve this result more?
Yes, in addition to everything else suggested, you can just use
k1,k2,k3 instead of k[1], k[2], k[3],
and then you don't have to do some substitutions. This may save a
little time.
Perhaps more important is this possible change, which implements a rule
of thumb in
computer algebra "efficiency hacks". If you can encode your problem
into a simpler domain,
it will run faster. For example, removing sqrt, or removing exp or
changing a ratio of
polynomials to a polynomial (perhaps in 1/x)...
In this case, we may not need all the properties of the sqrt.... mostly
the derivative properties.
/* Introduce r(a,b,c) , which is sqrt(m^2+a^2+b^2+c^2) but we won't
tell Maxima that exactly, yet. */
rk: r(k1,k2,k3)$
gradef(r(a,b,c), a/rk, b/rk, c/rk)$
f :(1/2)*exp(-(1/4)*(r(k1,k2,k3)-m)^2/(sigma^2*rho[0,
0]))/(%pi^(3/2)*sigma*sqrt(rho[0, 0]*d)*r(k1,k2,k3))$
...
at the end,
facsum (subst (r = lambda([a,b,c], sqrt(m^2+a^2+b^2+c^2)), Z0[3]), m) ;
now this is probably not as fast as could be because at some points we
are missing
a simplification that r(a,b,c)^2 --> something without r.
This kind of hacking is perhaps unreasonable if you are just trying to
"do the same program"
in two different systems. To compare, you'd have to consider if Maple
could be made
faster, too.
There are indeed things that Maple does more efficiently, but a cleverly
devised benchmark
can turn it on its head and make it excruciatingly slow (compared to
Maxima).
Maple advocates tend to suppress conference papers that point out
problems in
Maple's model of computation.
RJF
The program looks something like this, now. (though it is probably buggy..)
kill(all)$
N:3$
showtime:true$
/* r(a,b,c) is sqrt(m^2+a^2+b^2+c^2) but we won't tell Maxima that */
rk: r(k1,k2,k3)$
gradef(r(a,b,c), a/rk, b/rk, c/rk)$
f :(1/2)*exp(-(1/4)*(r(k1,k2,k3)-m)^2/(sigma^2*rho[
0]))/(%pi^(3/2)*sigma*sqrt(rho[0]*d)*r(k1,k2,k3))$
f:rat(f)$
C[0] : 1/(2*%pi^(3/2)*sqrt(rho[ 0]*d)*m*sigma)$
Z[0] : f/C[0]$
for j: 0 thru N do(
d11 : diff(Z[j],k1,2),
d22 : subst([Sigma(1)=Sigma(2), Sigma(2)=Sigma(1), k1=k2, k2=k1], d11),
d33 : subst([Sigma(2)=Sigma(3), Sigma(3)=Sigma(2), k2=k3, k3=k2], d22),
Z[j+1]:ratsubst((m^2+k1^2+k2^2+k3^2),rk^2,
Sigma(1)*d11+Sigma(2)*d22+Sigma(3)*d33),
Z0[j+1] :subst([k1=0,k2=0,k3=0,r=lambda([a,b,c],m)], Z[j+1]),
Z0[j+1] :subst(Sigma= lambda([v],sigma^2*rho[v]), Z0[j+1]),
C[j+1] :Z0[j+1]/((j+1)!))$
/* now define r: */
facsum (subst (r = lambda([a,b,c], sqrt(m^2+a^2+b^2+c^2)), Z0[3]), m) ;