jacek 12/12/2017
• Programowanie
• 2217 czytań
0 komentarzy

Być może niedościgniemy najlepszych w ilości uzyskanych cyfr Ludolfiny, ale to dobry test dla komputera i nas samych. Przy okazji poznamy kawał dobrej matematyki

#!/bin/bin/nv python
#-*- coding: utf8 -*-

'''Chudnovsky algorithm
https://en.wikipedia.org/wiki/Chudnovsky_algorithm'''

from decimal import Decimal as Dec, getcontext as gc

def PI(maxK=70, prec=2001, disp=2000): # parameter defaults chosen to gain 1000+ digits within a few seconds
    gc().prec = prec
    K, M, L, X, S = 6, 1, 13591409, 1, 13591409 
    for k in range(1, maxK+1):
        M = (K**3 - (K<<4)) * M / k**3 
        L += 545140134
        X *= -262537412640768000
        S += Dec(M * L) / X
        K += 12
    pi = 426880 * Dec(10005).sqrt() / S
    pi = Dec(str(pi)[:disp]) # drop few digits of precision for accuracy
    print("PI(maxK=%d iterations, gc().prec=%d, disp=%d digits) =\n%s" % (maxK, prec, disp, pi))
    return pi

Pi = PI()
print("\nFor greater precision and more digits (takes a few extra seconds) - Try")
print("Pi = PI(317,4501,4500)")
print("Pi = PI(353,5022,5020)")

kolejny kod wykorzystuje formułę Machina:

#!/usr/bin/nv python
# -*- coding: utf-8 -*-

'''
Python 2.7.3

Today, computers can be used to calculate millions of digits of pi by using
infinite series formulas.  There are several different ways to do this,
but one of the most efficient methods that is relatively straightforward
is to combine Machin’s formula:

http://lostmathlessons.blogspot.com/2016/03/using-your-computer-to-calculate-pi.html
'''
# Pi Calculator
# Python 2.7.3
# After running, type "pi(n)" where n is the number of decimals for pi.  For
#  example, if you would like to calculate 100 decimals for pi, type "pi(100)".
# import python libraries
from decimal import Decimal, getcontext
from time import time, strftime
import datetime
# arccot function using power formula arccot = 1/x - 1/(3x^3) + 1/(5x^5) ...
def arccot(x, digits):
    # set precision and starting values
    getcontext().prec = digits
    total = 0
    n = 1

    # loop while new term is large enough to actually change the total
    while Decimal((2 * n - 1) * x ** (2 * n - 1)) < Decimal(10 ** digits):
        # find value of new term
        term = ((-1)**(n - 1)) * 1 / Decimal((2 * n - 1) * x ** (2 * n - 1))
        # add the new term to the total
        total += term
        # next n
        n += 1
    # return the sum
    return total

# pi function
def pi(decimals):
    # start timer
    timestart = time()
    # find pi using Machin's Formula pi = 4 * (4 * arccot(5) - arccot(239))
    #  and the power formula for arccot (see arccot function above)
    print "pi = "+ str(Decimal(4*(4*arccot(5,decimals+3)-arccot(239,decimals+3))).quantize(Decimal(10)**(-decimals)))
    # display elapsed time
    timeelapsedint = round(time()-timestart,2)
    timeelapsedstr = str(datetime.timedelta(seconds = round(timeelapsedint, 0)))
    print "runtime: " + timeelapsedstr + " or " + str(timeelapsedint) + " seconds."
pi(1000)

Używamy modułu wieloprocesorowego Pythona do wykonania długiego zadania w oddzielnym procesie. Ten pakiet obsługuje procesy odradzania przy użyciu API podobnego do modułu wątków. Moduł przetwarzania wieloprocesowego działa inaczej w systemach Windows i Linux. W systemie Linux wywoływane jest wywołanie fork(), natomiast w systemie Windows używany jest moduł pickle. Ma to pewne ważne konsekwencje. Podzielimy przykład kodu na programy Linux i Windows.

