我见过 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 的特异性需要数小时才能完成。有什么方法可以优化它吗?
当然,这里有 5 个主要的低效率,
这是一个将所有这些都考虑在内的改进
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}")
瞬间计算出所有越来越准确的分数:
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())
瞬间生成 3126535/995207:
from math import pi
from fractions import Fraction
print(Fraction(pi).limit_denominator(10**6))