題面
https://www.luogu.com.cn/problem/P4464
題解
前置知識
-
Pollard-rho https://www.luogu.com.cn/problemnew/solution/P4718
-
Miller-rabin https://www.cnblogs.com/xh092113/p/12288424.html
-
莫比烏斯反演：《具體數學：計算機科學基礎》4.9節
-
伯努利數：《具體數學：計算機科學基礎》6.5節
本題要求\({\sum_{i=1}^{n}}(i,n)^x[i,n]^y\)。上手即反演：
\[={\sum\limits_{i=1}^{n}}(i,n)^{x-y}i^yn^y \]
\[=n^y{\sum\limits_{i=1}^{n}}i^y(i,n)^{x-y} \]
\[=n^y{\sum\limits_{i=1}^{n}}i^y{\sum\limits_{d|(i,n)}}g(d) \]
其中，
\[g(d)={\sum\limits_{p|d}}{\mu(p)}({\frac{d}{p}})^{x-y} \]
\[原式=n^y{\sum\limits_{d|n}}g(d){\sum\limits_{i=1}^{n}}i^y[d|i] \]
\[=n^y{\sum\limits_{d|n}}g(d)d^y{\sum\limits_{i=1}^{\frac{n}{d}}}i^y \]
利用伯努利數，有
\[{\sum\limits_{i=1}^{\frac{n}{d}}}i^y={\frac{1}{y+1}}{\sum\limits_{i=0}^{y}}{\tbinom{y+1}{i}}B_i(\frac{n}{d})^{y+1-i} \]
所以原式=
\[n^y*{\frac{1}{y+1}}*{\sum\limits_{i=0}^{y}}{\tbinom{y+1}{i}}B_i{\sum\limits_{d|n}}g(d)d^y(\frac{n}{d})^{y+1-i} \]
\[=n^y*{\frac{1}{y+1}}*{\sum\limits_{i=0}^{y}}{\tbinom{y+1}{i}}B_in^{y+1-i}{\sum\limits_{d|n}}g(d)d^{i-1} \]
前面兩項均為常數；由於伯努利數滿足
\[B_i=\begin{cases} 1(i=0) \\ 1-{\sum\limits_{j=0}^{i-1}}B_j{\tbinom{i+1}{j}}{\frac{1}{i+1}} \end{cases} \]
的遞推式，所以所有的\(B_i\)可以\(O(y^2)\)預處理出來，也變成了常數。所以問題的瓶頸就落在了計算\({\sum_{d|n}}g(d)d^{i-1}\)上。由於前面已經循環過i=0到y,並且還有T組數據，所以這里必須在最多\(O(\log n)\)內解決。不過，我們能夠看出\(g(d)d^{i-1}\)是一個積性函數。所以，設它是\(h(d)\)，並設n的素因數分解式是\({\prod_{i=1}^{k}}p_k^{\alpha_{k}}\)，就可以完成下列轉化：
\[{\sum\limits_{d|n}}h(d) \]
\[={\prod\limits_{l=1}^{k}}{\sum\limits_{j=0}^{\alpha_l}}h(p_l^{j}) \]
現在我們再拆開h。
\[h(d)=g(d)d^{i-1} \]
\[={\sum\limits_{e|d}}{\mu(e)}(\frac{d}{e})^{x-y}d^{i-1} \]
對於質數p和自然數j，發現
\[h(p^j)=\begin{cases} 1(j=0) \\ p^{j(x-y+i-1)}-p^{j(x-y+i-1)+y-x} \end{cases}(j{\geq}1) \]
而且，\(x-y+i-1\)和\(y-x\)都處於\([-3001,6000]\)的范圍之中，這意味着我們可以對於所有的n的質因子p，預處理出p的-3001到6000次方，這樣的預處理總時間是\(O(T\log n*y)\)，可以接受。在j>1的時候，只需要乘上\(p^{x-y+i-1}\)，就可以從\(h(p^{j-1})\)轉移過來了。
於是，計算\({\sum_{d|n}}h(d)\)的時間復雜度降為了\(O({\sum}{\alpha_i})\)即\(O(\log n)\)。
過程中，需要對n做素因數分拆；由於n過大，所以需要使用Pollard-rho配合Miller-rabin實現，其時間復雜度約為\(O(n^{\frac{1}{4}})\)。
總時間復雜度\(O(Tn^{\frac{1}{4}}+Ty\log n+y^2)\)。
代碼
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ld long double
#define rg register
#define M 3000
#define Mod 1000000007
namespace ModCalc{
inline void Inc(ll &x,ll y,ll mod){
x += y;if(x >= mod)x -= mod;
}
inline void Dec(ll &x,ll y,ll mod){
x -= y;if(x < 0)x += mod;
}
inline void Adjust(ll &x,ll mod){
x = (x % mod + mod) % mod;
}
inline void Tms(ll &x,ll y,ll mod){ //O(1)快速乘
x = x * y - (ll)((ld)x * y / mod + 1e-3) * mod;
Adjust(x,mod);
}
inline ll Add(ll x,ll y,ll mod){
Inc(x,y,mod);return x;
}
inline ll Sub(ll x,ll y,ll mod){
Dec(x,y,mod);return x;
}
inline ll Check(ll x,ll mod){
Adjust(x,mod);return x;
}
inline ll Mul(ll x,ll y,ll mod){
Tms(x,y,mod);return x;
}
}
using namespace ModCalc;
inline ll read(){
ll s = 0,ww = 1;
char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-')ww = -1;ch = getchar();}
while('0' <= ch && ch <= '9'){s = 10 * s + ch - '0';ch = getchar();}
return s * ww;
}
inline ll power(ll a,ll n,ll mod){
ll x = a,s = 1;
while(n){
if(n & 1)Tms(s,x,mod);
Tms(x,x,mod);
n >>= 1;
}
return s;
}
ll B[M+5],b[M+5];
ll jc[M+5],iv[M+5]; //階乘和階乘的逆元
inline void prepro(){ //預處理jc,iv
jc[0] = 1;
for(rg ll i = 1;i <= M + 1;i++)jc[i] = jc[i-1] * i % Mod;
iv[M+1] = power(jc[M+1],Mod - 2,Mod);
for(rg ll i = M;i >= 0;i--)iv[i] = iv[i+1] * (i + 1) % Mod;
}
inline void calcB(){
B[0] = 1;
for(rg ll i = 1;i <= M;i++){
ll cur = 0;
for(rg ll j = 0;j < i;j++)Inc(cur,jc[i+1] * iv[j] % Mod * iv[i+1-j] % Mod * B[j] % Mod,Mod);
cur = cur * jc[i] % Mod * iv[i+1] % Mod;
B[i] = Sub(1,cur,Mod);
}
}
inline ll test(ll p,ll n){
if(n == p)return 1;
ll m = n - 1;
while(!(m & 1))m >>= 1;
ll x = power(p,m,n);
if(x == 1)return 1;
for(;m < n - 1;m <<= 1){
ll y = Mul(x,x,n);
if(y == 1)return x == n - 1;
x = y;
}
return 0;
}
ll pri[9] = {2,3,5,7,11,13,17,19,23};
inline ll MR(ll n){ //Miller-rabin
if(n < 2)return 0;
if(!(n & 1))return n == 2;
for(rg ll i = 0;i < 9;i++)if(!test(pri[i],n))return 0;
return 1;
}
inline ll f(ll x,ll c,ll mod){
return Add(Mul(x,x,mod),c,mod);
}
inline ll gcd(ll a,ll b){
return b ? gcd(b,a % b) : a;
}
inline ll PR(ll n){ //Pollard-rho
ll c = 1ll * rand() % (n - 1) + 1,x = 0,y = 0,prod = 1;
for(rg ll step = 1;;step <<= 1){
for(rg ll i = 1;i <= step;i++){
y = f(y,c,n);
Tms(prod,abs(x-y),n);
if(i % 127 == 0){
ll d = gcd(prod,n);
if(d > 1)return d;
}
}
ll d = gcd(prod,n);
if(d > 1)return d;
x = y;
}
}
vector<ll>fac; //PR算法分拆出的無序質因子
vector<pair<ll,ll> >Pfac; //整理好后的質因子和對應的次數
inline void divide(ll n){ //分拆n
if(n == 1)return;
if(MR(n)){
fac.push_back(n);
return;
}
ll d = PR(n);
while(d == n)d = PR(n);
divide(d);
divide(n / d);
}
ll pw[60+5][3*M+5]; //n的所有素因子的-3001~6000次方
int main(){
srand((int)time(NULL));
prepro();
calcB();
ll T = read();
while(T--){
ll n = read(),x = read(),y = read();
fac.resize(0);
Pfac.resize(0);
divide(n);
sort(fac.begin(),fac.end());
ll p = 0,v = 0;
for(rg ll i = 0;i < fac.size();i++){
if(fac[i] != p){
if(p)Pfac.push_back(make_pair(p,v));
p = fac[i],v = 1;
}
else v++;
}
if(p)Pfac.push_back(make_pair(p,v));
for(rg ll i = 0;i < Pfac.size();i++){ //預處理pw
ll p = Pfac[i].first,vp = power(p % Mod,Mod - 2,Mod);
pw[i][3001] = 1;
for(rg ll j = 1;j <= 6000;j++)pw[i][j+3001] = pw[i][j+3000] * (p % Mod) % Mod;
for(rg ll j = -1;j >= -3001;j--)pw[i][j+3001] = pw[i][j+3002] * vp % Mod;
}
ll ans = 0;
for(rg ll i = 0;i <= y;i++){
ll b = B[i] * jc[y] % Mod * iv[i] % Mod * iv[y+1-i] % Mod;
b = b * power(n % Mod,y + 1 - i,Mod) % Mod;
for(rg ll l = 0;l < Pfac.size();l++){
ll s = 1,w = pw[l][x-y+i-1+3001];
ll cur = Sub(1,pw[l][y-x+3001],Mod);
for(rg ll j = 1;j <= Pfac[l].second;j++){
cur = cur * w % Mod;
Inc(s,cur,Mod);
}
b = b * s % Mod;
}
Inc(ans,b,Mod);
}
ans = ans * power(n % Mod,y,Mod) % Mod;
cout << ans << endl;
}
return 0;
}