[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.
More information about the erlang-questions
mailing list