欧拉计划第932题：2025-数

问题描述

对于年份2025，有下面等式成立。

$$ 2025 = (20 + 25)^2 $$

给定两个正整数 a 和 b，如果将它们拼接在一起，结果满足：

$$ a ※ b = (a + b)^2 $$

则称这个拼接数 a※b 为一个 2025-数。

据此定义，3025 和 81 都是 2025-数。

注意：9801 不是 一个 2025-数，因为 98 和 1 拼接后得到的是 981，而不是 9801。

设 T(n) 表示所有位数不超过 n 的 2025-数之和。

已知：T(4) = 5131

求：T(16)

第一步，根据题意，暴力求解

给定一个数n，拆分为a和b，判断是否满足(a+b)^2==n，遍历所有情况。

def split_number(n: int):
    s = str(n)
    result = []

    # 拆分位置从第1位到第 n-1 位
    for i in range(1, len(s)):
        a, b = s[:i], s[i:]
        if b[0] != '0': # b不能有前导零
            result.append((int(a), int(b)))
    return result


def is_2025_number(n):
    for a, b in split_number(n):
        if (a+b)**2 == n:
            print(a, b)
            return True
    return False


def T(n):
    return sum(i for i in range(11, 10**n) if is_2025_number(i))


assert(T(4) == 5131)

print(T(7))

这个算法效率非常差，计算T(7)耗时25秒，T(8)需要7分钟。

把8位数之内的所有2025-数打印出来，便于以后寻找规律：

8,1
20,25
30,25
88,209
494,209
494,1729
744,1984
2450,2500
2550,2500
5288,1984
6048,1729

第二步，改进拆分函数

原来的整数拆分为a和b的算法使用字符串运算，如果全部采用整数计算，可以稍快一点。T(7)从25秒提升至15秒，但并没有本质的提升。

def split_number(n: int):
    result = []
    power = 10
    for _ in range(1, len(str(n))):
        a, b = divmod(n, power)
        if b >= power // 10: # b不能有前导零
            result.append((a, b))
        power *= 10
    return result

第三步，发现拆分的规律

可以发现一个规律，如果n的位数为偶数，则在中间分隔，例如：5288,1984。如果是奇数，则在中间偏左的位置分隔，例如：494,1729。

改进拆分的函数，T(7)可提升至4秒，但计算T(16)仍不可能。

