GCD序列

定义序列$g(n)$，初值为$g(4)=13$，之后有 $$g(n)=g(n-1)+\gcd(n,g(n-1)).$$ 从$n=4$到$n=20$，序列中项依次是$13, 14, 16, 17, 18, 27, 28, 29, 30, 31, 32, 33, 34, 51, 54, 55, 60,\ldots$。

可以验证$g(1000)=2524, g(1000000)=2624152$。

求$g(10^{15})$。

本题难度:

解答

若$\gcd(n,g(n-1))=1$，那么$g(n)=g(n-1)+1$，观察所给的例子可以注意到这是该数列中的普遍现象，尤其是可以注意到$g(1000)$与$1000$的数量级相当，$g(1000000)$同理。

所以本题的关键是找到产生跳跃的位置，即满足$\gcd(n,g(n-1))\neq1$的n的位置。

设当前位置是n，下一个跳跃点是m，亦即m是最小的满足m>n，且$\gcd(m,g(m-1))\neq1$的数，特别地，由此定义可得$g(m-1)=g(n)+m-n-1$。

记$k=m-n-1, g=g(n)$，那么需要找到最小的k，使得$\gcd(n+k+1,g+k)>1$，亦即$\gcd(n-g+1,g+k)>1$。

枚举$n-g+1$的素因数，对每个素因数p，找出最小的$k_p$，使得$g+k_p$是p的倍数，再取$\min_pk_p$即得k。

重复这一过程直到n达到或超过$10^{15}$，最终结果是$2744233049300770$。

import sympy,math

target=10**15

n,g=4,13
while True:
    d=min((g+p-1)//p*p for p in sympy.primefactors(g-n-1))
    m=n+d-g+1
    if m>target:
        g+=target-n
        break
    g=d+math.gcd(m,d)
    n=m

print(g)