2020杭電多校第六場 7. A Very Easy Math Problem

題意：

求

\[\sum_{a_1=1}^n\sum_{a_2=1}^n...\sum_{a_x=1}^n(\prod_{j=1}^xa_j^k)f(gcd(a_1,a_2,...,a_x))*gcd(a_1,a_2,...,a_x) \]

其中\(f(x)=0\)，如果\(x\)有平方因子；否則\(f(x)=1\)。

\(t\)組詢問，資料範圍為\(t\leq 10^4,1\leq k\leq 10^9,1\leq x\leq 10^9,1\leq n\leq 2\times 10^5\)。

思路：

先做一些分析，如何求出\(f\)

對於一個數\(x\)，通過算術基本定理可知一定可以寫成若干個質數相乘的情況，如果某一個質數的冪次\(\geq 2\)，那麼他一定是有完全平方因子的，同時它的莫比烏斯函式等於\(0\)。

否則它的莫比烏斯函式等於\(1\)或\(-1\)。

那麼\(f=\mu^2\)。

接著來推公式吧。

\[\sum_{a_1=1}^n\sum_{a_2=1}^n...\sum_{a_x=1}^n(\prod_{j=1}^xa_j^k)f(gcd(a_1,a_2,...,a_x))*gcd(a_1,a_2,...,a_x) \]

列舉\(gcd=d\)，那麼有

\[\sum_{a_1=1}^n\sum_{a_2=1}^n...\sum_{a_x=1}^n(\prod_{j=1}^xa_j^k)f(d)*d[gcd(a_1,a_2,...,a_x)==d] \]

提取\(d\)，即列舉\(a_id\)：

\[\sum_{d=1}^n\sum_{a_1=1}^{\frac{n}{d}}\sum_{a_2=1}^{\frac{n}{d}}...\sum_{a_x=1}^{\frac{n}{d}}(\prod_{j=1}^xa_j^k)d^{kx}f(d)*d[gcd(a_1,a_2,...,a_x)==1] \]

反演：

\[\sum_{d=1}^n\sum_{a_1=1}^{\frac{n}{d}}\sum_{a_2=1}^{\frac{n}{d}}...\sum_{a_x=1}^{\frac{n}{d}}(a_1a_2...a_x)^kd^{kx}\sum_{j|gcd}f(d)*d \]

提取\(j\)，跟上面提取\(d\)一樣的原理：

\[\sum_{d=1}^n\sum_{j=1}^{\frac{n}{d}}\sum_{a_1=1}^{\frac{n}{dj}}\sum_{a_2=1}^{\frac{n}{dj}}...\sum_{a_x=1}^{\frac{n}{dj}}(a_1a_2...a_x)^k(dj)^{kx}f(d)*d \]

整理一下：

\[\sum_{d=1}^n\sum_{j=1}^{\frac{n}{d}}\mu(j)(dj)^{kx}f(d)d(1^k+2^k+...+(\frac{n}{dj})^k)^x \]

這一步將乘積展開了，可能有點迷惑，所以做一些解釋。

如果只有一個求和\(\sum_{i=1}^ni\)，那麼就是\((1+2+3+...+n)\)。

如果再加一個求和變成\(\sum_{i1=1}^n\sum_{i2=1}^ni1\times i2=\sum_{i1=1}^ni1\sum_{i2=1}^ni2\)。

那麼是不是就可以寫成\((1+2+3+...+n)^2\)了。

這時候如果只有一組詢問就可以寫了，但是多組，繼續化簡。

令\(dj=T\)

\[\sum_{T=1}^n\sum_{j|T}\mu(j)(T)^{kx}f(\frac{T}{j})\frac{T}{j}(1^k+...+(\frac{n}{T}^k))^x \]

令

\[g(T)=\sum_{j|T}\mu(j)(T)^{kx}f(\frac{T}{j})\frac{T}{j} \]

很明顯可以用倍數法預處理。

這樣在每次詢問的時候就可以對\(n\)分塊處理，時間複雜度\(O(T\sqrt{n}+nlogn)\)。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

template<typename T>inline void read(T &x){
	x=0; static int p;p=1;
	static char c;c=getchar();
	while(!isdigit(c)){if(c=='-')p=-1;c=getchar();}
	while(isdigit(c)) {x=(x<<1)+(x<<3)+(c-48);c=getchar();}
	x*=p;
}

const int maxn = 2e5+10;
const int mod = 1e9+7;
bool vis[maxn];
int primes[maxn], cnt;
int mu[maxn];
int f[maxn];
ll k, x;
ll sum[maxn];

ll dj[maxn];
ll ndj[maxn];

inline ll qmi(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b&1) res = (res*a)%mod;
        a = (a*a)%mod;
        b >>= 1;
    }
    return res%mod;
}

inline void init(int n)
{
    mu[1] = 1;
    for(register int i = 2; i <= n; i++)
    {
        if(!vis[i])
        {
            primes[++cnt] = i;
            mu[i] = -1;
        }
        for(register int j = 1; primes[j] <= n/i; j++)
        {
            vis[primes[j]*i] = 1;
            if(i%primes[j] == 0) break;
            else mu[i*primes[j]] = -mu[i];
        }
    }
    for(register int i = 1; i <= n; i++)
        f[i] = mu[i]*mu[i];
}

int n;
inline void solve()
{
    read(n);

    ll ans = 0;
    for(int l = 1, r; l <= n; l = r+1)
    {
        r = n/(n/l);
        ans += ((sum[r]-sum[l-1])%mod+mod)%mod*ndj[n/l]%mod;
        ans = (ans%mod+mod)%mod;
    }

    printf("%lld\n", ans);
}

int main()
{
    init(200000+3);
    int _;read(_);
    read(k); read(x);
    for(int i = 1; i <= 2e5; i++)
    {
        dj[i] = qmi(i, (k*x)%(mod-1));
        ndj[i] = qmi(i, k);
    }
    for(int i = 1; i <= 2e5; i++)
    {
        ndj[i] = ndj[i]+ndj[i-1];
        ndj[i] = (ndj[i]%mod+mod)%mod;
    }
    for(int i = 1; i <= 2e5; i++)
        ndj[i] = qmi(ndj[i], x)%mod;

    for(int i = 1; i <= 200000; i++)
        for(int j = 1; j <= 200000/i; j++)
        {
            int T = i*j;
            sum[T] += mu[j]*dj[T]%mod*f[i]%mod*i%mod;
            sum[T] = (sum[T]%mod+mod)%mod;
        }
        
    for(int i = 1; i <= 200000; i++)
    {
        sum[i] += sum[i-1];
        sum[i] = (sum[i]%mod+mod)%mod;
    }

    while(_--)solve();
    return 0;
}