題意：$T$ 組數據，每組給定 $n,a$，求滿足下列條件的項鏈數量：

• 有 $n$ 個珠子。
• 每個珠子上有三個 $[1,a]\cap \mathbb{Z}$ 的數，且三個數 $\gcd$ 為 $1$。
• 相鄰兩個珠子不同。

珠子旋轉、翻轉同構，項鏈旋轉同構。對 $10^9+7$ 取模。

$n\le 10^{14},a\le 10^7$

有毒……

顯然是道縫合怪題，先考慮求有多少個不同的珠子。

設 $f(k)$ 表示 $k$ 個 $[1,a]$ 的數 $\gcd$ 為 $1$ 的方案數，用 Polya 定理數一下可知

$$ans=\frac{f(3)+3f(2)+2f(1)}{6}$$

然后就是個簡單的反演

$$f(k)=\sum _{i=1}^a\mu (i){?\frac{a}{i}?}^k$$

整除分塊算就可以了（其實直接暴力也可以？）。

然后第二步，數項鏈。設第一步得到的珠子種數為 $k$，眾所周知，$F(n)$ 為 $n$ 個珠子不循環同構的方案，答案就是 $F$ 和 $\phi$ 的狄利克雷卷積除以 $n$。

考慮 $F$ 怎么算。

考慮一個 $n$ 個珠子的合法方案，我們刪掉編號為 $1$ 的點。如果兩邊的點不同，就對應了一種 $F(n?1)$ 的方案，乘上 $1$ 號珠子的方案 $(k?2)$。如果相同，任意刪掉一個后就對應一種 $F(n?2)$ 的方案，乘上 $1$ 號珠子的方案 $(k?1)$。所以

$$F(n)=(k?2)F(n?1)+(k?1)F(n?2)$$

注意邊界是 $F(1)=0,F(2)=k(k?1)$,倒退出來 $F(0)=m$。

列出生成函數：

$$F=(k?2)xF+(k?1)x^{2}F+k?k(k?2)x$$

這里減 $k(k?2)x$ 是為了平衡 $(k?2)xF$ 對一次項的影響。

即

$$F=\frac{?m(m?2)x+m}{(1+x)(1?(m?1)x)}$$

直接待定系數拆掉

$$F=\frac{a}{1+x}+\frac{b}{1?(k?1)x}$$

列出方程：

$$a?a(k?1)x+b+bx=?k(k?2)x+k$$

解得 $a=k?1,b=1$

所以

$$F=\frac{k?1}{1+x}+\frac{1}{1?(k?1)x}$$

$$F(n)=(k?1{)}^n+(?1{)}^n(k?1)$$

然后就可以算了。

這里有點卡，就不用 Polya 模板那個垃圾做法，而是把因子預處理出來 dfs，可以做到 $\mathrm{O}(d(n))$ 的優秀復雜度。

然后因為 $n$ 很大，最后除 $n$ 的時候又會有喜聞樂見的除 $0$ 問題。所以要對 $(10^9+7{)}^2$ 取模，最后再判一下。

總復雜度 $\mathrm{O}(n+T(\sqrt{a}+d(n)\log ))$

#include <iostream> #include <cstdio> #include <cstring> #include <cctype> #define int long long using namespace std; const int mod=1e9+7,MOD=mod*mod; inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;} inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;} inline int mul(int a,int b){return (a*b-(int)((long double)a/MOD*b+0.5)*MOD+MOD)%MOD;} inline int qpow(int a,int p) {int ans=1;while (p){if (p&1) ans=mul(ans,a);a=mul(a,a),p>>=1;}return ans; } int INV6=833333345000000041ll; int n,a; const int N=1e7,MAXN=N+5; bool np[MAXN]; signed pl[MAXN],cnt,mu[MAXN]; void init() {np[1]=1,mu[1]=1;for (int i=2;i<=N;i++){if (!np[i]) mu[pl[++cnt]=i]=-1;for (int j=1,x;(x=i*pl[j])<=N;j++){np[x]=1;if (i%pl[j]==0) break;mu[x]=-mu[i];}}for (int i=2;i<=N;i++) mu[i]+=mu[i-1]; } int calc() {int ans=0;for (int l=1,r;l<=a;l=r+1){r=a/(a/l);int t=mu[r]-mu[l-1];(t<0)&&(t+=MOD);ans=add(ans,mul(t,add(mul(a/l,mul(a/l,a/l)),3*mul(a/l,a/l)%MOD)));}return ans; } int k,ans; int lis[20005],siz[20005],tot; void dfs(int pos,int phi,int res) {if (pos>tot){int t=qpow(k,res);if (res&1) t=dec(t,k);else t=add(t,k);ans=add(ans,mul(phi,t));return;}int t=qpow(lis[pos],siz[pos]);dfs(pos+1,phi,res*t);t/=lis[pos];dfs(pos+1,phi*=(lis[pos]-1),res*t);t/=lis[pos];if (!t) return;for (;t;dfs(pos+1,phi*=lis[pos],res*t),t/=lis[pos]); } signed main() {init();int T;cin>>T;while (T--){cin>>n>>a;k=dec(mul(calc()+2,INV6),1);tot=ans=0;int t=n;for (int i=1;i<=cnt;i++)if (t%pl[i]==0){lis[++tot]=pl[i],siz[tot]=0;while (t%pl[i]==0) t/=pl[i],++siz[tot];}if (t>1) lis[++tot]=t,siz[tot]=1;dfs(1,1,1);if (n%mod==0) n/=mod,ans/=mod;ans=mul(ans,qpow(n,mod-2));cout<<ans%mod<<'\n';}return 0; }

總結

以上是生活随笔為你收集整理的【SDOI2013】项链【莫比乌斯反演】【Polya定理】【递推式求通项】【数论】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。