欧拉计划第80题：平方根数字展开

众所周知，如果一个自然数的平方根不是整数，那么就一定是无理数，也即无限不循环小数。

2的平方根为1.41421356237309504880，它的前一百个数字（包括小数点前的1）的和是475。

考虑前一百个自然数的平方根，求其中所有无理数的前一百个数字的总和。

---

解：

平方根展开式

网上搜到如下的的展开式。

试着编写程序，由于不知道迭代多少次，可以精确到100位小数，就循环了2000次，为了计算分数的精确小数表示，用了decimal。设置100位小数可能会有舍入误差，这里设置102的精度比较安全。

import decimal
import math

decimal.getcontext().prec = 102


def root(x, repeat):
    a, b = 1, 1
    for _ in range(repeat):
        a, b = (x - 1) * b, a + 2 * b
    return a + b, b


def sum_digits(frac, digits=100):
    a, b = frac
    a_div_b = decimal.Decimal(a) / decimal.Decimal(b)
    digits = str(a_div_b).replace('.', '')[:digits]
    return sum(map(int, digits))


assert sum_digits(root(2, 1000)) == 475

total = 0
for n in range(1, 101):
    sq = round(math.sqrt(n))
    if sq * sq == n:
        continue
    total += sum_digits(root(n, 2000))

print(total)
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

另一种展开式

还有一种展开式，收敛速度更快一些。我只循环了200次。

import decimal
import math

decimal.getcontext().prec = 102


def root(x, repeat):
    a = int(math.sqrt(x))
    r = x - a * a
    m, n = 1, 1
    for _ in range(repeat):
        m, n = n * r, 2 * a * n + m
    return a * n + m, n


def sum_digits(frac, digits=100):
    a, b = frac
    a_div_b = decimal.Decimal(a) / decimal.Decimal(b)
    digits = str(a_div_b).replace('.', '')[:digits]
    return sum(map(int, digits))


assert sum_digits(root(2, 1000)) == 475

total = 0
for n in range(1, 101):
    sq = round(math.sqrt(n))
    if sq * sq == n:
        continue
    total += sum_digits(root(n, 200))

print(total)
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

decimal里支持sqrt

实际上decimal里直接可以进行sqrt()运算。

import decimal
import math

decimal.getcontext().prec = 102


def sum_digits(n):
    sq = str(decimal.Decimal(n).sqrt())
    digits = sq.replace('.', '')[:100]
    return sum(map(int, digits))


assert sum_digits(2) == 475

total = 0
for n in range(1, 101):
    sq = round(math.sqrt(n))
    if sq * sq == n:
        continue
    total += sum_digits(n)

print(total)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22