1.題目

求

\[\sum_{i=1}^n \sum_{j=1}^m \gcd(F_i,F_j) \]

其中 \(F_k\) 表示斐波那契數列的第 \(k\) 項，對 \(10^9 + 7\) 取模。

多組資料。

2.題解

莫比烏斯反演板子題，但是太菜了，多做了一次反演，複雜度變為 \(tn\sqrt{n}\) 。實際是 \(t\sqrt{n}\)

直接推式子吧。

首先需要知道性質，\(\gcd(F_i,f_j)=F_{\gcd(i,j)}\)

這個性質是一道板子題，為洛谷上的斐波那契公約數，證明簡單，本文略過。

\[Ans=\sum_{i=1}^n\sum_{j=1}^m \gcd(F_i,F_j) \]\[=\sum_{i=1}^n\sum_{j=1}^mF_{\gcd(i,j)} \]

我們發現 \(\gcd(i,j)\)

只有可能在 \(1\sim\min(n,m)\) 於是我們可以考慮去列舉這個 \(\gcd(i,j)\) ，然後乘上所對應的值，這樣既為答案。

也就是說，寫成這樣(假設 \(n \leq m\))：

\(f(k)\) 表示的是公約數為 \(k\) 的數量。

\[Ans=\sum_{k=1}^n F(k)f(k) \]

問題關鍵在於求 \(f(k)\) 。

\[Ans =\sum_{k=1}^n F(k) \sum_{i=1}^n\sum_{j=1}^m [(i,j)=k] \]

容易發現這就是一個嵌入式反演的變形，那麼直接上莫比烏斯反演。

\[=\sum_{k=1}^n F(k) \sum_{k|i}^n\sum_{k|j}^m[(i,j)=k] \]\[=\sum_{k=1}^n F(k) \sum_{i=1}^{n/k}\sum_{j=1}^{m/k}[(ik,jk)=k] \]

發現可以將 \(k\)

約掉，也就是：

\[=\sum_{k=1}^n F(k) \sum_{i=1}^{n/k}\sum_{j=1}^{m/k}[(i,j)=1] \]

變為經典反演形式，開始進行反演。

\[Ans=\sum_{k=1}^n F(k) \sum_{i=1}^{n/k}\sum_{j=1}^{m/k} \sum_{d|(i,j)} \mu(d) \]\[Ans=\sum_{k=1}^n F(k) \sum_{i=1}^{n/k}\sum_{j=1}^{m/k} \sum_{{d|i,}{d|j}} \mu(d) \]

然後改變列舉變數。

\[Ans=\sum_{k=1}^n F(k) \sum_{d=1}^{n/k} \mu(d) \sum_{d|(n/k)}\sum_{d|(m/k)}1 \]

也就是:

\[Ans=\sum_{k=1}^n F(k) \sum_{d=1}^{n/k} \mu(d) \lfloor \dfrac{n}{k} \rfloor\lfloor \dfrac{m}{k} \rfloor \]

然後交換求和順序，以及內部改為列舉因數，最外層列舉 \(d\) ，就有：

\[Ans=\sum_{d=1}^n \lfloor \dfrac{n}{d} \rfloor \lfloor \dfrac{m}{d} \rfloor \sum_{k|d} F_k \mu(\dfrac{k}{d}) \]

然後就預處理字首和，然後套路整除分塊回答。

時間複雜度為 \(t\sqrt{n}+nlogn\)

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6+99 ,mod = 1e9+7;
int e[N+4],p[N+4],mu[N+4],tn;
void mobius(int n){
 e[1]=1;mu[1]=1;
 for(int i=2;i<=n;i++){
  if(!e[i]){mu[i]=-1;p[++tn]=i;}
   for(int j=1;j<=tn;j++){
	if(i*p[j]>n) break;
	mu[p[j]*i]=(i%p[j]==0 ? 0 :-mu[i]); 
	e[p[j]*i]=1;
	if(i%p[j]==0) break;
  }
 }
}
int g[N+4],T,n,m,fib[N+4];
signed main(){
 freopen("fibonacci.in","r",stdin);
 freopen("fibonacci.out","w",stdout);	
 ios::sync_with_stdio(0);cin.tie(0),cout.tie(0);
 mobius(N);
 cin>>T;
 fib[1]=1,fib[2]=1;
 for(int i=3;i<=N;i++){
  fib[i]=fib[i-1]+fib[i-2];
  fib[i]%=mod;
 }		
 for(int i=1;i<=N;i++){
  for(int j=i;j<=N;j+=i){
   g[j]=(g[j]+fib[i]*mu[j/i]%mod+mod)%mod;	
  }	
 }
 for(int i=1;i<=N;i++)
  g[i]=(g[i]+g[i-1])%mod;	
 while(T--){
  cin>>n>>m;
  int ans=0;
  for(int l=1,r=0;l<=min(n,m);l=r+1){
   r=min(n/(n/l),m/(m/l));	
   ans=(ans+(n/l)*(m/l)%mod*(g[r]-g[l-1])%mod+mod)%mod;
   	
  }	
  cout<<ans<<"\n";
 }
 return 0;	
}