r1ngs's blog

bbp算法求PI的小数点后第n位

2019-03-16

上周学长给我出了一道ACM题,以前没有打过ACM,但是通过这道题也是学到了东西的,干脆就记录一下

基本定理

$\left \lfloor \frac{p}{q} \right \rfloor mod\;m=\left \lfloor \frac{p\;mod\;mq}{q} \right \rfloor​$,这里的p和q都是整数。

啥?这个定理叫什么名字,怎么没有听说过?

这大概就是我和大佬的差距吧,自己去推了一下,不难,但是要想到这种做法,运用到后面的快速幂算法里,我觉得还是很不容易的:

令$p=kmq+r$,$0\leqslant r<mq$

$\left \lfloor \frac{p\;mod\;mq}{q} \right \rfloor=\left \lfloor \frac{r}{q}\right \rfloor=\left \lfloor \frac{p-kmq}{q} \right \rfloor=\left \lfloor \frac{p}{q}\right \rfloor-km$

由于$0\leqslant r<mq$,则$0\leqslant \frac{r}{ q}<m$所以$0\leqslant\left \lfloor \frac{p}{q}\right \rfloor-km<m$

故$\left \lfloor \frac{p}{q}\right \rfloor-km=\left \lfloor \frac{p}{q}\right \rfloor\;mod\;m​$

求小数点后第n位

如果要求$\frac{1}{a}​$小数点后n位可以这样求:$\left \lfloor \frac{10^{n}}{a}\right \rfloor\;mod\;10=\left \lfloor \frac{10^n\;mod\;10a}{a} \right \rfloor​$,而分子运用快速模幂算法,就能很快算出来

所以$\frac{1}{a}+\frac{1}{b}+\frac{1}{c}$小数点后第n位只需要将三者小数点后第n位相加然后求和,然后考虑第n位后面的低位是否存在进位的问题就可以了,我实现了用函数递归调用的求法,一直向低位探索直到不存在进位

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def my_pow(a, b, c):
a = a % c
res = 1

while b != 0:
if b & 1:
res = res * a
b >>= 1
a = (a ** 2) % c

return res % c

def add(a, b, c):

i = 10000
a_0 = my_pow(10, i, 10 * a) / a
b_0 = my_pow(10, i, 10 * b) / b
c_0 = my_pow(10, i, 10 * c) / c
s_0 = (a_0 + b_0 + c_0)%10

i += 1
a_t = my_pow(10, i, 10*a)/a
b_t = my_pow(10, i, 10*b)/b
c_t = my_pow(10, i, 10*c)/c
s_t = a_t+b_t+c_t

if s_t >= 20:
return (s_0 + 2) % 10
elif s_t == 18:
return (s_0 + if_18(a, b, c, i))%10
elif s_t == 19:
return (s_0 + if_19(a, b, c, i))%10
elif 10 <= s_t < 18:
return (s_0 + 1) % 10
elif s_t == 8:
return (s_0+if_8(a, b, c, i))%10
elif s_t == 9:
return (s_0+if_9(a, b, c, i))%10
else:
return s_0


def if_8(a, b, c, i):
a_0 = my_pow(10, i, 10 * a) / a
b_0 = my_pow(10, i, 10 * b) / b
c_0 = my_pow(10, i, 10 * c) / c
s_0 = a_0 + b_0 + c_0

i += 1
a_t = my_pow(10, i, 10*a)/a
b_t = my_pow(10, i, 10*b)/b
c_t = my_pow(10, i, 10*c)/c
s_t = a_t+b_t+c_t

if s_t>=20:
return 1
elif s_t == 18:
return (s_0+if_18(a, b, c, i))/10
elif s_t == 19:
return (s_0+if_19(a, b, c, i))/10
else :
return 0

def if_9(a, b, c, i):
a_0 = my_pow(10, i, 10 * a) / a
b_0 = my_pow(10, i, 10 * b) / b
c_0 = my_pow(10, i, 10 * c) / c
s_0 = a_0 + b_0 + c_0

i += 1
a_t = my_pow(10, i, 10*a)/a
b_t = my_pow(10, i, 10*b)/b
c_t = my_pow(10, i, 10*c)/c
s_t = a_t+b_t+c_t

if s_t>=10:
return 1
elif s_t == 8:
return (s_0+if_8(a, b, c, i))/10
elif s_t == 9:
return (s_0+if_9(a, b, c, i))/10
else :
return 0

