python/math/pi_continued_fraction_convergents.py

#!/usr/bin/env python
from functools import lru_cache
import math

"""
Continued Fraction Representation of Pi

This script computes the convergents
(rational representation) of some continued
fraction representations of Pi.

Happy Pi Day!

2019-03-14
"""

class GeneralContinuedFraction(object):
    """Abstract class representing a 
    general continued fraction.
    This class generates terms in the
    general continued fraction representation
    of a number:

        x = b0 +   a1
                 ------------
                 b1 +   a2
                      --------
                      b2 +  a3
                           ----
                            b3 + ...
    """
    def __init__(self,nterms):
        self.nterms = nterms

    def __len__(self):
        return self.nterms

    def a(self,k):
        raise NotImplementedError

    def b(self,k):
        raise NotImplementedError

    @classmethod
    def str_convergent(cls,n):
        n,d = cls.convergent(n)
        return "%d / %d"%(n,d)

    @classmethod
    def dec_convergent(cls,n):
        n,d = cls.convergent(n)
        return "%0.16f"%(n/d)

    @classmethod
    def tex_convergent(cls,n):
        n,d = cls.convergent(n)
        return "\dfrac{ %d }{ %d }"%(n,d)

    @classmethod
    @lru_cache(maxsize=64)
    def convergent(cls,n):
        """Compute and return the nth convergent for

        x = b0 + a1/( b1 + a2/( b2 + a3/( b3 + a4/( ... 

        by applying the recurrence formula:

        x0 = A0/B0 = b0
        x1 = A1/B1 = (b1 b0 + a1)/(b1)
        x2 = A2/B2 = (b2 (b1 b0 + a1) + a2 b0)/(b2 b1 + a2)

        ...

        where An is the numerator and Bn is the denominator,
        called continuants, given by the recursion:

        A_{n} = b_{n} A_{n-1} + a_{n} A_{n-2}
        B_{n} = b_{n} B_{n-1} + a_{n} B_{n-2}

        for n > 0 with init values

        A_{-1}  = 1
        A_{0}   = b0
        B_{-1}  = 0
        B_{0}   = 1
        """
        def A(k):
            if k==-1:
                return 1
            elif k==0:
                return cls.b(0)
            else:
                return cls.b(k)*A(k-1) + cls.a(k) * A(k-2)

        def B(k):
            if k==-1:
                return 0
            elif k==0:
                return 1
            else:
                return cls.b(k)*B(k-1) + cls.a(k) * B(k-2)

        return A(n), B(n)



class SimpleContinuedFraction(GeneralContinuedFraction):
    """A simple continued fraction is equivalent
    to a general continued fraction with all
    denominators (a_i) set to 1.
    """
    @classmethod
    def a(cls,k):
        return 1



class FourOverPi_ContinuedFraction(GeneralContinuedFraction):
    """Class representing Pi
    as a general continued fraction
    """
    @classmethod
    def a(cls,k):
        if k>0:
            return (2*k-1)*(2*k-1)

    @classmethod
    def b(cls,k):
        if k==0:
            return 1
        else:
            return 2


class Pi_ContinuedFraction(FourOverPi_ContinuedFraction):

    @classmethod
    def str_convergent(cls,n):
        # switch numerator and denominator
        d,n = cls.convergent(n)
        # multiply numerator by 4
        n *= 4
        return "%d / %d"%(n,d)

    @classmethod
    def dec_convergent(cls,n):
        # switch numerator and denominator
        d,n = cls.convergent(n)
        # multiply numerator by 4
        n *= 4
        return "%0.16f"%(n/d)

    @classmethod
    def tex_convergent(cls,n):
        # switch numerator and denominator
        d,n = cls.convergent(n)
        # multiply numerator by 4
        n *= 4
        return "\dfrac{ %d }{ %d }"%(n,d)




class Pi_ContinuedFraction2(GeneralContinuedFraction):
    """Class representing Pi
    as a general continued fraction
    (redux)
    """
    @classmethod
    def a(cls,k):
        if k>0:
            return (2*k-1)*(2*k-1)

    @classmethod
    def b(cls,k):
        if k==0:
            return 3
        else:
            return 6



class Pi_Simple(SimpleContinuedFraction):
    terms = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2]

    @classmethod
    def b(cls,k):
        return cls.terms[k]



if __name__=="__main__":

    print("%40s"%("Pi (pretty bad)"))

    print("-"*80)
    for i in range(1,24):
        print("%-60s\t%-18s"%(
            Pi_ContinuedFraction.str_convergent(i),
            Pi_ContinuedFraction.dec_convergent(i),
        ))

    print("\n\n")






    print("%40s"%("Pi (a bit worse)"))

    print("-"*80)
    for i in range(1,24):
        print("%-60s\t%-18s"%(
            Pi_ContinuedFraction2.str_convergent(i),
            Pi_ContinuedFraction2.dec_convergent(i),
        ))

    print("\n\n")




    print("%40s"%("Pi (simplest and best)"))
    print("-"*80)
    for i in range(1,14):
        print("%-60s\t%-18s"%(
            Pi_Simple.str_convergent(i),
            Pi_Simple.dec_convergent(i),
        ))

    print("\n\n")