MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Text in graphics
  • Next by Date: Re: Summary: List element manipulation
  • Previous by thread: some Limits
  • Next by thread: Emmathfnt