Jest jedno ograniczenie naszego programu. Używa obiektu Queue do wymiany wartości między procesami. Ta kolejka ma limit rozmiaru wartości, którą może otrzymać. Rozmiar jest określany przez potok OS. Po przekroczeniu tego limitu program kończy się impasem. Największe Pi obliczone na Windows 7 miało 8156 cyfr. W systemie Linux udało nam się obliczyć ponad 60 000 cyfr Pi.

Osobiście uważam, i widzę z powyższych doświadczeń, że poniższy kod jest najmniej wydajny. Tzn. Powyższe dwa kody znacznie szybciej wyliczą zdecydowanie większa ilość cyfr liczby Pi w znacznie krótszym czasie. Ale program ma okienka ... bez komentarza.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
ZetCode Tkinter e-book

This script produces a long-running task of calculating 
a large Pi number, while keeping the GUI responsive.

Author: Jan Bodnar
Last modified: January 2016
Website: www.zetcode.com

http://zetcode.com/articles/tkinterlongruntask/
"""

from tkinter import (Tk, BOTH, Text, E, W, S, N, END, 
    NORMAL, DISABLED, StringVar)
from tkinter.ttk import Frame, Label, Button, Progressbar, Entry
from tkinter import scrolledtext

from multiprocessing import Queue, Process
import queue 
from decimal import Decimal, getcontext

DELAY1 = 80
DELAY2 = 20

class Example(Frame):
  
    def __init__(self, parent, q):
        Frame.__init__(self, parent)   
         
        self.queue = q 
        self.parent = parent
        self.initUI()
        
        
    def initUI(self):
      
        self.parent.title("Pi computation")
        self.pack(fill=BOTH, expand=True)
        
        self.grid_columnconfigure(4, weight=1)
        self.grid_rowconfigure(3, weight=1)
        
        lbl1 = Label(self, text="Digits:")
        lbl1.grid(row=0, column=0, sticky=E, padx=10, pady=10)
        
        self.ent1 = Entry(self, width=10)
        self.ent1.insert(END, "4000")
        self.ent1.grid(row=0, column=1, sticky=W)
        
        lbl2 = Label(self, text="Accuracy:")
        lbl2.grid(row=0, column=2, sticky=E, padx=10, pady=10)

        self.ent2 = Entry(self, width=10)
        self.ent2.insert(END, "100")
        self.ent2.grid(row=0, column=3, sticky=W)        
        
        self.startBtn = Button(self, text="Start", 
            command=self.onStart)
        self.startBtn.grid(row=1, column=0, padx=10, pady=5, sticky=W)
        
        self.pbar = Progressbar(self, mode='indeterminate')        
        self.pbar.grid(row=1, column=1, columnspan=3, sticky=W+E)     
        
        self.txt = scrolledtext.ScrolledText(self)  
        self.txt.grid(row=2, column=0, rowspan=4, padx=10, pady=5,
            columnspan=5, sticky=E+W+S+N)
       
        
    def onStart(self):
        
        self.startBtn.config(state=DISABLED)
        self.txt.delete("1.0", END)
        
        self.digits = int(self.ent1.get())
        self.accuracy = int(self.ent2.get())
        
        self.p1 = Process(target=self.generatePi, args=(self.queue,))
        self.p1.start()
        self.pbar.start(DELAY2)
        self.after(DELAY1, self.onGetValue)
        
       
    def onGetValue(self):
        
        if (self.p1.is_alive()):
            
            self.after(DELAY1, self.onGetValue)
            return
        else:    
        
            try:
                self.txt.insert('end', self.queue.get(0))
                self.txt.insert('end', "\n")
                self.pbar.stop()
                self.startBtn.config(state=NORMAL)
            
            except queue.Empty:
                print("queue is empty")
            
            
    def generatePi(self, queue):
        
        getcontext().prec = self.digits
        
        pi = Decimal(0)
        k = 0
        n = self.accuracy
        
        while k < n:
            pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1)) - \
                (Decimal(2)/(8*k+4)) - (Decimal(1)/(8*k+5))- \
                (Decimal(1)/(8*k+6)))
            k += 1
            print (self.p1.is_alive())
            
        queue.put(pi)    
        print("end")    


def main():
    
    q = Queue()
  
    root = Tk()
    root.geometry("400x350+300+300")
    app = Example(root, q)
    root.mainloop()  


if __name__ == '__main__':
    main()  

Bardzo szybki programik. Efekty widać od razu. Warto skierować dane wynikowe do pliku.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# http://rosettacode.org/wiki/Pi
def calcPi():
    q, r, t, k, n, l = 1, 0, 1, 1, 3, 3
    while True:
        if 4*q+r-t < n*t:
            yield n
            nr = 10*(r-n*t)
            n  = ((10*(3*q+r))//t)-10*n
            q  *= 10
            r  = nr
        else:
            nr = (2*q+r)*l
            nn = (q*(7*k)+2+(r*l))//(t*l)
            q  *= k
            t  *= l
            l  += 2
            k += 1
            n  = nn
            r  = nr
 
import sys
pi_digits = calcPi()
i = 0
for d in pi_digits:
    sys.stdout.write(str(d))
    #i += 1
    #if i == 60: print(""); i = 0

Z obliczeniami Pi jest taki problem, że kolejne cyfry zajmują miejsce w pamięci komputerów, w tym dysków twardych. Poniższe pokazuje jak szybko zapełniamy tymi cyframi wszelkie pamięci. Więc nim się weźmiesz za liczenie na swoim komputerze pomyśl, czy to w ogóle jest możliwe: A teraz porównanie kilku algorytmów na Pi:

from decimal import *

#Sets decimal to 25 digits of precision
#http://blog.recursiveprocess.com/2013/03/14/calculate-pi-with-python/
getcontext().prec = 25

def factorial(n):
    if n<1:
        return 1
    else:
        return n * factorial(n-1)

def plouffBig(n): #http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6)))
        k += 1
    return pi

def bellardBig(n): #http://en.wikipedia.org/wiki/Bellard%27s_formula
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) -Decimal(1)/(4*k+3))
        k += 1
    pi = pi * 1/(2**6)
    return pi

def chudnovskyBig(n): #http://en.wikipedia.org/wiki/Chudnovsky_algorithm
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))
        k += 1
    pi = pi * Decimal(10005).sqrt()/4270934400
    pi = pi**(-1)
    return pi
print ("\t\t\t Plouff \t\t Bellard \t\t\t Chudnovsky")
for i in range(1,20):
    print ("Iteration number ",i, " ", plouffBig(i), " " , bellardBig(i)," ", chudnovskyBig(i))

I jeszcze dwa skrypty Python3 dal obliczania Pi:

"""
Python3 program to calculate Pi using python long integers, and the
Chudnovsky algorithm

