Just for fun

[Dru Lee Parsec]

Just for fun I coded the Leibniz series for calculating pi in Java and compared it to the static variable Math.PI that's contained in the java language. I think it's a pretty sweet algorithm. Quite simple in it's beauty.


public class doPi {

public static void main(String[] args) {
double pi = 0;
double subpi = 0;
double iteration = 10000000;
for (double n = 0; n < iteration; n++ ) {
subpi += ( Math.pow(-1,n)/((2*n) +1) );
}
pi = subpi * 4;
System.out.println(pi + " is pi");
System.out.println(Math.PI + " is Math.PI" );
}
}


Anyone want to code it up in another language?

[GnuVince]

I made a program a while back that infinitely calculates pi. It's written in O'Caml. Only glitch: the period after the leading 3 does not appear, but who cares?


open Big_int

let _ =
let k = ref (big_int_of_int 2) in
let a = ref (big_int_of_int 4) in
let b = ref unit_big_int in
let a1 = ref (big_int_of_int 12) in
let b1 = ref (big_int_of_int 4) in

(* Moo *)
let a2 = ref zero_big_int in
let b2 = ref zero_big_int in

while true do
let p = (square_big_int !k) in
let q = (add_big_int (mult_big_int (big_int_of_int 2) !k) unit_big_int) in
k := (add_big_int !k unit_big_int);


a2:= !a;
b2:= !b;
a := !a1;
b := !b1;
a1:= (add_big_int (mult_big_int p !a2) (mult_big_int q !a1));
b1:= (add_big_int (mult_big_int p !b2) (mult_big_int q !b1));

let d = ref (div_big_int !a !b) in
let d1 = ref (div_big_int !a1 !b1) in

while (eq_big_int !d !d1) do
print_string (string_of_big_int !d);
flush stdout;
a := mult_big_int (mod_big_int !a !b) (big_int_of_int 10);
a1:= mult_big_int (mod_big_int !a1 !b1) (big_int_of_int 10);
d := div_big_int !a !b;
d1:= div_big_int !a1 !b1;
done;
done



To compile: ocamlopt nums.cmxa pi.ml -o pi

[GnuVince]

[vince@vincent: ~/prog/ocaml/contest/pi]% java doPi
3.1415925535897915 is pi
3.141592653589793 is Math.PI
[vince@vincent: ~/prog/ocaml/contest/pi]% ./pi


A flaw?

[Dru Lee Parsec]

A flaw?

I don't think so. I think it's more of a case of not enough iterations. In fact if the iterations are scaled down to only 50,000 then you get this output:


3.1415726535897814 is pi
3.141592653589793 is Math.PI
......^


So I think it's simply a by product of the series. The more iterations, the more accurate it gets. But this particular series needs LOTS of iterations to become relativly accurate.

I mean, the original code looped 10 million times. The actual series definition has n going from 0 to infinity.

[GnuVince]

Ah ok. You might want to try my program :) I'm going to write yours in O'Caml now.

[Dru Lee Parsec]

GnuVince: Do you have a web link to the algorithm you used for your program? Honestly, I can't make heads or tails out of what O'Caml is doing from looking at the code.

[GnuVince]

I don't have a link, but I have a Ruby version (which is quite readable). Here it is:


k, a, b, a1, b1 = 2, 4, 1, 12, 4

while TRUE
# Next approximation
p, q, k = k*k, 2*k+1, k+1
a, b, a1, b1 = a1, b1, p*a+q*a1, p*b+q*b1
# Print common digits
d = a / b
d1 = a1 / b1
while d == d1
print d
$stdout.flush
a, a1 = 10*(a%b), 10*(a1%b1)
d, d1 = a/b, a1/b1
end
end



You must use big int's by the way.

[Dru Lee Parsec]

In Ruby is the percent symbol % a modulus operator?
i.e.
1%3 = 1
2%3 = 2
3%3 = 0
4 % 3 = 1
5 % 3 = 2
6 % 3 = 0
and so on . . .

[GnuVince]

Yes.

[Strike]

Any alternating series (such as the one in Dru's first post) has a maximum possible error from the final answer of the next term in the series.

So, to get precision to 0.000 000 001 (1/100 000 000), you need to do n iterations where 2n+1 ?= 100 000 000, or n = 50 000 000. It doesn't converge very quickly. There are otherways of doing it, though, as GnuVince showed.

[GnuVince]

Strike: Yes, my technique is faster, but uses integers (well, big integers) instead of floats

[Dru Lee Parsec]

Well, my code looks right, but the first time I go into the while (d == d1) part of the loop I have these values


p is 4.0 q is 5.0 a is 12.0 a1 is 108.0
b is 4.0 b1 is 36.0


d = a/b which is 3 and d1 = a1/b1 which is also 3 so we print out a 3. However, the modulus (the remainder) of a%b and a1%b1 is zero. So this code:


a = 10*(a%b);
a1 = 10 * (a1%b1);
d = a/b;
d1 = a1/b1;


sets everything to zero because a and a1 become zero. So from that point it loops forever.

Of course all of this assumes that in Ruby this line:

d, d1 = a/b, a1/b1


is equivalent to this


d = a/b
d1 = a1/b1


Do you see where the problem is?

[GnuVince]

Using a comma in Ruby makes a simultaneous assignement. IIRC, I had this problem with my O'Caml program. So I added a2 and b2, and using them fixed my problem. I can probably try to use a2 and b2 in the Ruby program (which could help you a lot)

[file13]

GnuVince: Do you have a web link to the algorithm you used for your program? Honestly, I can't make heads or tails out of what O'Caml is doing from looking at the code.

i keep telling him to COMMENT! ;)

yeah i remember this beast....still gives me headaches looking at it.... :grr:

[GnuVince]

I *did* comment!


(* Moo *)


Doesn't that explain everything?

[file13]

:D

why yes! it's all clear now!

;)