Kalendarz

<< Grudzień 2019 >>
Po Wt Śr Cz Pi So Ni
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31          

OstatnieFoto

Partnerzy

Astro-Miejsca


URANIA

100 lat IAU

IAU

Centrum Nauki Kepler

Planetarium Wenus

ERC

Centrum Nauk Przyrodniczych

Orion,serwis,astronomii,PTA

POLSA

Astronomia Nova

Astronarium

forum astronomiczne

IPCN

Portal AstroNet

Puls Kosmosu

Forum Meteorytowe

kosmosnautaNET

kosmosnautaNET

Nauka w Polsce

astropolis

astromaniak

PTMA

PTR

heweliusz

heweliusz

ESA

Astronomers Without Borders

Hubble ESA

Space.com

Space Place

Instructables

Tu pełno nauki

Konkursy

Olimpiady Astronomiczne
Olimpiada Astronomiczna przebiega w trzech etapach.
Zadania zawodów I stopnia są rozwiązywane w warunkach pracy domowej. Zadania zawodów II i III stopnia mają charakter pracy samodzielnej. Zawody finałowe odbywają się w Planetarium Śląskim. Tematyka olimpiady wiąże ze sobą astronomię, fizykę i astronomiczne aspekty geografii. Olimpiady Astronomiczne


Urania Postępy Astronomii - konkurs dla szkół


astrolabium

Organizatorem konkursu astronomicznego jest Fundacja dla Uniwersytetu Jagiellońskiego a patronat nad akcją sprawuje Obserwatorium Astronomiczne im. Mikołaja Kopernika będące instytutem Wydziału Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego w Krakowie.
Zobacz szczegóły »

astrolabium

konkurs, astronomiczny

AstroSklepy

Astro Schopy
Uniwersał

Planeta Oczu

Astrocentrum

Aktualnie online

> Gości online: 4

> Użytkowników online: 0

> Łącznie użytkowników: 1
> Najnowszy użytkownik: jacek

Odwiedziny gości

Dziś:331
Wczoraj:3,390
W tym tygodniu:14,148
W tym miesiącu:22,236
W tym roku:701,430
Wszystkich:13,733,217

Ankieta

Gdzie jest Nowa Kelpera?

Lew

LMC

Rak

Wężownik

Smok

Rak

Pogoda F

Wszystko o Nas

Logo SA GW, autor Jacek Patka

Forum Astronomiczne PL


BOINC

Classroom

FB

Słoneczny panel

>Dziś jest:

Wschód słońca: 7:57
Zachód słońca: 15:51
>Dzień trwa:
7 Godzin 53 minut
Jest krótszy od najdłuższego dnia o: 10:43
Dane dla:
Żagań
Szerokość: 51°37 N
Długość: 15°19 E
Imieniny obchodzą:
Cecylia, Walerian, Mścigniew, Wolimir, Drogosława

Księżyc


Data: 15-12-2019 07:05:03

faza

Słońce

Na niebie


Mapa Nieba

TheSkyLive

CALSKY

Położenie ISS
The current position of the ISS
tranzyty ISS

Misja KEPLER

ZOONIVERSE odkrywanie planet

EPUP
4106 planet

Astropogoda

Pogoda


sat24, chmury, pogoda

III Prawo Keplera




Czytelnia


vademecum, miłośnika, astronomii, dwumiesięcznik, astronomia

Urania, numery archiwalne,przedwojenne

Light Pollution

M-WiFi

gwiazdy,zmienne,poradnik,gazeta,pdf,astronomia,pomiary

vademecum, miłośnika, astronomii, dwumiesięcznik, astronomia

astronomia amatorska

KTW'

Astronautilius

KTW'

kreiner, ziemia i wszechświat

kreiner, ziemia i wszechświat

poradnik, miłośnika, astronomii, książka, Tomasz, Rożek

poradnik, miłośnika, astronomii, książka, Rudż, Przemysław

atlas, nieba, książka, astronomia

atlas, księżyca, książka, astronomia

Poradnik Miłośnika Astronomii

Mądre Książki

Obliczenia liczby Pi

Obliczanie liczby Pi 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:
Z biegiem lat uzyskiwano coraz lepsze przybliżenia wartości π sięgające kilkuset miejsc po przecinku. W 1853 William Rutherford podał liczbę Pi z dokładnością 440 miejsc po przecinku. Rekordzistą w ręcznych obliczeniach liczby Pi jest William Shanks, któremu w 1874 udało się uzyskać 707 miejsc po przecinku. Zajęło mu to 15 lat. Później okazało się, że 180 ostatnich cyfr obliczył błędnie (wynik, który uznano za prawidłowy uwzględnia 527 miejsc po przecinku). W 1946 roku Ferguson podał wartość π do 620. miejsca po przecinku. W końcowych obliczeniach wspomagał się już kalkulatorem. Od 1949, kiedy to przy pomocy komputera ENIAC obliczono 2038 miejsc po przecinku, dokładniejsze aproksymacje liczby π uzyskiwano już tylko przy użyciu komputerów. We wrześniu 1999 roku obliczono π z dokładnością 2,0615·1011 miejsc po przecinku. Dokonał tego Takahasi przy pomocy komputera Hitachi SR8000 - Link.
31 grudnia 2009 r. Fabrice Bellard ogłosił, że udało mu się obliczyć π z dokładnością do 2 700 miliardów cyfr. Obliczenia ze sprawdzeniem zajęły 131 dni, do obliczeń użyto komputera z procesorem Intel Core i7 (2,93 GHz) i 6 GB RAM. Sam zapis binarny liczby zajmuje około 1,12 TB.
W roku 2010 obliczono cyfrę będącą na 2 000 000 000 000 000 miejscu po przecinku w rozwinięciu dziesiętnym liczby pi i wynosi ona zero. Obliczenia trwały 23 dni na 1000 maszynach.
W październiku 2011 Alexander J. Yee i Shigeru Kondo uzyskali dokładność ok. 10 bilionów (1013) miejsc po przecinku. Obliczenia zajęły 371 dni.
W październiku 2014 anonimowa osoba o nicku houkouonchi uzyskała dokładność ok. 13,3 bilionów miejsc po przecinku. Obliczenia zajęły 208 dni, a sprawdzanie 182 godziny.
W listopadzie 2016 Peter Trueb uzyskał dokładność ok. 22,5 bilionów miejsc po przecinku przy pomocy programu y-cruncher. Obliczenia zajęły 105 dni, a sama liczba zajęła ok. 120 TB miejsca.
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]))
Brak komentarzy. Może czas dodać swój?

Dodaj komentarz

Zaloguj się, aby móc dodać komentarz.

Oceny

Tylko zarejestrowani użytkownicy mogą oceniać zawartość strony
Zaloguj się , żeby móc zagłosować.

Brak ocen. Może czas dodać swoją?
21,818,014 unikalne wizyty