【bzoj 2627】JZPTAB - 乱搞数学题

　　我怎么就又手贱点开了一道数学题呢？？？
　　
【题意】
　　n≤1018,x,y≤3000，100组询问，求

∑k=1ngcd(k,n)xlcm(k,n)y

【解法】
　　首先一阵瞎展开瞎推。
　　

answer　　　　　　　　=∑k=1ngcd(k,n)xlcm(k,n)y=∑k=1n∑d[gcd(k,n)=d]dx(kn/d)y=ny∑k=1n∑d[gcd(k/d,n/d)=1]dx(k/d)y=ny∑d∣ndx∑k=1n/d[gcd(k,n/d)=1]ky

　　
　　后面那个玩意比较漂亮，稍微提出来搞一搞（差不多是bzoj3601吧）
　　

f(n)　　　　　　=∑k=1n[gcd(k,n)=1]ky=∑k=1n∑p∣k,p∣nμ(p)ky=∑p∣nμ(p)py∑k=1n/pky　　

　　我们意识到y次方和的y+1次多项式可以通过O(y2)预处理伯努利数然后每次O(y)得到多项式的系数，因此y次方和直接写成多项式形式，带回answer的式子，得到

　　

answer　　　　　　=ny∑d∣ndx∑p∣(n/d)μ(p)py∑k=0y+1ak(ndp)k=ny∑k=0y+1ak∑d∣n(nd)x∑p∣dμ(p)py(dp)k=nx+y∑k=0y+1ak∑d∣n∑p∣dμ(p)py−kdk−x　　

　　把后面的三项卷积给拎出来单独搞搞
　　

g(x,y)=∑d∣ndx∑p∣dμ(p)py

　　不妨从组合的角度来看，注意到这个式子中n的每个质因子的贡献都是独立的。
　　对于一个因子pi，其次数为ci，如果d中选择了次数为j的pi，其对d的贡献为pji；这个因子对μ(p)p最高只会贡献一个幂次，p对其选或不选时这个因子的总贡献为1−pi。
　　通过这样分析，我们可以得到所有因子的贡献是
　　

∏i⎛⎝1+(1−pyi)∑j=1cipjxi⎞⎠

　　而这就是g(x,y)的值。如果我们已经得到了n的质因子分解，那么每个g(x,y)可以在O(ylogn)的预处理下O(∑ci)=O(logn)地算出来。

　　至此，不考虑质因子分解的复杂度，我们已经在总预处理O(y2)，单次查询O(ylogn)预处理，然后O(ylogn)计算答案。
　　质因子分解的时候直接用pollard-rho算法就可以了。
　　总复杂度O(y2+Tn√4+Tylogn)。

　　有些细节要注意一下。。。比如求pxi的等比数列的时候直接算是比用公式算复杂度更优的，因为这样的复杂度相当于所有因子的次幂的和，总复杂度是O(logn)的，而用公式算是对每个因子都要O(logn)，就变成了O(log2n)。。。还有一些别的细节，稍微优化一下就跑得很快了。。。

丑的一b的代码：

