• Register

Kako sympy izračuna število pi?

+11 votes
1,013 views
Odgovor na to vprašanje je vreden 1xP!

Pravilni odgovor: najdete source definicije (def) kjer se računa Pi in tudi povezavo na teoretično ozadje (npr. na Wikipediji!

Če noben pravilno ne ugotovi v cca 1 dnevu, bom dal namig spodaj z tag-i!
asked Mar 14, 2016 by janko.slavic (78,700 points)

4 Answers

+9 votes
 
Best answer

Po brskanju po datotekah v Sympy, natančneje:(Anaconda3\Lib\site-packages\sympy\mpmath\libmp\libelefun.py) sem našel, da se št. pi izračuna po algoritmu bratov Chudnovsky.

Torej:

def ln10_fixed(prec):
    """
    Computes ln(10). This is done with a hyperbolic Machin-type formula.
    """
    return machin([(46, 31), (34, 49), (20, 161)], prec, True)


"""
For computation of pi, we use the Chudnovsky series:

             oo
             ___        k
      1     \       (-1)  (6 k)! (A + B k)
    ----- =  )     -----------------------
    12 pi   /___               3  3k+3/2
                    (3 k)! (k!)  C
            k = 0

where A, B, and C are certain integer constants. This series adds roughly
14 digits per term. Note that C^(3/2) can be extracted so that the
series contains only rational terms. This makes binary splitting very
efficient.

The recurrence formulas for the binary splitting were taken from
ftp://ftp.gmplib.org/pub/src/gmp-chudnovsky.c

Previously, Machin's formula was used at low precision and the AGM iteration
was used at high precision. However, the Chudnovsky series is essentially as
fast as the Machin formula at low precision and in practice about 3x faster
than the AGM at high precision (despite theoretically having a worse
asymptotic complexity), so there is no reason not to use it in all cases.

"""

# Constants in Chudnovsky's series
CHUD_A = MPZ(13591409)
CHUD_B = MPZ(545140134)
CHUD_C = MPZ(640320)
CHUD_D = MPZ(12)

def bs_chudnovsky(a, b, level, verbose):
    """
    Computes the sum from a to b of the series in the Chudnovsky
    formula. Returns g, p, q where p/q is the sum as an exact
    fraction and g is a temporary value used to save work
    for recursive calls.
    """
    if b-a == 1:
        g = MPZ((6*b-5)*(2*b-1)*(6*b-1))
        p = b**3 * CHUD_C**3 // 24
        q = (-1)**b * g * (CHUD_A+CHUD_B*b)
    else:
        if verbose and level < 4:
            print("  binary splitting", a, b)
        mid = (a+b)//2
        g1, p1, q1 = bs_chudnovsky(a, mid, level+1, verbose)
        g2, p2, q2 = bs_chudnovsky(mid, b, level+1, verbose)
        p = p1*p2
        g = g1*g2
        q = q1*p2 + q2*g1
    return g, p, q

@constant_memo
def pi_fixed(prec, verbose=False, verbose_base=None):
    """
    Compute floor(pi * 2**prec) as a big integer.

    This is done using Chudnovsky's series (see comments in
    libelefun.py for details).
    """
    # The Chudnovsky series gives 14.18 digits per term
    N = int(prec/3.3219280948/14.181647462 + 2)
    if verbose:
        print("binary splitting with N =", N)
    g, p, q = bs_chudnovsky(0, N, 0, verbose)
    sqrtC = isqrt_fast(CHUD_C<<(2*prec))
    v = p*CHUD_C*sqrtC//((q+CHUD_A*p)*CHUD_D)
    return v
answered Mar 16, 2016 by jan_pestotnik (640 points)
selected Mar 16, 2016 by janko.slavic
Tule je še koda od Sympy:

class Pi(with_metaclass(Singleton, NumberSymbol)):
    r"""The `\pi` constant.

    The transcendental number `\pi = 3.141592654\dots` represents the ratio
    of a circle's circumference to its diameter, the area of the unit circle,
    the half-period of trigonometric functions, and many other things
    in mathematics.

    Pi is a singleton, and can be accessed by ``S.Pi``, or can
    be imported as ``pi``.

    Examples
    ========

    >>> from sympy import S, pi, oo, sin, exp, integrate, Symbol
    >>> S.Pi
    pi
    >>> pi > 3
    True
    >>> pi.is_irrational
    True
    >>> x = Symbol('x')
    >>> sin(x + 2*pi)
    sin(x)
    >>> integrate(exp(-x**2), (x, -oo, oo))
    sqrt(pi)

    References
    ==========

    .. [1] http://en.wikipedia.org/wiki/Pi
    """

    is_real = True
    is_positive = True
    is_negative = False
    is_irrational = True
    is_number = True
    is_algebraic = False
    is_transcendental = True

    __slots__ = []

    def _latex(self, printer):
        return r"\pi"

    @staticmethod
    def __abs__():
        return S.Pi

    def __int__(self):
        return 3

    def _as_mpf_val(self, prec):
        return mpf_pi(prec)

    def approximation_interval(self, number_cls):
        if issubclass(number_cls, Integer):
            return (Integer(3), Integer(4))
        elif issubclass(number_cls, Rational):
            return (Rational(223, 71), Rational(22, 7))

    def _sage_(self):
        import sage.all as sage
        return sage.pi
pi = S.Pi
Odlično!

To je morebiti najbolj popolni odgovor... torej sympy st tukaj opre na mpmath itd...

PS: upazil sem, da je Gregor vprašal tudi na stack in je odgovoril avtor kode na mpmath:
http://stackoverflow.com/questions/36004680/how-does-sympy-calculate-pi/36018371#36018371
+1 vote
Python izračuna pi s funkcijo evalf (pi.evalf(n= 'št.decimalnih mest') ali s funkcijo N(pi, 'število decimalnih mest').

Natančneje: evalvira dano formulo na natančno vrednost zapisano na n decimalnih mest.
answered Mar 14, 2016 by iztok.jerman (230 points)
Dober poskus, vendar je to odgovor na vprašanje kako izpišem pi z željenim številom decimalnih mest.
+3 votes

SymPy ne shranjuje Pi vrednosti kot konstanto/float (ker je neskončno število , ampak Pi shrani kot objekt, kateri vsebuje konstanto... Medtem, ko na primer numpy vrednosti Pi so le približno zaokrožene vrednosti...

SymPy uses mpmath in the background, which makes it possible to perform computations using arbitrary-precision arithmetic. That way, some special constants, like e, pi, oo (Infinity), are treated as symbols and can be evaluated with arbitrary precision

Tudi jaz sem mnenja, da ga izračunamo s  funkcijo evalf.

pi.evalf()

 

answered Mar 14, 2016 by urskamlakar (4,080 points)
Dober poskus, gre v pravo smer, vendar me zanima "kako" določi poljubno število mest? kako to naredi numerično? kaj je zadaj?
Pri funkciji evalf je privzeto numerično vrednotenje z natančnostjo 15 decimalnih mest. Bolj natančno pa ne bi vedela..
+4 votes

Sympy upprablja mpmath v ozadju in definicija števila pi je : 

V tem primeru izračuna pi na 1000 mest natančno

Še povezava na stran kjer sem dobil tole - http://stackoverflow.com/questions/9004789/1000-digits-of-pi-in-python

answered Mar 15, 2016 by jans (3,290 points)
Še vsa teorija števila pi - https://en.m.wikipedia.org/wiki/Pi
Super! Zelo blizu; lahko najdete še definicijo funkcije v sympy oz mpmath?
Mpmath je brezplačna knjižnica v Pythonu in se uporablja za realna in kompeksna števila z veliko natančnostjo. Brez težav izračunamo rezultat na 10 mest natančno ali pa na 1000. http://mpmath.org

Definicija ki jo dobim pod help v jupyter notebooku - help(mpmath.pi)

class constant(_constant)
 |  Represents a mathematical constant with dynamic precision.
 |  When printed or used in an arithmetic operation, a constant
 |  is converted to a regular mpf at the working precision. A
 |  regular mpf can also be obtained using the operation +x.
 |  
 |  Method resolution order:
 |      constant
 |      _constant
 |      _mpf
 |      mpnumeric
 |      builtins.object
Zaslužite si 1xP, vendar to ni povsem pravi odgovor. Sympy je open source in vsa koda je dostopna; kdo prvi najde dejansko kodo, ki izračuna še vedno dobi 1xP! Ne iščem torej načelne kode, ampak katera koda se izvede! Namig: sledite source kodi sympy.pi!
...