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

jm <>
Thu Apr 3 06:30:20 CEST 2008

```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.

```