rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

abs() or cabs() are much slower than just sqrt(re^2+im^2) #3831

Closed rtoy closed 3 months ago

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-08 18:15:41 Created by s_i_m on 2010-01-07 12:57:47 Original: https://sourceforge.net/p/maxima/bugs/1868


Doing some calculations I noticed that there is something wrong with imlementation of abs() and cabs(). Basically, in the code below I can calculate and plot sqrt(re^2 + im^2), but I cannot do it with abs(), nor with cabs(): it never finishes and Maxima does not react to the keyboard interrupt either.

I run this version of maxima (on Debian sid):

Maxima 5.20.1 http://maxima.sourceforge.net using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (a.k.a. GCL)

========================================================

eps : [1,2];

mu : [1,1];

n : sqrt(eps*mu);

vc : [0,0.5];

epst : eps*(1-vc^2)/(1-n^2*vc^2);

mut : mu*(1-vc^2)/(1-n^2*vc^2);

a : vc*(n^2-1)/(1-n^2*vc^2);

kxe : -%i*sqrt(-k0^2*eps*mut+eps/epst*(kz+k0*a)^2+ky^2);

kxh : -%i*sqrt(-k0^2*mu*epst+mu/mut*(kz+k0*a)^2+ky^2);

betae : k0*kxe*epst/(k0^2*epst*mut-(kz+k0*a)^2);

betah : k0*kxh*mut/(k0^2*epst*mut-(kz+k0*a)^2);

alpha : -ky*(kz+k0*a)/(k0^2*epst*mut-(kz+k0*a)^2);

a1 : part(alpha,1); a2 : part(alpha,2);

be1 : part(betae,1); be2 : part(betae,2);

bh1 : part(betah,1); bh2 : part(betah,2);

denom : (a1-a2)^2+(be1+be2)*(bh1+bh2);

Ree : -((a1-a2)^2-(be1-be2)*(bh1+bh2))/denom;

plot3d(sqrt(realpart(subst(k0=%i,Ree))^2+imagpart(subst(k0=%i,Ree))^2),[ky,-3,3],[kz,-3,3]);

plot3d(cabs(subst(k0=%i,Ree)),[ky,-3,3],[kz,-3,3]);

plot3d(abs(subst(k0=%i,Ree)),[ky,-3,3],[kz,-3,3]);

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-08 18:15:42 Created by crategus on 2010-01-10 23:35:11 Original: https://sourceforge.net/p/maxima/bugs/1868/#93b7


Fixed in revision 1.28 in rpart.lisp. An extra expand in the routine absarg has been removed. This expand has caused an enormous growing of the expression of this example.

With this change the time to get the rectform and the absolute value are comparable. The plot for the absolute value is faster than the plot with the sqrt function.

Example with sbcl on a slow notebook:

(%i7) plot3d(sqrt(realpart(test)^2 +imagpart(test)^2),[ky,-3,3],[kz,-3,3]); Evaluation took 119.8990 seconds (122.6900 elapsed) using 3259.698 MB. (%o7) false

(%i8) plot3d(abs(test),[ky,-3,3],[kz,-3,3]); Evaluation took 54.8320 seconds (56.3300 elapsed) using 628.454 MB. (%o8) false

Closing this bug report as fixed. Dieter Kaiser