See: http://www.craig-wood.com/nick/articles/pi-chudnovsky/ for more
info

Nick Craig-Wood 
"""

import math
from time import time

def sqrt(n, one):
    """
    Return the square root of n as a fixed point number with the one
    passed in.  It uses a second order Newton-Raphson convgence.  This
    doubles the number of significant figures on each iteration.
    """
    # Use floating point arithmetic to make an initial guess
    floating_point_precision = 10**16
    n_float = float((n * floating_point_precision) // one) / floating_point_precision
    x = (int(floating_point_precision * math.sqrt(n_float)) * one) // floating_point_precision
    n_one = n * one
    while 1:
        x_old = x
        x = (x + n_one // x) // 2
        if x == x_old:
            break
    return x

def pi_chudnovsky(one=1000000):
    """
    Calculate pi using Chudnovsky's series

    This calculates it in fixed point, using the value for one passed in
    """
    k = 1
    a_k = one
    a_sum = one
    b_sum = 0
    C = 640320
    C3_OVER_24 = C**3 // 24
    while 1:
        a_k *= -(6*k-5)*(2*k-1)*(6*k-1)
        a_k //= k*k*k*C3_OVER_24
        a_sum += a_k
        b_sum += k * a_k
        k += 1
        if a_k == 0:
            break
    total = 13591409*a_sum + 545140134*b_sum
    pi = (426880*sqrt(10005*one, one)*one) // total
    return pi

if __name__ == "__main__":
    print(pi_chudnovsky(10**100))
    for log10_digits in range(1,7):
        digits = 10**log10_digits
        one = 10**digits

        start =time()
        pi = pi_chudnovsky(one)
        #print(pi)
        print("chudnovsky: digits",digits,"time",time()-start)

Poniższy wymaga zainstalowania biblioteki gmpy2 - Link na stronie autora skryptu znajdziemy linki do biblioteki oraz instrukcje. W Linux'ie dla Python'a najprościej zainstalować ją poprzez sudo pip3 install gmpy2. Gdyby pojawiły się problemy z instalacją, zapewne zabraknie bibliotek, to należy je doinstalować. Wcześniej przeczytać wygenerowany opis.

Skrypt ten ma już poprawiny zawarty w nim błąd. Więc nieco się różni od oryginału. Bibliotek GMPY2 pozwala na realizację obliczeń na dużych liczbach z odpowiednią precyzją. System i Python mają swoje ograniczenia. Wielkie liczby ... trzeba je obejść.

"""
Python3 program to calculate Pi using python long integers, binary
splitting and the Chudnovsky algorithm

