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];