2627: JZPKIL

题意：求
\[ \sum_{i=1}^n (n,i)^x [i,n]^y,\ [i,n] = lcm(i,n) \]
\(n \le 10^{18},\ x,y\le 3000\)

本题带来了一种新技巧，n太大，转化成一个积性函数然后求这个积性函数，质因子分解利用积性，这养只与质因子的数量和指数有关

官方题解清橙上有

首先套路推♂倒
\[ n^y \sum_{d \mid n} d^x \sum_{e \mid \frac{n}{d}} \mu(e) e^y \sum_{i=1}^{\frac{n}{de}} i^{y+1-k} \\ \]
这里使用\(D=de\)然后整除分块并不好做

代入伯努利数，交换求和顺序
\[ \frac{n^y}{y+1} \sum_{k=0}^y \binom{y+1}{k} B_k^+ \sum_{d \mid n} d^x \sum_{e \mid \frac{n}{d}} \mu(e) e^y \sum_{i=1}^{\frac{n}{de}} i^{y+1-k} \\ \]
后面那块就是\((id^x *( (\mu \cdot id^y) * id^z)(n) = ( (\mu \cdot id^y) * id^x * id^z)(n)\)

考虑求这个积性函数

用Pollard-rho质因子分解n，然后求

\(( (\mu \cdot id^y) * id^x * id^z)(p^c)\)

这个很好求，\(\mu\)那里只能是\(1,p\)，剩下的直接\(O(c)\)计算就行了

质因子个数和指数都是log级的，所以复杂度大约是\(O(y^2 + n^{\frac{1}{4}}logn+ylog^2n)\)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N = 3005, P = 1e9+7, mo = P, inv2 = (P+1)/2;
inline int read(){char c=getchar(); int x=0,f=1;while(c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}while(c>='0'&&c<='9') {x=x*10+c-'0';c=getchar();}return x*f;
}namespace qwq {ll mul(ll a, ll b, ll P) {ll t = (a * b - (ll) ((long double) a / P * b + 0.5) * P);return t<0 ? t+P : t;}ll Pow(ll a, ll b, ll P) {ll ans = 1; a %= P;for(; b; b>>=1, a = mul(a, a, P)) if(b&1) ans = mul(ans, a, P);return ans;}bool witness(ll a, ll n, ll u, int t) {ll x = Pow(a, u, n), y = x;while(t--) {x = mul(x, x, n);if(x == 1 && y != 1 && y != n-1) return true;y = x;}return x != 1;}int p[10] = {2, 3, 5, 7, 11, 13, 17, 19, 23};bool m_r(ll n) {if(n == 2) return true;if(n <= 1 || ~n&1) return false;ll u = n-1, t = 0;while(~u&1) u >>= 1, t++;for(int i=0; i<9; i++) {if(p[i] == n) return true;if(witness(p[i], n, u, t)) return false;}return true;}ll gcd(ll a, ll b) {return !b ? a : gcd(b, a%b);}ll _rho(ll n, ll c) {int k = 2; ll x = rand()%n, y = x, d = 1;for(int i=1; d==1; i++) {x = (mul(x, x, n) + c) %n;d = gcd(n, x>y ? x-y : y-x);if(i == k) y = x, k <<= 1;}return d;}ll d[N];void p_r(ll n) {if(n == 1) return;if(m_r(n)) {d[++d[0]] = n; return;}ll t = n;while(t == n) t = _rho(n, rand() % (n-1) + 1);p_r(t); p_r(n/t);}
} using qwq::d;ll Pow(ll a, ll b) {ll ans = 1; a %= P;for(; b; b>>=1, a=a*a%P)if(b&1) ans=ans*a%P;return ans;
}int b[N];
ll inv[N], fac[N], facInv[N];
inline ll C(int n, int m) {return fac[n] * facInv[m] %P * facInv[n-m] %P;}
void init(int n) {inv[1] = fac[0] = facInv[0] = 1;for(int i=1; i<=n+1; i++) {if(i != 1) inv[i] = (P - P/i) * inv[P%i] %P;fac[i] = fac[i-1] * i %P;facInv[i] = facInv[i-1] * inv[i] %P;}b[0] = 1; b[1] = P-inv2;for(int m=2; m<=n; m++) if(!(m&1)) {for(int k=0; k<=m-1; k++) b[m] = (b[m] - C(m+1, k) * b[k]) %P;b[m] = b[m] * Pow(C(m+1, m), P-2) %P;if(b[m] < 0) b[m] += P;}b[1] = inv2;
}struct meow{ll p; int k;} a[N]; int m;ll p1[N], p2[N];
void solve(ll n, int x, int y) {d[0] = 0; qwq::p_r(n); sort(d+1, d+d[0]+1);m = 0; a[++m] = (meow){d[1], 1};for(int i=2; i<=d[0]; i++) {if(d[i] == d[i-1]) a[m].k++;else a[++m] = (meow){d[i], 1};}//for(int i=1; i<=m; i++) printf("factor %d  %lld %d\n", i, a[i].p, a[i].k);ll ans = 0;for(int k=0; k<=y; k++) { ll t = C(y+1, k) * b[k] %mo; int z = y + 1 - k;for(int i=1; i<=m; i++) {ll p = a[i].p, t1 = 0, t2 = 0; int c = a[i].k;p1[0] = 1; ll u = p1[1] = Pow(p, x); p2[0] = 1; ll v = p2[1] = Pow(p, z);for(int j=2; j<=c; j++) p1[j] = p1[j-1] * u %mo, p2[j] = p2[j-1] * v %mo;for(int j=0; j<=c; j++) t1 = (t1 + p1[j] * p2[c-j]) %mo;for(int j=0; j<c; j++)  t2 = (t2 + p1[j] * p2[c-1-j]) %mo;t2 = (mo - Pow(p, y)) * t2 %mo; t = t * (t1 + t2) %mo;}ans = (ans + t) %mo;}ans = ans * Pow(n, y) %mo * Pow(y+1, mo-2) %mo;printf("%lld\n", ans);
}int x, y; ll n;int main() {freopen("in", "r", stdin);init(3001);int T = read();while(T--) {scanf("%lld", &n); x = read(); y = read();solve(n, x, y);}
}