def if_19(a, b, c, i):
a_0 = my_pow(10, i, 10 * a) / a
b_0 = my_pow(10, i, 10 * b) / b
c_0 = my_pow(10, i, 10 * c) / c
s_0 = a_0 + b_0 + c_0

i += 1
a_t = my_pow(10, i, 10*a)/a
b_t = my_pow(10, i, 10*b)/b
c_t = my_pow(10, i, 10*c)/c
s_t = a_t+b_t+c_t

if s_t>=10:
return 2
elif s_t == 8:
return (s_0+if_8(a, b, c, i))/10
elif s_t == 9:
return (s_0+if_9(a, b, c, i))/10
else :
return 1

def if_18(a, b, c, i):
a_0 = my_pow(10, i, 10 * a) / a
b_0 = my_pow(10, i, 10 * b) / b
c_0 = my_pow(10, i, 10 * c) / c
s_0 = a_0 + b_0 + c_0

i += 1
a_t = my_pow(10, i, 10*a)/a
b_t = my_pow(10, i, 10*b)/b
c_t = my_pow(10, i, 10*c)/c
s_t = a_t+b_t+c_t

if s_t>=20:
return 2
elif s_t == 18:
return (s_0+if_18(a, b, c, i))/10
elif s_t == 19:
return (s_0+if_19(a, b, c, i))/10
elif s_t == 8:
return (s_0+if_8(a, b, c, i))/10
elif s_t == 9:
return (s_0+if_9(a, b, c, i))/10
else :
return 1

if __name__ == '__main__':
a = 123
b = 234
c = 456

print '1/a的10000位小数点为'+str(my_pow(10, 10000, 10*a)/a)
print '1/b的10000位小数点为'+str(my_pow(10, 10000, 10*b)/b)
print '1/c的10000位小数点为'+str(my_pow(10, 10000, 10*c)/c)
print '\n'
print '1/a的10001位小数点为'+str(my_pow(10, 10001, 10*a)/a)
print '1/b的10001位小数点为'+str(my_pow(10, 10001, 10*b)/b)
print '1/c的10001位小数点为'+str(my_pow(10, 10001, 10*c)/c)
print '\n'
print '1/a的10002位小数点为'+str(my_pow(10, 10002, 10*a)/a)
print '1/b的10002位小数点为'+str(my_pow(10, 10002, 10*b)/b)
print '1/c的10002位小数点为'+str(my_pow(10, 10002, 10*c)/c)
print '\n'
print '1/a+1/b+1/c的10000位小数点为'+str(add(a, b, c))

bbp算法

ACM题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=6217

BBP公式提出了求$\pi$小数点后任意位的算法,只是求的是16进制表示的$\pi$

那么依照上面的我们求$\left \lfloor16^n\pi\right\rfloor\;mod\;16$就可以了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def my_pow(a, b, c):
a = a % c
res = 1

while b != 0:
if b & 1:
res = res * a
b >>= 1
a = (a ** 2) % c

return res % c


def S(d, j):
sum = 0

for k in range(d):
sum+=my_pow(16, d-k, 16*(8*k+j))/(8*k+j)

return sum

def accuracy(d, j, a):
sum = 0

for k in range(d+a):
sum += my_pow(16, d - k+a , (16**a) * (8 * k + j)) / (8 * k + j)

return sum

def cal(d, a):
sum = (4*S(d, 1) - 2*S(d, 4)-S(d, 5)-S(d, 6))
aa = (4*accuracy(d, 1, a) - 2*accuracy(d, 4, a)-accuracy(d, 5, a)-accuracy(d, 6, a))


if aa<0:
sum -= (-aa)/(16**a)+1
else:
sum += aa/(16**a)

return sum % 16

if __name__ == '__main__':
a = 5
i = 1
print '第%d位为:' % (i)+hex(cal(i, a))[2:].strip('L')
i = 11
print '第%d位为:' % (i) + hex(cal(i, a))[2:].strip('L')
i = 111
print '第%d位为:' % (i) + hex(cal(i, a))[2:].strip('L')
i = 1111
print '第%d位为:' % (i) + hex(cal(i, a))[2:].strip('L')
i = 11111
print '第%d位为:' % (i) + hex(cal(i, a))[2:].strip('L')

虽然可以使用上面的递归计算来实现精确进位,但通常只考虑后面几位就能有较高的准确率了,代码中的a即表示精确的位数

Tags: 随笔

1000000