How julia calculates pi?

Offtopic
1

I think I have found a simple formula to calculate pi but my formula is much slower than BigFloat (pi).
How julia calculates pi?

I also ask you: how do I know if my formula already exists?

2

For BigFloat(pi), Julia uses the MPFR C library’s mpfr_const_pi function. More details about Julia’s π-handling are in this 2017 post on the Julia Blog.

6 Likes
3

One way would be to scour the literature to see what algorithms are out there, you could probably start with the Wikipedia page on pi.

Another way, if you don’t have the training to do an exhaustive literature search, is to write your method up as clean as you can and post it on something like r/maths. Caution, though, if you lack that training the likelihood that your method is genuinely novel is vanishingly small. The way you’re formulating it here (“I came up with an algorithm, is it new?”) is more likely to get a supportive response than the more common “Here’s the new fastest way to do it and I’ve named it for myself”.

2 Likes
4

Can I post the formula here?

2 Likes
5

sure. why not?

6

@gustaphe said to post on reddit math.
It is not a question of the julia language.
(I signed up on reddit a few minutes ago and got blocked)

7

Post some Julia code that implements it and the math that the Julia code is based on. I don’t think anyone will complain. Maybe if you can, edit the topic to “offtopic”. As long as you’re Julia adjacent and honestly asking questions not pushy or confrontational I don’t think anyone is going to mind. Some people might even give you tips about how to improve the speed.

7 Likes
8

First one(π):

a = BigFloat(1) 
# you can also use other numbers 
# and you will converge to multiples of pi (or -pi, or 0)
a = a + sin(a) # 1.84147098480789650665250232163
a = a + sin(a) # 2.8050617093497299136094750235101
a = a + sin(a) # 3.1352763328997160003524035699574
a = a + sin(a) # 3.141592611590653149601161318361
a = a + sin(a) # and so on..... -> π

Second one (π/2 = 1.5707963267948966192313216916397):

a = BigFloat(1) 
# you can also use other numbers 
a = a + cos(a) # 1.5403023058681397174009366074441
a = a + cos(a) # 1.5707916010242611729993559629573
a = a + cos(a) # 1.5707963267948966016412881420516
a = a + cos(a) # and so on..... -> π/2

Third one (π):

a = BigFloat(3)  
# you can also use other numbers
a = a - tan(a) # 3.14254654307427780529563541053391
a = a - tan(a) # 3.14159265330047681544988577171991
a = a - tan(a) # and so on..... -> pi

Do these formulas exist?

3 Likes
9

These are what you’d call “fixed point iterations”, whether an iteration converges to a fixed point can be analyzed in terms of various properties of “iterative maps” which are studied in dynamical systems theory.

4 Likes
10

Thank you.
So these formulas already exist?

11

They are one example of a well known family of processes, yes.

6 Likes
12

Sind discussion here:

2 Likes
13

For those who are curious how it is implemented, you can read the source here: src/const_pi.c ¡ master ¡ mpfr / mpfr ¡ GitLab

3 Likes
14

Sorry if I came across as rude. I just meant that you’ll get better answers if you choose a specific forum, it could have happened that nobody here knew the answer, while someone on r/maths could have pointed you in the right direction. You were in luck this time :slight_smile:

2 Likes
15

Interesting! It’s the Brent-Salamin formula.

3 Likes
16

This gets my hopes high that there will be an application for the geothmetic meandian one day

3 Likes
17

Indeed! Adapting a previous reference to this concept on this forum:

F(x) = [(x[1]+x[2])/2, √(x[1]x[2]), (x[1]-x[2])/2, x[1]^2 - x[2]^2]

Fⁿ(x,n) = ∘(fill(F,n)...)(x)

Q(n) = [Fⁿ([1,1/√BigFloat(2)],i) for i=1:n]

Pi(n) = (P = Q(n); 4*(P[n][1])^2/(1 - sum(2^(j)*P[j][4] for j=2:n)))

This function has impressive convergence:

julia> Pi(2) - π
-0.001013403067624990151312114303679572423729161862157485331457898962233648251246769

julia> Pi(3) - π
-7.376250956313298951296807109882732176029503026415405107486381639065704112971525e-09

julia> Pi(4) - π
-1.831306084776389098159438945396588851876138902814342684477531468625531008649221e-19

julia> Pi(5) - π
-5.472109145689941832748533178964178489138195886835923610924029609286031340771535e-41

julia> Pi(6) - π
3.109020679834000065139086670608143845601760131197061298608533249260693049138474e-76

julia> Pi(7) - π
3.109020679834000065139086670608143845601760131197061298608533249260693049138474e-76

(It exhausts the BigFloat precision after 6 iterations.)

4 Likes
18

Sorry for spamming the thread; here’s Brent’s actual algorithm in seven lines:

## Brent's method for calculation pi:
# Ref.: https://maths-people.anu.edu.au/~brent/pd/rpb028.pdf
# See page 10.
F(x) = [(x[1]+x[2])/2, √(x[1]x[2]), x[3] - x[4]*((x[1]-x[2])^2)/4, 2*x[4]]
compose(F,n) = ∘(fill(F,n)...)
function Pi(n)
    A, B, T, X = compose(F,n)([1,inv(sqrt(big(2))),1/4,1])
    return (A+B)^2/(4T)
end
# Compare:
Diff_to_π = [Pi(k) - big(π) for k in 1:5]

I’m done now!

8 Likes
19

Hi,

I once looked into various ways of computing pi, I implemented some of those in Pi.jl. This was a while a ago, so the code might not be up to date with julia 1.0.

2 Likes