【bzoj 2627】JZPTAB - 乱搞数学题

　　我怎么就又手贱点开了一道数学题呢？？？
　　
【题意】
　　$n\leq 10^{18},x,y\leq 3000$，100组询问，求

$$\sum_{k=1}^{n} \gcd(k,n)^x lcm(k,n)^y$$

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

$$\begin{align*}answer 　　& = \sum_{k=1}^{n} \gcd(k,n)^x lcm(k,n)^y \\ 　　& = \sum_{k=1}^{n} \sum_{d} [gcd(k,n)=d] d^x (kn/d)^y \\ 　　& = n^y \sum_{k=1}^{n} \sum_{d} [gcd(k/d,n/d)=1]d^x (k/d)^y \\ 　　& = n^y \sum_{d\mid n} d^x \sum_{k=1}^{n/d}[gcd(k,n/d)=1]k^y \end{align*}$$
　　
　　后面那个玩意比较漂亮，稍微提出来搞一搞（差不多是bzoj3601吧）
　　 $$\begin{align*} f(n) 　　& = \sum_{k=1}^ n [gcd(k,n)=1]k^y \\ 　　& = \sum_{k=1}^ n \sum_{p\mid k,p \mid n} \mu(p)k^y \\ 　　& = \sum_{p \mid n} \mu(p) p^y \sum_{k=1}^{n/p} k^y 　　\end{align*}$$

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

　　

$$\begin{align*}answer 　　& = n^y \sum_{d\mid n} d^x \sum_{p\mid (n/d)}\mu(p)p^y\sum_{k=0}^{y+1}a_k\left(\frac{n}{dp}\right)^k \\ 　　& = n^y \sum_{k=0}^{y+1}a_k \sum_{d\mid n} \left(\frac{n}{d}\right)^x\sum_{p\mid d}\mu(p)p^y\left(\frac{d}{p}\right)^k \\ 　　& = n^{x+y}\sum_{k=0}^{y+1}a_k \sum_{d\mid n}\sum_{p\mid d}\mu(p)p^{y-k}d^{k-x} 　　\end{align*}$$

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

$$g(x,y) = \sum_{d\mid n}d^x \sum_{p\mid d} \mu(p)p^{y}$$

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

$$\prod_i\left( 1+(1-p_i^y)\sum_{j=1}^{c_i}p_i^{jx} \right)$$
　　而这就是 $g(x,y)$的值。如果我们已经得到了 $n$的质因子分解，那么每个$g(x,y)$可以在 $O(y\log n)$的预处理下 $O(\sum c_i)=O(\log n)$地算出来。

　　至此，不考虑质因子分解的复杂度，我们已经在总预处理$O(y^2)$，单次查询$O(y\log n)$预处理，然后$O(y\log n)$计算答案。
　　质因子分解的时候直接用pollard-rho算法就可以了。
　　总复杂度$O(y^2+T\sqrt[4]{n}+Ty\log n)$。

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

丑的一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;
}

// beginning

namespace 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;
}