P3700 [CQOI2017] 小 Q 的表格 题解

· · 题解

UPD20231121

现在看当年写的这个题解，感觉真的太稚嫩了，趁没人看到这个黑历史重置一遍。

思路

蒟蒻没有什么数论基础，决定乱搞，找出了题目的一些性质。

我们在修改某一个位置的数的时候，表格中一定会有其他一些位置会同步修改。

先尝试找出会同步修改的位置。

一眼看不出来，于是就在 Excel 里面弄了张图，大概长这样：

修改任意绿色位置，所有绿色位置都会被修改。

然后我们发现，对于位置 (a,b)，如果 \gcd (a,b)=1，那么这个位置就应该修改。

分析条件 2：

b \times f(a,a+b) = (a+b) \times f(a,b)

即修改位置 (a,b) 的时候，位置 (a,a+b) 也会被修改。

然后我们来看到条件 1：

f(a,b) = f(b,a)

结合条件 2 可以转换为：

b \times f(a+b,a) = (a+b) \times f(a,b)

那么假如我们修改了位置 (i,j)，那么对于 \forall t \in \mathbb{N} ^ {\star}，位置 (i,t \times i + j)，(i + t \times j,j) 也要修改。

这个东西看起来像什么？

\gcd(i,j) = \gcd(i,t \times i + j) = \gcd(i + t \times j,j)

所以说，假如修改位置 (a,b)，那么对于任意位置 (i,j)，如果 \gcd(a,b)=\gcd(i,j)，那么位置 (i,j) 会随位置 (a,b) 一起修改。

(所以这道题的题面也是互质表的一种画法，使用 Excel 自带的复制功能可以很方便地画出上文的那张表)

好了，我们已经知道修改任一位置，应该同时修改的位置有哪些了。但是到底要怎么求出数和呢？

条件 1,2 可以转换为：

简单概括一下，就是说对于同一行的两个位置，它们的数之比等于列数之比；对于同一列的两个位置，它们的数之比等于行数之比。初始时，f(i,j)=i \times j。

还能发现，当我们把某位置的值修改为原来的 t 倍，那么所有同时修改的位置的值都应该变为原来的 t 倍。题面还提到了保证所有位置的数修改后仍然为整数，难得一见的良心出题人！

那么我们考虑未修改时 k \in [1,n] 的答案，可以算出是：

\sum_{i=1}^{k} i^{3}

这个可以用自然数幂方求和公式 O(1) 算出：

\sum_{i=1}^{k} i^{3} = \frac{k^{2} \times (k+1)^{2}}{4}

然后对于修改操作，由于是一次性修改 \gcd 相等的所有位置，那么正常来说考虑莫比乌斯反演。

但是我不正常，毕竟当时考试做到这道题的时候我还不会莫反，我直接考虑使用欧拉函数！

这个修改操作就是给定 t,k，修改 \forall i,j \in [1,k] ,gcd(i,j)=t 的位置 (i,j)，边界都是 k，而两个量边界条件相等时可以直接使用欧拉函数前缀和求得，时间复杂度直接从 O(\sqrt{k}) 降到了 O(1)。酷炫！

但是这道题我们要求的是：

\sum_{i=1}^{k} \sum_{j=1}^{k} [\gcd(i,j)=1] ij

那么我们就是要求这个函数（需要注意 n=1 时 f(n) 不为整数，向下取整会变成 0，需要注意在使用该函数后加上缺失的值）：

g(n)=\frac{\varphi(n) \times n}{2}

那么原式可以化为：

\sum_{i=1}^{k}g(i) \times 2

由于每次修改后，所有修改的位置都会变为原来的值的倍数，那么对于每个 \gcd 我们只需要记录修改后的值相对于初始值的倍数是多少，最后暴力修改即可。

总复杂度 O(n+m^{2})，可以通过。但由于 \gcd 的位置分配性质，若没有刻意构造数据，复杂度会在 O(n+m \sqrt{m}) 级别，甚至更小，所以这个方法跑的飞快。

code

#include<bits/stdc++.h>
using namespace std;
const int mn=1048576,mod=1e9+7;
int m,n,k;
long long a,b,x;

vector<int> prime,chl;
int ppt;
int answ;
long long pgcd,pgcd2;
long long mpt,ans[4*mn],pans;
long long sumphin[4*mn],phi[4*mn],phin[4*mn];
long long ch[4*mn];
bool isphi[4*mn],isch[4*mn];
bitset<4*mn> isprime;

int getgcd(int x1,int x2){
    if(x1<x2) swap(x1,x2);
    while(x2>0){
        x1%=x2;
        swap(x1,x2);
    }
    return x1;
}

void euler(){
    phi[1]=1;sumphin[1]=1;phin[1]=1;sumphin[1]=1;
    isprime.set();
    for(int i=2;i<=n;i++){
        if(isprime[i]){
            prime.push_back(i);
            phi[i]=i-1;
        }
        for(int i1=0,pi;i1<prime.size();i1++){
            pi=i*prime[i1];
            if(pi>n) break;
            isprime[pi]=false;
            if(i%prime[i1]==0){
                phi[pi]=phi[i]*prime[i1];
                break;
            }
            else phi[pi]=phi[i]*phi[prime[i1]];
        }
    }

    for(int i=2;i<=n;i++){
        phin[i]=phi[i]*i/2;
        phin[i]%=mod;
        phin[i]*=i;
        phin[i]%=mod;
        sumphin[i]=sumphin[i-1]+(phin[i])*2;
        sumphin[i]%=mod;
    }
    return;
}

int main(){
    cin>>m>>n;

    for(int i=1;i<=n;i++){
        mpt+=2*i-1;
        mpt%=mod;
        ans[i]=ans[i-1]+mpt*i;
        ans[i]%=mod;
    }

    euler();

    for(int i=1;i<=m;i++){
        cin>>a>>b>>x>>k;

        pgcd=getgcd(a,b);
        pgcd2=pgcd*pgcd;
        ch[pgcd]=x/((a*b)/pgcd2);
        ch[pgcd]%=mod;
        if(!isch[pgcd]){
            chl.push_back(pgcd);
            isch[pgcd]=true;
        }

        pans=ans[k];
        for(int i1=0;i1<chl.size();i1++){
            answ=k/chl[i1];
            pans-=((sumphin[answ]*chl[i1])%mod)*chl[i1];
            pans+=sumphin[answ]*ch[chl[i1]];
            pans%=mod;
        }
        if(pans<0) pans+=mod;
        cout<<pans<<endl;
    }
    return 0;
}