def split_number(n: int):
    n_len = len(str(n))          
    if n_len % 2 == 0: # 偶数情况，中间拆分
        power = 10 ** (n_len // 2)   
    else:  # 左边比右边的位数少1位
        power = 10 ** (n_len // 2 + 1)   

    a, b = divmod(n, power)
    if b >= power // 10:
        return a, b
    else: 
        return 0, 0


def is_2025_number(n):
    a, b = split_number(n)
    return  (a+b)**2 == n


def T(n):
    return sum(i for i in range(11, 10**n) if is_2025_number(i))


assert(T(4) == 5131)
print(T(7))

第四步，根据a,b寻找n

枚举位数分别为 p、q 的所有整数 a、b，检查ab拼接后的数是否是2025-数。

find_2025_number(p,q) 负责找出所有p位和q位数组成的 (a,b) 对； T(m) 遍历可能的 (p,q) 位数组合，汇总这些结果。这里q=p，或q=p+1。

T(7)计算只需0.7秒。T(8)需7秒。T(9)需78秒。 但要计算T(16)仍不可能。

'''
寻找所有p位数和q位数组成的2025-数。
'''
def find_2025_number(p, q):
    results = []
    ten_p, ten_q = 10 ** p, 10 ** q

    for a in range(ten_p // 10, ten_p):
        for b in range(ten_q // 10, ten_q):
            if (a + b) ** 2 == a * ten_q + b:
                results.append((a+b)**2)
                print(f"{a},{b}")

    return results


def T(m):
    n = m // 2

    total = 0
    for p in range(1, n + 1):
        # 直接处理 q = p
        vals = find_2025_number(p, p)
        total += sum(vals)

        # 处理 q = p + 1的情况
        q = p + 1
        if q + n <= m:
            vals = find_2025_number(p, q)
            total += sum(vals)

    return total


assert(T(4) == 5131)
print(T(7))

第五步，进行数学推导

根据题意，a, b都是正整数，假设q是b的位数，则有：

$$ a\cdot 10^q + b = (a+b)^2 $$$$ a\cdot 10^q + b = a^2 + 2ab + b^2 $$

整理为关于b的一元二次方程：

$$ b^2 + (2a-1)b + (a^2 - a\cdot 10^q)=0 $$

解这个方程：

$$ b=\frac{1-2a\pm\sqrt{4a(10^q-1)+1}}{2} $$

因为 b>0，所以舍弃负数根。

$$ b=\frac{1-2a+\sqrt{4a(10^q-1)+1}}{2} $$

因为b为正整数，所以根号内的值必为正整数，而且是完全平方数。

令：

$$ D^2 = 4a(10^q-1)+1 $$

这里的D为正整数。

$$ b=\frac{1-2a+D}{2} $$$$ a + b=\frac{D+1}{2} $$

因为(D+1)被2整除，所以D为奇数。

对于n位的2025-数，

$$ (a+b)^2 < 10^n $$

即

$$ (\frac{D+1}{2})^2 < 10^n $$

求T(n)时，D的上界可以由下面不等式确定。

$$ D<2\sqrt{10^n}-1 $$

import math

def T(n):
    total = 0
    Dmax = int(2 * math.sqrt(10**n)) - 1
    for D in range(3, Dmax + 1, 2):
        dd1 = D * D - 1
        if dd1 % 36:
            continue

        a_plus_b = (D + 1) // 2
        # b的位数肯定不大于(a+b)的位数
        qmax = math.ceil(math.log10(a_plus_b))
        for q in range(1, qmax + 1):
            ten_q = 10 ** q
            t = 4 * (ten_q - 1)
            if dd1 % t == 0:
                a = dd1 // t
                b = a_plus_b - a
                if b < ten_q // 10:
                    continue
                total += a_plus_b ** 2
                print(f'{a},{b}  q={q}  D={D}  a+b={a_plus_b}')
    return total

assert(T(4) == 5131)
assert(T(5) == 93340)
print(T(16))

这个程序比以前快了无数倍，大约30秒可以计算出T(16)的值。

其实程序里的q就等于qmax，可以去除关于q的循环，还可以提高速度，但我没有找到 b的位数==(a+b)的位数 证明方法。

第六步，换一种推导方式

还可以整理为关于 a 的一元二次方程：

$$a^2 + (2b - 10^q)a + (b^2 - b) = 0$$

解出a的两个根：

$$ a = \frac{10^q - 2b \pm \Delta}{2} $$

这里的Δ满足：

$$ \Delta^2 = 10^{2q} - 4b(10^q - 1) $$

方程如果有解，则平方根号内的数必定大于等于0，即：

$$ 10^{2q} - 4b(10^q - 1) \ge 0 $$$$ b \le \frac{10^{2q}}{4(10^q - 1)} = \frac{10^q}{4}\cdot\frac{10^q}{10^q - 1}. $$

最后的源代码：

import math
'''
已知b是q位数，寻找可能的a，构成2025-数。
'''
def find_2025_number(q):
    results = []
    ten_q = 10 ** q

    bmax = (int)((ten_q / 4) * (ten_q / (ten_q - 1)))
    for b in range(ten_q // 10, bmax + 1):
        delta2 = ten_q ** 2 - 4 * b * (ten_q - 1)
        delta = math.isqrt(delta2)
        if delta * delta == delta2:
            a1 = (ten_q - 2*b - delta) // 2
            if a1 > 0 and (a1 + b) ** 2 == a1 * ten_q + b:
                    results.append((a1+b)**2)
                    print(f"{a1},{b}  delta={delta}")
            a2 = (ten_q - 2*b + delta) // 2
            if (a2 + b) ** 2 == a2 * ten_q + b:
                    results.append((a2+b)**2)
                    print(f"{a2},{b}  delta={delta}")

    return results


def T(n):
    total = 0
    for q in range(1, (n+1)//2+1): # 要注意n是奇数的情况
        vals = find_2025_number(q)
        total += sum([x for x in vals if len(str(x))<=n]) #超出位数的需舍弃掉
    return total

assert(T(4) == 5131)
assert(T(5) == 93340)
print(T(16))

3秒可计算出T(16)。