PyPrintPi på en Raspberry Pi (12 / 23 steg)

Steg 12: Kodning Gauss formel

I vår föregående steg funnit vi formeln:

Π/4=12*ARCTAN(1/18) + 8*arctan(1/57) - 5*arctan(1/239)

Vi kan enkelt förvandla detta till en funktion i python om vi antar att vi har definierat en arctan funktion:

def gauss_pi_method(precision):

Return 4*(12*arctan(18) + 8*arctan(57) - 5*arctan(239))

Nu måste vi skapa en arctan funktion. Den snabbaste metoden för att genomföra skulle vara att använda formeln vi såg i steg 9:

ARCTAN(x) = x - (x³/3) + (x⁵/5) - (x⁷/7) + (x⁹/9) - (x¹¹/11)...

men ser så vi ska bara räkna ut arctans med siffror i form av 1 / x det är vettigare att omdefiniera formeln som:

ARCTAN(1/x) = (1 / x)-(1/3 x ³) + (1/5x⁵) - (1/7x⁷) + (1/9x⁹) - (1/11 x ¹¹)...

Detta ger oss följande funktion i Python:

def arctan(x): #defins en funktion för att definiera pi med formeln vi härrör
noggrannhet = Decimal(Decimal(1)/Decimal(10**(getcontext().prec-3)))
resultat = Decimal(Decimal(1)/Decimal(x)) #result är en variabel används för att lagra våra nuvarande värdet på pi som ett decimaltal, 1 och x förvandlas till decimaler att undvika avrundningsfel
tecken =-1 #sign håller spår om att addera eller subtrahera andelen
denominatorval = 3 #
fract = 0 #fraction håller värdet av den del vi vill addera eller subtrahera till resultatet, det används senare för att se om funktionen kan sluta.

samtidigt sant:
fract = Decimal(Decimal(1)/(denominatorval*(Decimal(x)**denominatorval)))
resultatet += sign*(fract)
om fract < noggrannhet:
Break
tecknet * = -1
denominatorval += 2

returnera resultat

Denna funktion är lite mer avancerad än vår tidigare funktioner eftersom det inte kräver att användaren kan ange hur många iterationer för att köra. I stället man tittar på värdet av getcontext () .prec och söker efter det minsta värde som Python programmet kan skilja från noll (till exempel om getcontext () .prec = 2 så är det lägsta värde som Python kan skilja från 0 0,1).

Sedan jämförs värdet av den förra mandatperioden sekvensen: om det är mindre än det minsta värdet python kan skilja från 0 då det är meningslöst att fortsätta så programmet slutar.

Om vi sammanställt de två funktionerna får vi följande program (gauss_pi_method.py):

från decimal import *

def gauss_pi_method():
Return 4*(12*arctan(18) + 8*arctan(57) - 5*arctan(239))

def arctan(x): #defins en funktion för att definiera pi med formeln vi härrör

noggrannhet = Decimal(Decimal(1)/Decimal(10**(getcontext().prec-3)))
resultat = Decimal(Decimal(1)/Decimal(x)) #result är en variabel används för att lagra våra nuvarande värdet på pi som ett decimaltal, 1 och x förvandlas till decimaler att undvika avrundningsfel
tecken =-1 #sign håller spår om att addera eller subtrahera andelen
denominatorval = 3 #
fract = 0 #fraction håller värdet av den del vi vill addera eller subtrahera till resultatet, det används senare för att se om funktionen kan sluta.

samtidigt sant:
fract = Decimal(Decimal(1)/(denominatorval*(Decimal(x)**denominatorval)))
resultatet += sign*(fract)
om fract < noggrannhet:

Break

tecknet * = -1

denominatorval += 2

returnera resultat

medan 1:
#Asks användaren hur många decimaler som de vill ha det svara för att ges till
decimaler = int (input ("Ange hur många decimaler vill du ha svaret ges till:")) + 1

getcontext () .prec = decimaler + 4 #sets precision med 4 fler decimaler än bad om.

PiEst = gauss_pi_method() #computes pi

getcontext () .prec = decimaler #this kodblock rundor pi antal decimaler som bad om
PiEst = + PiEst

#next line skriver ut resultaten
skriva ut ("Pi till" + str(decimals-1) + "decimaler är:" + str(PiEst))

Om du kör detta för att hitta 10 000 siffror av π det tar under en minut. På min dator tog det 17 sekunder. Det är mycket bättre än den tidigare bästa programmet vi hade, den som är baserad på metoden polygon som tog 43 sekunder på min dator för att beräkna 100 siffror av π.

