P4240 毒瘤之神的考验

题目背景

Salamander的家门口是一条长长的公路。

又是一年春天将至，Salamander发现路边长出了一排毒瘤！

Salamander想带一些毒瘤回家，但是，这时毒瘤当中钻出来了一个毒瘤之神！

毒瘤之神：你想要带毒瘤走吗？想要带走毒瘤，就必须回答我的问题！如果答不出来的话，你还是乖乖回家吧！

题目描述

毒瘤之神会问TTT次，每次给定n,m，Salamander需要回答出(sum_{i=1}^nsum_{j=1}^mvarphi(ij) mod 998244353)。

Salamander这么辣鸡当然不会做啦，于是把问题丢给了你。

输入格式

第一行包含一个正整数T。

接下来T行，每行包含两个正整数，用空格隔开，表示这次询问的n,m。

输出格式

包含T行每行一个整数表示答案。

输入输出样例

输入 #1

3
1 1
2 2
3 3

输出 #1

1
5
19

说明/提示

对于40%的数据，T=1,(n,m≤10^5)；

对于50%的数据，T≤1000,(n,m≤10^5)

对于另外10%的数据，T≤10000 , (n = m ≤ 10^5)

对于100%的数据，(T≤10^4,n,m≤10^5)。

体验推式子的乐趣

[largesum_{i=1}^nsum_{j=1}^mphi(ij)\ ]

其中 (large phi(ij) = frac{phi(i)phi(j)gcd(i,j)}{phi(gcd(i,j))})

证明的话可以 ， 考虑(phi(i) *phi(j)) 比 (phi(ij)) 多了什么东西。

由

[large phi(x) = xprod_{i=1}^k(1 - frac{1}{p}) ]

可知 ， 多的部分就是 i，j都有的质因子的(displaystyle prod_{i=1}(1 - frac{1}{p}))

再都乘上一个 gcd 就得到了

[large phi(ij) = frac{phi(i)phi(j)gcd(i,j)}{phi(gcd(i,j))} ]

然后。。。

[large{egin{aligned}Ans=&largesum_{i=1}^nsum_{j=1}^mphi(ij)\ =&sum_{i=1}^nsum_{j=1}^mfrac{phi(i)phi(j)gcd(i,j)}{phi(gcd(i,j))}\ =&sum_{d=1}^nsum_{i=1}^nsum_{j=1}^m[gcd(i,j)==d]frac{phi(i)phi(j)d}{phi(d)}\ =&sum_{d=1}^nsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]frac{phi(id)phi(jd)d}{phi(d)}\ =&sum_{d=1}^nsum_{i=1}^{n/d}sum_{j=1}^{m/d}sum_{p | i , p | j }^{n/d}mu(p)frac{phi(id)phi(jd)d}{phi(d)}\ =&sum_{d=1}^nfrac{d}{phi(d)}sum_{p=1}^{n/d}mu(p)sum_{i=1}^{n/dp}sum_{j=1}^{m/dp}phi(idp)phi(jdp)\ =&sum_{Q=1}^nsum_{d|Q}frac{d}{phi(d)}mu(Q/d)sum_{i=1}^{n/Q}sum_{j=1}^{m/Q}phi(iQ)phi(jQ) end{aligned}} ]

再设

[huge F[x] = sum_{k|x}mu{(frac{x}{k})}*frac{k}{phi(k)}\ huge G[x][y] = sum_{i=1}^yphi(i*x)\ huge S[x][y][k] = sum_{Q=1}^ksum_{d|Q}frac{d}{phi(d)}mu(Q/d)sum_{i=1}^{x}sum_{j=1}^{y}phi(iQ)phi(jQ) ]

依次处理出 ， F,G,S数组 ， 由于直接开的话肯定开不下， 写俩指针动态申请就行了 ，

然后 ， S 预处理的范围可以自己定 ， 最好是 35

对于 ，n / x > 35 可以直接暴力 ， 这样的 x < n / 35 比 1e4 还小 ， 没毛病。

之后的部分可以整除分块 ， x , y 就是整除分块的值， 总之就是搞一搞就好了。

#include<iostream>
#include<cstdio>
using namespace std;
#define int long long
const int mod = 998244353 , N = 1e5+10 , B = 35;
int n , m , tot; // S[B][B] 数组越界
int F[N] , *G[N] , *S[B+1][B+1] , inv[N] , mu[N] , phi[N] , prime[N] , vis[N];
void Init(int n)
{
    inv[1] = mu[1] = phi[1] = 1;
    for(int i = 2 ; i <= n ; ++i)
    {
        if(!vis[i]) prime[++tot] = i , mu[i] = -1 , phi[i] = i - 1;
        for(int j = 1 ; j <= tot && i * prime[j] <= n ; ++j)
        {
            vis[i*prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                mu[i*prime[j]] = 0; phi[i*prime[j]] = prime[j] * phi[i];
                break;
            }
            mu[i*prime[j]] = -mu[i]; phi[i*prime[j]] = phi[i] * phi[prime[j]]; //智障错误!!!!
        }
    }

    for(int i = 2 ; i <= n ; ++i) inv[i] = mod - (mod / i) * inv[mod % i] % mod;

    for(int i = 1 ; i <= n ; ++i)
        for(int j = 1 ; i * j <= n ; ++j)
            F[i*j] = (1LL * F[i*j] + 1LL * i * inv[phi[i]] % mod * mu[j] % mod) % mod;
    for(int i = 1 ; i <= n ; ++i)
    {
        G[i] = new int [n / i + 5];
        G[i][0] = 0;
        for(int j = 1 ; j <= n / i ; ++j)
            G[i][j] = (1LL * G[i][j-1] + phi[i*j]) % mod;
    }

    for(int i = 1 ; i <= B ; ++i)
        for(int j = 1 ; j <= B ; ++j)
        {
            int len = n / max(i , j);
            S[i][j] = new int [len + 5];
            S[i][j][0] = 0;
            for(int k = 1 ; k <= len ; ++k)
               S[i][j][k] = (1LL * S[i][j][k-1] + 1LL * F[k] * G[k][i] % mod * G[k][j] % mod) % mod;
        }
    return ;
}

void solve(int n , int m)
{
    if(n > m) swap(n , m);
    long long ans = 0;
    for(int i = 1 ; i <= m / B ; ++i)
        ans = (ans + 1LL * F[i] * G[i][m / i] % mod * G[i][n / i] % mod) % mod;
    for(int i = m / B + 1 , r ; i <= n ; i = r + 1)
    {
        r = min(n , min(n / (n / i) , m / (m / i)));
        ans = (ans + (S[n / i][m / i][r] - S[n / i][m / i][i - 1]) % mod) % mod;
    }
    ans = (ans % mod + mod) % mod;
    cout << ans << "
";
    return ;
}

signed main()
{
    Init(1e5); int T; cin >> T;
    while(T --)
    {
        cin >> n >> m;
        solve(n , m);
    }
    return 0;
}