[erlang-questions] how: big float aithmetics (was PI calculation)

[jm]

Zvi wrote:
> Jeff,
> 
> thank you for the solution. But I think, that mixing big integers with
> floats limits precision by the float.

I should have removed all the floating point operations. It is possible 
though. The real problem is that factor of 0.5. If you write out the 
equation, starting with,


      1             4                 (I - 0.5)
pi = - \Sigma ---------   where X = -----------
      N         1 + X^2                   N



rearranging you get

                        1
pi = 4 \Sigma ---------------------
                      1
               N + N( - (I - 0.5))^2
                      N

multiple though by 1/N

pi                     1
--  =  4 \Sigma -------------------
N                N^2 + (I - 0.5)^2


Multiply through top and by the square of a scaling factor SF


                      SF^2
pi = 4N \Sigma --------------------------------
               SF^2 * N^2 + (I * SF - SF / 2)^2


Notice all the constants. It now becomes

calc_pi(fixedpoint_serial_decr, N) ->
	DecimalPlaces = 10,
         %% fudge doesn't work properly should also create an Epsilon
	SFig = round(math:pow(10,DecimalPlaces)),
	SF = SFig * SFig,
       	SF2 = SF * SF,
	SF2N2 =	SF2 * N	* N,
	HalfSF = SF div	2,
	Numerator = SF2*SF2,  %% another fudge
         calc_pi1(N, N, Numerator, SF, SF2N2, HalfSF, 0).

calc_pi1(0, N, _SF2, _SF2N2, _HalfSF, Sum) -> 4*Sum*N;
calc_pi1(I, N,  SF2,  SF2N2,  HalfSF, Sum) ->
   %% io:format("I ~p N ~p Sum ~p~n", [I, N, Sum]),
   X = I * SF - HalfSF,
   V = SF2 div (SF2N2 + X*X)
   calc_pi1(I-1, N, SF2, SF2N2, HalfSF, Sum + V).


modifing calc_pi:test/0 and running gives
26> calc_pi:test().
[{math,{4127,3.14159}},
  {list_creation,{116696,3.14159}},
  {lc,{1117037,3.14159}},
  {lc_no_pow,{813246,3.14159}},
  {map,{736112,3.14159}},
  {mapfoldl,{601372,3.14159}},
  {serial_decr,{272162,3.14159}},
  {serial,{276182,3.14159}},
  {parallel,{275360,3.14159}},
 
{fixedpoint_serial_decr,{3828596,31415926535898765717959767164129964000000}}]
27> 3828596 - 272162.
3556434

Yes a bit slower. However...

Using the value of Pi given at 
http://www.research.att.com/~njas/sequences/A000796

this seems to equate to 14 places of accuracy,

 > Pi41 = 31415926535897932384626433832795028841971.
 > (31415926535898765717959767164129964000000 - Pi41)/math:pow(10,41).
8.33333e-15


And the going gets tough,

 > {_,P8} = timer:tc(calc_pi,calc_pi,[fixedpoint_serial_decr, 100000000]).
{752066430,31415926535897932467959747166535600000000}
 > (Pi41 - P8)/math:pow(10,41).
-8.33333e-19


each digit takes more and more processing.

> The more generic approach is to represent "big floats" using tuples, like
> {BigInt, Exponent}.

Won't be as fast as specifically optimised code but I can understand the 
reasoning.

> I think it works, but I have trouble to implement division using integer
> "div" operator:


Perhaps, try the older versions of the Numerical Recipes books ( 
http://www.nr.com/oldverswitcher.html ). I seem to remember something in 
one of them about arbitrary precision arithmetic.

Also try http://citeseer.ist.psu.edu/smith91fortran.html or Knuth 1981


That's about all I can think of at the moment.


Jeff.