Project Euler 817

题目

Digits in Squares

Define \(m = M(n, d)\) to be the smallest positive integer such that when \(m^2\) is written in base \(n\) it includes the base \(n\) digit \(d\). For example, \(M(10,7) = 24\) because if all the squares are written out in base \(10\) the first time the digit \(7\) occurs is in \(24^2 = 576\). \(M(11,10) = 19\) as \(19^2 = 361=2A9_{11}\).

Find \(\displaystyle \sum_{d = 1}^{10^5}M(p, p - d)\) where \(p = 10^9 + 7\).

解决方案

令\(p = 10^9 + 7,D=10^5\)，把任意正整数 \(x\) 写成 \(p\) 进制：\(x = x_0 + x_1 p + x_2 p^2 + \cdots, 0\le x_i < p\)，其中 \(x_0\) 是个位，\(x_1\) 是 \(p\) 位（相当于十位），等等。

注意 \(p-d\) 非常接近 \(p\)，即是一个接近进制上界的大数字。

我们要找的数字 \(p-d\) 满足：\(p-d \ge p-D\)，这是一个非常大的数字。如果它第一次出现的位置在 \(p^2\) 位或更高位，则 \(m^2\) 至少要达到数量级 \(p^2\cdot(p-D) \sim p^3\)，从而\(m \ge \sqrt{(p-d)p^2} \approx p\sqrt{p}\)

这对应的 \(m\) 量级约 \(3\times 10^{13}\)，远远大于我们即将构造出的 \(m\)（量级大约 \(10^9\) 到 \(10^{10}\)）。因此真正的最小值一定会让 \(p-d\) 尽可能出现在低位，只需要考虑个位或 \(p\) 位。

若 \(p-d\) 出现在 \(m^2\) 的个位，则：\(m^2 \equiv p-d \pmod p\)这是模 \(p\) 的二次同余是否可解问题。

对素数 \(p\)，有 Euler 准则：

\[ \left(\frac{a}{p}\right)= \begin{cases} 1 & \text{if } a \in (\mathbb{F}_p^\times)^2,\\ -1 & \text{if } a \notin (\mathbb{F}_p^\times)^2, \end{cases} \Longleftrightarrow a^{\frac{p-1}{2}} \equiv \left(\frac{a}{p}\right)\pmod p \]

其中：\((\mathbb{F}_p^\times)^2={x^2 \bmod p : x\in\mathbb{F}_p^\times}\) 表示模 \(p\) 的非零平方剩余集合

因此对 \(a = p-d\)：若 \(a^{\frac{p-1}{2}}\equiv 1\pmod p\)，则 \(a\) 是二次剩余，有解，否则无解（个位不可能出现该数字）。

本题\(p\) 满足 \(p \equiv 3 \pmod 4\)，此时若 \(a\) 是二次剩余，则\(\sqrt a \equiv a^{\frac{p+1}{4}} \pmod p\)。得到一个根 \(r\) 后，另一个根是 \(p-r\)。最小正整数解显然取：\(M(p, p-d)=\min(r,p-r)\)，因为这意味着直接让 \(m< p/2\)，并且 \(m^2\) 的个位已经出现目标数字，无法再更小。

当 \(a=p-d\) 是二次非剩余时，\(m^2 \bmod p\) 永远不等于 \(a\)，所以它不可能在个位第一次出现。下一优位置就是 p 位。设 \(m^2\) 写成：\(m^2 = y p^2 + a p + r, 0\le r < p\)

则其 \(p\) 进制下第二位数字就是 \(a\)（也就是 \(p-d\)）。等价于：\(y p^2 + ap \le m^2 < y p^2 + (a+1)p\)

记区间上下界：\(L_y = y p^2 + ap, U_y = L_y + p\)，那么我们要找的就是最小 \(m\)，使得存在某个 \(y\ge 0\)：\(m^2 \in [L_y, U_y)\)，并且希望 \(m\) 最小，因此自然从 \(y=0,1,2,\dots\) 依次尝试。

由于平方函数单调递增。对于区间左端点 \(L_y\)，令\(m_0 = \left\lceil \sqrt{L_y}\right\rceil\)，则：\(m_0\) 是使得 \(m^2 \ge L_y\) 的最小 \(m\)。若 \(m_0^2 < U_y\)，则区间内存在平方（且就是最早的那个平方），否则此 \(y\) 不行，继续增大 \(y\)。由此完成\(p-d\)出现在第二位的情况。

代码

from tools import *
P = 10 ** 9 + 7
K = 10 ** 5
def solve(p, k):
    p2 = p * p
    ans = 0

    euler_exp = (p - 1) // 2
    sqrt_exp = (p + 1) // 4

    for d in range(1, k + 1):
        a = p - d
        if pow(a, euler_exp, p) == 1:
            r = pow(a, sqrt_exp, p)
            if r > p - r:
                r = p - r
            ans += r
        else:
            y = 0
            while True:
                L = y * p2 + a * p
                m = int_sqrt(L)
                if m * m < L:
                    m += 1
                if m * m < L + p:
                    ans += m
                    break
                y += 1

    return ans

if __name__ == "__main__":
    print(solve(P, K))