利用马青公式输出π的后任意位数字

马青公式

π = 16 a r c t a n 1 5 4 a r c t a n 1 239
a r c t a n x = x x 3 3 + x 5 5 . . . . . .
π = 16 × ( 1 1 × 5 1 3 × 5 3 + . . . . . . ) 4 × ( 1 1 × 239 1 3 × 239 3 + . . . . . . )
级数中的分数,分母增长很快,但我们可以对一个分式,不断除以同一个低精度数( 5 2 239 2 ),就能得到所有分数的值。


# Python 实现

import pandas as pd
def f(number):    
    number1 = number+10

    # 算到小数点后number1位
    b = 10**number1

    # 求含4/5的首项
    x1 = b*4//5
    # 求含1/239的首项
    x2 = b// -239

    # 求第一大项
    he = x1+x2
    #设置下面循环的终点,即共计算n项
    number *= 2

    #循环初值=3,末值2n,步长=2
    for i in range(3,number,2):
        # 求每个含1/5的项及符号
        x1 //= -25
        # 求每个含1/239的项及符号
        x2 //= -57121
        # 求两项之和
        x = (x1+x2) // i
        # 求总和
        he += x

    # 求出π
    pai = he*4
    #舍掉后十位
    pai //= 10**10

    ############ 输出圆周率π的值
    result=str(pai)
    return result

最后为了方便统计出现的数字的次数,将其转换为series格式,并定义以下函数来可视化。

def show(n,plt_type='pie'):
    result = f(n)
    a = pd.Series([int(i) for i in result])
    b = a.value_counts()
    if plt_type=='pie':
        plt.figure(figsize=(10,8))
        plt.pie(b,labels=b.index,explode=[0.05]*10,shadow=True,autopct='%1.1f%%')
        plt.show()
    elif plt_type=='bar':
        plt.figure(figsize=(10,8))
        plt.bar(b.index,b.values)
        plt.show()
    else:
        print('Type Wrong')

利用此函数可输出饼图或柱形图
当n=100时:

show(100,'pie')
show(100,'bar')

这里写图片描述
这里写图片描述

当n=10000时:

show(10000,'bar')
show(10000,'pie')

这里写图片描述
这里写图片描述

由此来看 π 应该是正规数。突然想看看在其中会不会出现自己的生日,于是定义以下函数

import re
def find_(pattern,result):
    loc = result.find(pattern)
    return loc
result=f(300000)#产生300000位数
find_('199766',result)

最后return了140151,惊奇竟然真的在里面了!(虽然这个算法有误差hhh不过姑且当他正确吧)

猜你喜欢

转载自blog.csdn.net/qq_35719435/article/details/80848747