/*    I will chase the meteor for you, a thousand times over.    Please wait for me, until I fade forever.    Just 'coz GEOTCBRL.*/#include <bits/stdc++.h>using namespace std;#define fore(i,u)  for (int i = head[u] ; i ; i = nxt[i])#define rep(i,a,b) for (int i = a , _ = b ; i <= _ ; i ++)#define per(i,a,b) for (int i = a , _ = b ; i >= _ ; i --)#define For(i,a,b) for (int i = a , _ = b ; i <  _ ; i ++)#define Dwn(i,a,b) for (int i = ((int) a) - 1 , _ = (b) ; i >= _ ; i --)#define fors(s0,s) for (int s0 = (s) , _S = s ; s0 ; s0 = (s0 - 1) & _S)#define foreach(it,s) for (__typeof(s.begin()) it = s.begin(); it != s.end(); it ++)#define mp make_pair#define pb push_back#define pii pair<int , int>#define fir first#define sec second#define MS(x,a) memset(x , (a) , sizeof (x))#define gprintf(...) fprintf(stderr , __VA_ARGS__)#define gout(x) std::cerr << #x << "=" << x << std::endl#define gout1(a,i) std::cerr << #a << '[' << i << "]=" << a[i] << std::endl#define gout2(a,i,j) std::cerr << #a << '[' << i << "][" << j << "]=" << a[i][j] << std::endl#define garr(a,l,r,tp) rep (__it , l , r) gprintf(tp"%c" , a[__it] , " \n"[__it == _])template <class T> inline void upmax(T&a , T b) { if (a < b) a = b ; }template <class T> inline void upmin(T&a , T b) { if (a > b) a = b ; }typedef long long ll;const int maxn = 3017;const int maxm = 200007;const int mod = 1000000007;const int inf = 0x7fffffff;const double eps = 1e-3;typedef int arr[maxn];typedef int adj[maxm];//int mul_tot;inline int add(int a , int b) { a += b; if (a >= mod) a -= mod; return a; }inline int mul(int a , int b) { /*++ mul_tot; */return ((ll) a * b) % mod ; }inline int dec(int a , int b) { a -= b; if (a < 0) a += mod; return a; }inline int Pow(int a , int b) {    int t = 1;    while (b) {        if (b & 1) t = mul(t , a);        a = mul(a , a) , b >>= 1;    }    return t;}#define rd RD<int>#define rdll RD<long long>template <typename Type>inline Type RD() {    Type x = 0;    int flag = 0;    char c = getchar();    while (!isdigit(c) && c != '-')        c = getchar();    (c == '-') ? (flag = 1) : (x = c - '0');    while (isdigit(c = getchar()))        x = x * 10 + c - '0';    return flag ? -x : x;}inline char rdch() {    char c = getchar();    while (!isalpha(c)) c = getchar();    return c;}// beginningnamespace IntegerDivision {    typedef long double db;    inline ll mul(ll a , ll b , ll c) {        return ( a * b - (ll) ( a / (db) c * b + eps ) * c + c ) % c;    }    inline ll Pow(ll a , ll b , ll c) {        ll t = 1;        while (b) {            if (b & 1) t = mul(t , a , c);            a = mul(a , a , c) , b >>= 1;        }        return t;    }    inline bool st(ll n , ll p) {        if (n <= p) return 1;        if (n % p == 0) return 0;        ll x = n - 1;        int s = -1;        while ( !(x & 1) )            x >>= 1 , s ++;        x = Pow(p , x , n);        if (x == 1 || x == n - 1)            return 1;        while (s --) {            x = mul(x , x , n);            if (x == 1)                return 0;            if (x == n - 1)                return 1;        }        return 0;    }    inline bool miller_rabin(ll p) {        return st(p , 2) && st(p , 3) && st(p , 5) && st(p , 7) && st(p , 11)            && st(p , 17) && st(p , 19) && st(p , 23) && st(p , 29) && st(p , 31)            && st(p , 37);    }    ll *frac;    int tot;    ll gcd(ll a , ll b) { return b ? gcd(b , a % b) : a; }    ll rho(ll n) {        ll x1 = rand() % (n - 1) + 1 , x2 = mul(x1 , x1 , n) + 1 , t = 1;        while (t == 1) {            x1 = mul(x1 , x1 , n) + 1;            x2 = mul(x2 , x2 , n) + 1;            x2 = mul(x2 , x2 , n) + 1;            t = gcd(x1 - x2 + n , n);        }        return t;    }    void divide(ll x) {        if (x == 1) return;        else if (x % 2 == 0) frac[++ tot] = 2 , divide(x / 2);        else if (miller_rabin(x)) frac[++ tot] = x;        else {            ll t = rho(x);            divide(t) , divide(x / t);        }    }    inline void work(ll x , ll *a , int &t) {        frac = a , tot = 0;        divide(x);        t = tot;    }}arr C[maxn] , B , a;inline void init() {    const int N = 3002;    rep (i , 0 , N) {        C[i][0] = 1;        rep (j , 1 , i) C[i][j] = add(C[i - 1][j] , C[i - 1][j - 1]);    }    For (i , 0 , N) {        B[i] = (i == 0);        For (j , 0 , i)            B[i] = dec(B[i] , mul(C[i + 1][j] , B[j]));        B[i] = mul(B[i] , Pow(i + 1 , mod - 2));    }}ll N;int n , x , y;void input() {    static unsigned int rand_seed = 18230742;    N = rdll() , x = rd() , y = rd();    rand_seed += (N + (x << 16) + y) + (rand_seed << 2);    srand(rand_seed);}inline void get_poly(int m) {    rep (i , 0 , m + 1) a[i] = 0;    rep (k , 0 , m) a[m + 1 - k] = mul(C[m + 1][k] , B[k]);    int x = Pow(m + 1 , mod - 2);    rep (i , 0 , m + 1)        a[i] = mul(a[i] , x);    a[m] = add(a[m] , 1);    if (!m) a[0] = dec(a[0] , 1);}int p[20] , c[20] , pw[20][3007] , iv[20][3007];int tot;inline void get_frac() {    static ll frac[97]; int cnt; tot = 0;    IntegerDivision::work(N , frac , cnt);    sort(frac + 1 , frac + cnt + 1);    rep (i , 1 , cnt) if (frac[i] != frac[i - 1])        p[++ tot] = frac[i] % mod , c[tot] = 1;    else        c[tot] ++;    rep (i , 1 , tot) {        pw[i][0] = iv[i][0] = 1;        int t = max(max(x , y) , c[i]) + 1;        rep (j , 1 , t) pw[i][j] = mul(pw[i][j - 1] , p[i]);        int v = Pow(p[i] , mod - 2);        rep (j , 1 , t) iv[i][j] = mul(iv[i][j - 1] , v);    }}inline int geo(int i , int x) {    int c = ::c[i];    int p;    if (x < 0)        p = iv[i][-x];    else        p = pw[i][ x];    if (p == 1)        return c;    int ret = 0 , t = p;    rep (j , 1 , c) {        ret = add(ret , t);        t = mul(t , p);    }    return ret;}inline int calc(int x , int y) {// assert(y >= -1);    int ret = x < 0 ? Pow(Pow(n , mod - 2) , -x) : Pow(n , x);    rep (i , 1 , tot) {        int t1 = geo(i , -x);        int t2 = dec(1 , y >= 0 ? pw[i][y] : Pow(p[i] , mod - 2));        ret = mul(ret , add(1 , mul(t1 , t2)));    }    return ret;}void solve() {    get_poly(y);    get_frac();    n = N % mod;    if (!n) return (void) (puts("0"));    int ans = 0;    int t = 1;    rep (k , 0 , y + 1) {        int tmp = calc(x - k , y - k);        ans = add(ans , mul( mul(a[k] , t) , tmp ));        t = mul(t , n);    }    ans = mul(ans , Pow(n , y));    printf("%d\n" , ans);}int main() {    #ifndef ONLINE_JUDGE        freopen("data.txt" , "r" , stdin);    #endif    init();    rep (T , 1 , rd()) {        input();        solve();    }//  gout(mul_tot);    return 0;}