优化 pi 分数估计器

问题描述 投票:0回答:3

我见过 pi 估计值 22/7 和 355/113,我想我会制作一个程序来找到这样的分数。你给它一个特异性,它将数字除以它并找到最接近 pi 的值。

代码如下:

from math import pi, inf, ceil
from time import time, gmtime, strftime

print("\033c", end="")

specificity = int(input("Specificity? "))

fraction = [0, 0]

min_diff = inf

progress_print_interval = 0.05

start_time = time()

last_progress_print = time()

print()
print()

for i in range(specificity):
    for j in range(ceil(i / 3.1)):
        if not i or not j:
            continue
        if abs(i / j - pi) < min_diff:
            min_diff = abs(i / j - pi)
            fraction = [i, j]
        if (time() > last_progress_print + progress_print_interval):
            last_progress_print = time()

            total = i * (i / 3.1) + j / (i / 3.1)

            percent = total / (specificity * (specificity / 3.1))

            print("\033[F", end="")
            print("\033[F", end="")
            print(f"{(percent * 100):2f}%")
            print("Estimated time remaining:", strftime(f'%H:%M:%S', gmtime((time() - start_time) / percent - (time() - start_time))))

print("\033c", end="")
print("Done!")
print(f"{fraction[0]}/{fraction[1]}: {fraction[0] / fraction[1]}")

# best estimation so far: 833719/265381: 3.141592653581078

超过 1,000,000 的特异性需要数小时才能完成。有什么方法可以优化它吗?

python pi pypy
3个回答
1
投票

当然,这里有 5 个主要的低效率,

  1. 您使用两个嵌套循环来迭代分子和分母
  2. 没有算法优化。例如扩展或二进制搜索
  3. 您应该使用模块化算法来减少 Pi 的最佳近似值的搜索空间,这也可能导致大特异性值的性能下降。
  4. min_diff 没有初始化为足够大的值,所以第一个近似值总是被认为是最好的。
  5. ceil(i / 3.1) 不是最优上限

这是一个将所有这些都考虑在内的改进

import math

def pi_continued_fraction():
    """Generates the continued fraction expansion of Pi."""
    a = math.floor(math.pi)
    yield a
    while True:
        r = 1 / (math.pi - a)
        a = math.floor(r)
        yield a

def best_fraction(specificity):
    """Finds the best rational approximation of Pi up to a given specificity."""

    # use the continued fraction expansion of Pi to generate the best approximations
    cf = pi_continued_fraction()
    p, q = 0, 1
    for i in range(specificity):
        a = next(cf)
        p, q = q, a*q + p
        if math.gcd(p, q) != 1:
            continue
        approx = p / q
        if abs(math.pi - approx) < 1e-10:
            return p, q

    # If no match is found, perform bin search

    #  Use binary search to find the closest approximation
    lo_p, lo_q, hi_p, hi_q = 0, 1, 1, 0
    while True:
        mid_p, mid_q = lo_p + hi_p, lo_q + hi_q
        mid_approx = mid_p / mid_q
        if abs(math.pi - mid_approx) < 1e-10:
            return mid_p, mid_q
        elif mid_approx > math.pi:
            hi_p, hi_q = mid_p, mid_q
        else:
            lo_p, lo_q = mid_p, mid_q

# Example usage:
p, q = best_fraction(100000)
print(f"{p}/{q}: {p/q}")

0
投票

瞬间计算出所有越来越准确的分数:

