0399 无平方因子
* * * *
拉格朗日计划
* * * *
无平方因子

前15个斐波那契数是1,1,2,3,5,8,13,21,34,55,89,144,233,377,610。

可以看出其中8和144都不是无平方因子数,因此前13个无平方因子斐波那契数是1,1,2,3,5,13,21,34,55,89,233,377和610。

第200个无平方因子斐波那契数是971183874599339129547649988289594072811608739584170445。

用科学计数法可以写作9.7e53,同时这个数的最后十六位数字是1608739584170445。

求第$10^8$个无平方因子斐波那契数,并将结果写作如下形式:

该数的最后十六位数字,该数的科学计数法表示(四舍五入到小数点后一位),两者间用半角逗号隔开。

例如对于第200个无平方因子斐波那契数,按上面的格式应写作1608739584170445,9.7e53

注意: 本题中假定对于任意素数p,第一个能被p整除的斐波那契数不能同时被$p^2$整除(这是Wall猜想的一部分)。目前此假设对于不超过$3\times10^{15}$的素数已获得到验证,其它情况则未知。 若该假定最终被证伪,则本题所预期的答案并不一定是第$10^8$个无平方因子斐波那契数,而只能作为其下界。

本题难度:



解答

斐波那契数列的递推公式是线性递推式,因此容易验证其模p具有周期性。

记$F_n$为第n个斐波那契数,本题还用到这样一个结论(也就是题面中提及的Wall猜想的一部分):若n是第一个使得$F_n$能被素数p整除的数,那么$F_{np}$是第一个能被$p^2$整除的数,且这样的数以$pn$为周期出现。

接下来的策略就很直接,先筛出充分多的素数,对每个素数p,找出首个能被p整除的$F_n$(事实上n一定是$p$或$p\pm1$的约数,不过对本题而言可以直接枚举计算),再用筛法排除np的倍数(此外由于本题的数据量很大,因而需要使用bytearray实现)。

最终结果是$1508395636674243,6.5e27330467$。

注:为方便书写,代码中用到walrus算符、fstring等新特性,因此以下源码为Python 3。此外,本题数据量较大,需要数分钟运行,但并未输出进度。

from math import log10,sqrt,modf

S=10**8
T=3*S//2
mod=10**16

def g(p):
    f1=f2=1
    n=2
    while f2 and n*p<=T:
        n+=1
        f1,f2=f2,(f1+f2)%p
    return n if n*p<=T else 0

mul=lambda A,B:[[(sum(A[i][k]*B[k][j] for k in range(len(A))))%mod for i in range(len(A))] for j in range(len(A))]

def f(a): 
    b=format(a,'b')
    X=[[1,0],[0,1]]
    L=[[1,1],[1,2]]       
    for d in b[-2::-1]:
        if d=='1':
            X=mul(X,L)
        L=mul(L,L)
    return X[1][1] if b[-1] == '1' else X[1][0]

target=10**7
d=[0]*target
n=2
while n < target:
    for i in range(n*n,target,n):
        d[i]=d[i]+1
    n=n+1
    while n < target and d[n]>0:
        n=n+1

sieve=bytearray([1]*T)
sieve[0]=0

for p in range(2,target):
    if d[p]==0 and (a:=p*g(p))>0:
        for i in range(a,T,a):
            sieve[i]=0

t=[i for i,j in enumerate(sieve) if j==1][S-1]
a,b=modf(t*log10((1+sqrt(5))/2)-log10(sqrt(5)))
print(f"{f(t)},{10**a:.1f}e{int(b)}")