Men vi kan göra mycket bättre med samma funktion för att beräkna π om vi kunde beräkna arctan(1/x) snabbare.

Lyckligtvis, Euler kom upp med ett sätt att göra just detta:

ARCTAN(1/x) = (x / (1 + x²)) + ((2 * x) / (3*(1+x²)²)) + ((2 * 4 * x) / (3*5*(1+x²)³)) + ((2 * 4 * 6 * x) / (3*5*7*(1+x²)⁴)) +...

N: te termen i denna serie ges av funktionen f(n):

f(n) = f(n-1) * (2 * n) / ((2*n+1)*(1+x²))

där den första termen är (x / (1 + x²))

För att göra koden springa snabbare kan vi beräkna 1 + x² innan slingan så vi bara har att beräkna det en gång. Uppdaterade koden (gauss_pi_method_accelerated_arctan.py) är:

från decimal import *

def gauss_pi_method():
Return 4*(12*arctan(18) + 8*arctan(57) - 5*arctan(239))

def arctan(x): #defins en funktion för att definiera pi med formeln vi härrör

noggrannhet = Decimal(Decimal(1)/Decimal(10**(getcontext().prec-3)))
xSquaredPlusOne = Decimal (x ** 2 + 1)
resultat = Decimal(Decimal(x)/Decimal(xSquaredPlusOne)) #result är en variabel används för att lagra våra nuvarande värdet på pi som ett decimaltal, 1 och x förvandlas till decimaler att undvika avrundningsfel
fract = resultat
numurator_num = 2 # även num används till multiplie nämnaren
denominator_num = 3

samtidigt sant:
fract = (fract*numurator_num)/(Decimal(denominator_num)*xSquaredPlusOne)
resultatet += fract
om fract < noggrannhet:

Break
numurator_num += 2
denominator_num += 2

returnera resultat

medan 1:
#Asks användaren hur många decimaler som de vill ha det svara för att ges till
decimaler = int (input ("Ange hur många decimaler vill du ha svaret ges till:")) + 1

getcontext () .prec = decimaler + 4 #sets precision med 4 fler decimaler än bad om.

PiEst = gauss_pi_method() #computes pi

getcontext () .prec = decimaler #this kodblock rundor pi antal decimaler som bad om
PiEst = + PiEst

#next line skriver ut resultaten
skriva ut ("Pi till" + str(decimals-1) + "decimaler är:" + str(PiEst))

Denna kod tog strax över en halv sekund att beräkna π till 10.000 siffror, som är ungefär 30 gånger snabbare än innan!

Det finns ett snyggt trick att göra koden köra ännu snabbare. Hittills har vi använt decimal biblioteket i Python. Om vi gör beräkningar med hela heltal kommer det dock vara mycket snabbare. För att göra detta vi först multipliceras startvärdet med en stor makt på 10 och senare, när vi vill använda resultatet, vi delar det med den samma makten 10. Här är hur koden ska se ut (gauss_pi_method_fixed_point.py):

från decimal import *
Importera tid

def gauss_pi_method():
Return 4*(12*arctan(18) + 8*arctan(57) - 5*arctan(239))

def arctan(x):
en = 10 ** getcontext () .prec

xSquaredPlusOne = (x * x) + 1
fract = (x * en) / / xSquaredPlusOne
Total = fract
numurator_num = 2
medan 1:
nämnaren = (numurator_num + 1) * xSquaredPlusOne
fract * = numurator_num
fract = fract / / nämnare
om fract == 0:
Break
totala += fract
numurator_num += 2
återgå Decimal(Decimal(total)/Decimal(one))

medan 1:
#Asks användaren hur många decimaler som de vill ha det svara för att ges till
decimaler = int (input ("Ange hur många decimaler vill du ha svaret ges till:")) + 1
start_time = time.time()
getcontext () .prec = decimaler + 5 #sets precision med 4 fler decimaler än bad om.

PiEst = gauss_pi_method() #computes pi

getcontext () .prec = decimaler #this kodblock rundor pi antal decimaler som bad om
PiEst = + PiEst

end_time = time.time()
#next line skriver ut resultaten
skriva ut ("Pi till" + str(decimals-1) + "decimaler är:" + str(PiEst))

skriva ut ("algoritmen sprang i %s sekunder" % (time.time() - starttid))

Denna version av koden beräknas 10 000 decimaler av π i bara 0.26 sekunder, som är nästan dubbelt så snabb som den tidigare metoden!