Arbitrary precision addition

*To*: mathgroup at smc.vnet.net*Subject*: [mg25798] Arbitrary precision addition*From*: "Ersek, Ted R" <ErsekTR at navair.navy.mil>*Date*: Wed, 25 Oct 2000 03:53:52 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

Consider the lines below where (b) is given a value with 18 digits of precision, but the value we get for (a+b+ArcSin[3]) has a real part with 30 digits of precision. Instead the real part should have only 20 digits of precision. We also get an imaginary part with 48 digits of precision, but ArcSin[3] has an imaginary part with infinite precision. In[21]:= a=p-1/10; b=N[1/10,18]; a+b+ArcSin[3] {Precision at Re[%],Precision at Im[%]} Out[23]= 4.71238898038468985769396507492- 1.762747174039086050465218649959584618056320656523 I Out[24]= {30, 48} I asked a numerics developer from Wolfram Research about this and part of their response was, "This isn't really a bug, it is more of a historical design decision that numbers should be bounded so that the precision does not become negative." I still say this is a problem that can be corrected, so I recently coded an (Add) function which is designed to do the same thing as (Plus), but without the problem shown above. The result of using Add on the problem above is shown below. Notice the real part has the right amount of precision, and the imaginary part is exact. In[25]:= Add[a,b,ArcSin[3]] Precision[%] Out[25]= 4.712388980384689858 + I*Im[ArcSin[3]] Out[26]= 20 ----------------- I hope one of three things will come out of this. i) I will learn that my method can also give incorrect results. ii) I will learn that my method is very slow in some situations. iii) Wolfram Research will correct this problem in a future version. Maybe they will use a method similar to mine. The code for my (Add) function is provided below. ---------------------- Regards, Ted Ersek Down load Mathematica tips, tricks from http://www.verbeia.com/mathematica/tips/Tricks.html ----------------------- Attributes[Add]={Flat,Listable,NumericFunction,OneIdentity,Orderless}; AddReal[args__]:= Module[{pos,neg}, pos=Select[{args},Positive]; neg=Complement[{args},pos]; Plus@@{Plus@@pos,Plus@@neg} ] Add[args__]:= Module[{numbs,others,real1,cmplx,RealAccur,CmplxAccur,result}, numbs=Select[{args},NumericQ]; others=Complement[{args},numbs]; real1=Re[numbs]; cmplx=Im[numbs]; RealAccur=Min[Accuracy/@real1]+7; CmplxAccur=Min[Accuracy/@cmplx]+7; If[CmplxAccur===Infinity||MemberQ[cmplx, _?MachineNumberQ], cmplx=Plus@@cmplx, cmplx=AddReal@@(cmplx/.x_?ExactNumberQ:> Block[{Message},SetAccuracy[x,CmplxAccur]] ) ]; If[RealAccur===Infinity||MemberQ[real1, _?MachineNumberQ], real1=Plus@@real1, real1=AddReal@@(real1/.x_?ExactNumberQ:> Block[{Message},SetAccuracy[x,RealAccur]] ) ]; result=real1 + I*cmplx + Plus@@others; Expand[result]//. (I*t1_+I*t2_):>I(t1+t2) ] Protect[Add];