Overview

On March 14th, convenience stores and department stores are crowded with sweets sales, and I think many people are impressed that spring is approaching. People in the world are also paying attention to Pi Day, and I think that it will have a considerable economic effect. So, today I will write a program to calculate pi.

It is a simple story.

Gauss-Legendre's algorithm

Wikipedia's [Gauss-Legendre's Algorithm](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%EF%BC%9D%E3%83 % AB% E3% 82% B8% E3% 83% A3% E3% 83% B3% E3% 83% 89% E3% 83% AB% E3% 81% AE% E3% 82% A2% E3% 83% AB According to% E3% 82% B4% E3% 83% AA% E3% 82% BA% E3% 83% A0), it is relatively easy to implement.

I thought it was a great genius to think like this.

python implementation

However, I don't know how to calculate pi, so I looked at the sources of various people and it was completed. The feature is that Decimal is used to represent a large number of digits.

# -*- coding:utf-8 -*-
import math
import time
from decimal import *

import numpy as np


#1002 characters
PI_1000 = "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"


def get_pi(prec=1000, verbose=False):
    '''
Returns the circumference ratio of prec with the number of digits after the decimal point as a character string
    '''
    prec = prec+1+1 #accuracy"3."Since there is a part of, increase the accuracy by one, and if you do not have another digit, it may shift due to the roundness, so add one more
    N = 2*math.ceil(math.log2(prec)) #Number of repetitions According to wikipedia, log2(prec)It seems that the degree is enough, so just in case, double it

    getcontext().prec = prec 

    a = Decimal(1.0)
    b = Decimal(1.0) / Decimal(2.0).sqrt()
    t = Decimal(0.25)
    p = Decimal(1.0)


    start_time = time.clock() 

    #Start N trials
    for _ in range(1, N):
        a_next = (a + b)/Decimal(2)
        b_next = (a*b).sqrt()
        t_next = t - p * (a - a_next)**2
        p_next = Decimal(2)*p

        a = a_next
        b = b_next
        t = t_next
        p = p_next

    #Calculate pi
    pi = (a + b)**2 / (Decimal(4)*t)
    #Calculate execution time
    elapsed = time.clock()  - start_time

    #Display the execution result
    if verbose:
        print("accuracy:",prec)
        print("Pi:", pi)
        print("Execution time:%f" % (elapsed))

    return str(pi)[0:prec]

I ran it as below and confirmed that about 1000 digits were correct. (Some people calculate about 200 million digits, but this is fine.)

pi = get_pi(prec=1000, verbose=True)
print(pi)
#print(PI_1000)
#print(pi == PI_1000)

Accuracy: 1002
Pi: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410
270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925
409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244
065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079
227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313
783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019890
Execution time:0.014581

Statistics

[Introduction to Statistics](https://www.amazon.co.jp/ Introduction to Statistics-Basic Statistics-Statistics Department, Faculty of Liberal Arts, University of Tokyo / dp / 4130420658), Appearance of numbers appearing in pi There was a problem to check the rate, so I tried to solve it.

After looking at the graph, we also perform a chi-square test to see if the numbers are constant.

When 100 numbers are taken from the 123rd digit

The 123rd digit is appropriate.

start = 123
pi_123 = pi[2+start:2+start+100]

freq = np.zeros(10, dtype=np.int)
for a in pi_123:
    freq[int(a)] += 1
print("freq(100)=",freq)

#Ask using chisquare
print(chisquare(freq, 10*np.ones(10,dtype=np.int8)))

plt.xlim(-0.5, 9.5)
plt.bar(np.array([0,1,2,3,4,5,6,7,8,9]), freq, align='center')

Intuitively, looking at the graph, it feels unpleasant, but since the p-value is about 0.37, there is no contradiction even if the appearance rate is almost constant.

freq(100)= [ 9 12 11  7 15 13  7  4 12 10]
Power_divergenceResult(statistic=9.8000000000000007, pvalue=0.3669177991127523)

At the time of 1000 digits

freq = np.zeros(10, dtype=np.int)
for a in pi:
    if a != ".":
        freq[int(a)] += 1
print("freq(1001)=",freq)
print(chisquare(freq, 100.1*np.ones(10,dtype=np.int8)))

plt.xlim(-0.5, 9.5)
plt.bar(np.array([0,1,2,3,4,5,6,7,8,9]), freq, align='center')

When about 1000 pieces of data are collected, there is no contradiction even if the appearance rate is considered to be constant regardless of the number. The p-value is 0.85, which is good.

freq(1001)= [ 93 116 103 103  93  97  94  95 101 106]
Power_divergenceResult(statistic=4.7842157842157853, pvalue=0.85269838594638103)

To be worried

Actually, I wanted to know the statistical information of the sequence of numbers, but that is the next. I have to study transcendental numbers.

March 14th is Pi Day. The story of calculating pi with python