3
13/4
16/5
19/6
22/7
179/57
201/64
223/71
245/78
267/85
289/92
311/99
333/106
355/113
52163/16604
52518/16717
52873/16830
53228/16943
53583/17056
53938/17169
54293/17282
54648/17395
55003/17508
55358/17621
55713/17734
56068/17847
56423/17960
56778/18073
57133/18186
57488/18299
57843/18412
58198/18525
58553/18638
58908/18751
59263/18864
59618/18977
59973/19090
60328/19203
60683/19316
61038/19429
61393/19542
61748/19655
62103/19768
62458/19881
62813/19994
63168/20107
63523/20220
63878/20333
64233/20446
64588/20559
64943/20672
65298/20785
65653/20898
66008/21011
66363/21124
66718/21237
67073/21350
67428/21463
67783/21576
68138/21689
68493/21802
68848/21915
69203/22028
69558/22141
69913/22254
70268/22367
70623/22480
70978/22593
71333/22706
71688/22819
72043/22932
72398/23045
72753/23158
73108/23271
73463/23384
73818/23497
74173/23610
74528/23723
74883/23836
75238/23949
75593/24062
75948/24175
76303/24288
76658/24401
77013/24514
77368/24627
77723/24740
78078/24853
78433/24966
78788/25079
79143/25192
79498/25305
79853/25418
80208/25531
80563/25644
80918/25757
81273/25870
81628/25983
81983/26096
82338/26209
82693/26322
83048/26435
83403/26548
83758/26661
84113/26774
84468/26887
84823/27000
85178/27113
85533/27226
85888/27339
86243/27452
86598/27565
86953/27678
87308/27791
87663/27904
88018/28017
88373/28130
88728/28243
89083/28356
89438/28469
89793/28582
90148/28695
90503/28808
90858/28921
91213/29034
91568/29147
91923/29260
92278/29373
92633/29486
92988/29599
93343/29712
93698/29825
94053/29938
94408/30051
94763/30164
95118/30277
95473/30390
95828/30503
96183/30616
96538/30729
96893/30842
97248/30955
97603/31068
97958/31181
98313/31294
98668/31407
99023/31520
99378/31633
99733/31746
100088/31859
100443/31972
100798/32085
101153/32198
101508/32311
101863/32424
102218/32537
102573/32650
102928/32763
103283/32876
103638/32989
103993/33102
104348/33215
208341/66317
312689/99532
833719/265381
1146408/364913
3126535/995207
4272943/1360120
5419351/1725033
42208400/13435351
47627751/15160384
53047102/16885417
58466453/18610450
63885804/20335483
69305155/22060516
74724506/23785549
80143857/25510582
165707065/52746197
245850922/78256779
571845701/182024140
817696623/260280919
1881244168/598818617
2698940791/859099536
7279125750/2317017689
9978066541/3176117225
22655073873/7211333986
32633140414/10387451211
140510628197/44725922069
173143768611/55113373280
205776909025/65500824491
238410049439/75888275702
509453239292/162164002615
747863288731/238052278317
1257316528023/400216280932
3262496344777/1038484840181
4519812872800/1438701121113
5777129400823/1838917402045
10296942273623/3277618523158
231052542892506/73546308630589
241349485166129/76823927153747
251646427439752/80101545676905
261943369713375/83379164200063
272240311986998/86656782723221
282537254260621/89934401246379
292834196534244/93212019769537
303131138807867/96489638292695
313428081081490/99767256815853
323725023355113/103044875339011
334021965628736/106322493862169
344318907902359/109600112385327
354615850175982/112877730908485
364912792449605/116155349431643
375209734723228/119432967954801
385506676996851/122710586477959
395803619270474/125988205001117
406100561544097/129265823524275
416397503817720/132543442047433
426694446091343/135821060570591
436991388364966/139098679093749
447288330638589/142376297616907
884279719003555/281474976710656
(884279719003555, 281474976710656)

最后一行是

math.pi
所代表的精确分数。

代码(在线尝试!):

from math import pi
from fractions import Fraction

fraction = 0
while fraction != pi:
    a, b = fraction.as_integer_ratio()
    while Fraction(pi).limit_denominator(b) == fraction:
        b *= 2
    lo, hi = b//2 + 1, b
    while lo < hi:
        mid = (lo + hi) // 2
        if Fraction(pi).limit_denominator(mid) == fraction:
            lo = mid + 1
        else:
            hi = mid
    fraction = Fraction(pi).limit_denominator(lo)
    print(fraction)

print(pi.as_integer_ratio())

-5
投票

瞬间生成 3126535/995207:

from math import pi
from fractions import Fraction

print(Fraction(pi).limit_denominator(10**6))

在线尝试!

© www.soinside.com 2019 - 2024. All rights reserved.