Baom and Fibonacci 整除分块 + 杜教筛

Baom and Fibonacci 整除分块 + 杜教筛

题意

定义\(fib_i\)为斐波那契数列

计算

\[f(n) = \sum_{i=1}^n\sum_{j=1}^ngcd(fib(i),fib(j)) (mod \ 998244353) \]

\[1\leq 10^{10} \]

分析

\[f(n) = \sum \sum fib_{(i,j)} \\ = \sum_{d=1}^n fib(d)\sum \sum [(i,j)=d]\\ = \sum_{d=1}^n fib(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[(i,j) = 1]\\ = \sum_{d=1}^n fib(d)(2 \sum \phi(\lfloor\frac{n}{d}\rfloor) - 1) \]

于是矩阵快速幂计算fib + 杜教筛计算\(\sum \phi\) 即可

复杂度正确性在于杜教筛记忆化搜索

代码

struct mat{
    ll a[2][2];
};
 
 
mat mat_mul(mat x,mat y){
    mat res;
    memset(res.a,0,sizeof(res.a));
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        for(int k=0;k<2;k++)
        res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j] % MOD)%MOD;
    return res;
}
ll mat_pow(ll n){
    mat c,res;
    c.a[0][0]=c.a[0][1]=c.a[1][0]=1;
    c.a[1][1]=0;
    memset(res.a,0,sizeof(res.a));
    for(int i=0;i<2;i++) res.a[i][i]=1;
    while(n)
    {
        if(n&1) res=mat_mul(res,c);
        c=mat_mul(c,c);
        n >>= 1;
    }
    return res.a[0][1];
}
const int maxx=4700000;
const int mod=998244353;
map<ll,ll> P;
bool isP[maxx];
int prime[maxx];
int cnt;
ll phi[maxx];
ll inv=499122177;
void init(){
    phi[1]=1;
    for(int i=2;i<maxx;i++)
    {
        if(!isP[i]){prime[cnt++]=i;phi[i]=i-1;}
        for(int j=0;j<cnt&&(ll)i*prime[j]<maxx;j++)
        {
            isP[i*prime[j]]=true;
            if(i%prime[j])
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            else
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
        }
    }
    for(int i=1;i<maxx;i++)
        phi[i]=(phi[i]+phi[i-1])%mod;
}
ll work(ll x){
    if(x<maxx)return phi[x];
    if(P[x])return P[x];
    ll ans=((x+1)%mod*(x%mod)%mod)*inv%mod;
    for(ll i=2,last;i<=x;i=last+1)
    {
        last=x/(x/i);
        ans=(ans-(last-i+1)*work(x/i)%mod)%mod;
    }
    ans=(ans+mod)%mod;
    P[x]=ans;
    return ans;
}
ll getf(ll i,ll j){
	return (mat_pow(j + 2) - mat_pow(i + 1) + MOD) % MOD;
} 
 
int main(){
	init();
	ll n = rd();
	ll ans = 0;
	for(ll i = 1,j;i <= n;i = j + 1) {
		j = n / (n / i);
		ans += (2 * work(n / i) % MOD - 1 + MOD) % MOD * getf(i,j) % MOD;
		ans %= MOD;
		//cout << i << ' ' << j << ' ' << ans << '\n';
	}
	printf("%lld",ans);
}