See: http://www.craig-wood.com/nick/articles/pi-chudnovsky/ for more
info

Nick Craig-Wood 
https://www.craig-wood.com/nick/articles/pi-chudnovsky/
"""

import math
from gmpy2 import mpz, isqrt
from time import time

def pi_chudnovsky_bs(digits):
    """
    Compute int(pi * 10**digits)

    This is done using Chudnovsky's series with binary splitting
    """
    C = 640320
    C3_OVER_24 = C**3 // 24
    def bs(a, b):
        """
        Computes the terms for binary splitting the Chudnovsky infinite series

        a(a) = +/- (13591409 + 545140134*a)
        p(a) = (6*a-5)*(2*a-1)*(6*a-1)
        b(a) = 1
        q(a) = a*a*a*C3_OVER_24

        returns P(a,b), Q(a,b) and T(a,b)
        """
        if b - a == 1:
            # Directly compute P(a,a+1), Q(a,a+1) and T(a,a+1)
            if a == 0:
                Pab = Qab = mpz(1)
            else:
                Pab = mpz((6*a-5)*(2*a-1)*(6*a-1))
                Qab = mpz(a*a*a*C3_OVER_24)
            Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)
            if a & 1:
                Tab = -Tab
        else:
            # Recursively compute P(a,b), Q(a,b) and T(a,b)
            # m is the midpoint of a and b
            m = (a + b) // 2
            # Recursively calculate P(a,m), Q(a,m) and T(a,m)
            Pam, Qam, Tam = bs(a, m)
            # Recursively calculate P(m,b), Q(m,b) and T(m,b)
            Pmb, Qmb, Tmb = bs(m, b)
            # Now combine
            Pab = Pam * Pmb
            Qab = Qam * Qmb
            Tab = Qmb * Tam + Pam * Tmb
        return Pab, Qab, Tab
    # how many terms to compute
    DIGITS_PER_TERM = math.log10(C3_OVER_24/6/2/6)
    N = int(digits/DIGITS_PER_TERM + 1)
    # Calclate P(0,N) and Q(0,N)
    P, Q, T = bs(0, N)
    one_squared = mpz(10)**(2*digits)
    sqrtC = isqrt(10005*one_squared)
    return (Q*426880*sqrtC) // T

# The last 5 digits or pi for various numbers of digits
check_digits = {
        100 : 70679,
       1000 :  1989,
      10000 : 75678,
     100000 : 24646,
    1000000 : 58151,
   10000000 : 55897,
}

if __name__ == "__main__":
    digits = 100
    pi = pi_chudnovsky_bs(digits)
    print(pi)
    #raise SystemExit
    for log10_digits in range(1,9):
        digits = 10**log10_digits
        start =time()
        pi = pi_chudnovsky_bs(digits)
        print("chudnovsky_gmpy_mpz_bs: digits",digits,"time",time()-start)
        if digits in check_digits:
            last_five_digits = pi % 100000
            if check_digits[digits] == last_five_digits:
                print("Last 5 digits %05d OK" % last_five_digits)
            else:
                print("Last 5 digits %05d wrong should be %05d" % (last_five_digits, check_digits[digits]))