bzoj4815 [cqoi2017]小Q的表格

版权声明:虽然是个蒟蒻但是转载还是要说一声的哟 https://blog.csdn.net/jpwang8/article/details/79961504

Description


小Q是个程序员。
作为一个年轻的程序员,小Q总是被老C欺负,老C经常把一些麻烦的任务交给小Q来处理。每当小Q不知道如何解决时,就只好向你求助。
为了完成任务,小Q需要列一个表格,表格有无穷多行,无穷多列,行和列都从1开始标号。为了完成任务,表格里面每个格子都填了一个整数,为了方便描述,小Q把第a行第b列的整数记为f(a,b)。为了完成任务,这个表格要满足一些条件:
(1)对任意的正整数a,b,都要满足f(a,b)=f(b,a);
(2)对任意的正整数a,b,都要满足b×f(a,a+b),(a+b)×f(a,b)。
为了完成任务,一开始表格里面的数很有规律,第a行第b列的数恰好等于a×b,显然一开始是满足上述两个条件的。为了完成任务,小Q需要不断的修改表格里面的数,每当修改了一个格子的数之后,为了让表格继续满足上述两个条件,小Q还需要把这次修改能够波及到的全部格子里都改为恰当的数。由于某种神奇的力量驱使,已经确保了每一轮修改之后所有格子里的数仍然都是整数。为了完成任务,小Q还需要随时获取前k行前k列这个有限区域内所有数的和是多少,答案可能比较大,只需要算出答案mod1,000,000,007之后的结果。
对于 100%的测试点,1 ≤ m ≤ 10^4, 1 ≤ a,b,k ≤ n ≤ 4×10^6, 0 ≤ x < 10^18。

来自 https://www.luogu.org/problemnew/show/P3700

Solution


首先要扒一扒题目给出条件可知的性质
我们知道gcd(a+b,a)=gcd(a,a+b)=gcd(a,b)
可以发现改变了(x,y),所有gcd(x,y)=gcd(I,j)的(I,j)位置都会变,且令d=gcd(I,j),则f(I,j)=f(d,d)*i*j/d/d
我们可以只保留一维,对于所有gcd(i,j)=d的数对(i,j)都有f(i,j)=f(d)*ij
那么答案有

ans=d=1nf(d)i=1nj=1n[gcd(i,j)=d]i×j

ans=d=1nf(d)×d2i=1ndj=1nd[gcd(i,j)=1]i×j

到这里有两种做法,一种是根据i=1nj=1n[gcd(i,j)=1]i×ji=1nφ(i)×i2 来筛破φ,另一种比较麻烦,这里讲麻烦的方法

可以看到这里有令人兴奋的布尔表达式,二话不说套μ

ans=d=1nf(d)×d2i=1ndj=1ndijx|gcd(i,jμ(x)

拉出来
ans=d=1nf(d)×d2x=1ndμ(x)x2×S(ndx)2

g(n)=i=1nμ(i)i2×S(ni)2,可以发现g(n)是能递推的,具体说来就是g(n)=g(n1)+n3i|nμ(i)i

又令h(n)=i|nμ(i)i,不难发现这是一个积性函数,可以线性筛,顺带g也就递推出来了

注意到题目需要支持修改,但是询问的时间复杂度要套上n,因此选择十分优秀的分块使得查询O(1)修改O(√n)
考虑到gcd较小的数字更容易被修改,因此如果做后缀和会快一点(反正我卡不过bzoj但是卡过了洛谷)

Code


#include <stdio.h>
#include <algorithm>
#include <math.h>
#define rep(i,st,ed) for (register int i=st;i<=ed;++i)
#define drp(i,st,ed) for (register int i=st;i>=ed;--i)

typedef long long LL;
const int MOD=1000000007;
const int N=4000010;
const int Q=2010;

LL inv[N+5],f[N+5],g[N+5],h[N+5];
LL t[Q][Q],sum[Q];
int st[Q],ed[Q],bel[N],prime[N/10],n,m;

bool not_prime[N+5];

void mod(LL &x) {
    x-=x>=MOD?MOD:0;
}

inline LL read() {
    LL x=0,v=1; char ch=getchar();
    for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):(v),ch=getchar());
    for (;ch>='0'&&ch<='9';x=x*10+ch-'0',ch=getchar());
    return x*v;
}

inline int gcd(int x,int y) {
    return (!y)?(x):gcd(y,x%y);
}

void pre_work() {
    inv[1]=f[1]=g[1]=h[1]=1;
    rep(i,2,n) {
        f[i]=1;
        inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
            h[i]=(1-inv[i]+MOD)%MOD;
        }
        for (register int j=1,x=i*prime[j];i*prime[j]<=n&&j<=prime[0];x=i*prime[++j]) {
            not_prime[x]=1;
            if (i%prime[j]==0) {
                h[x]=h[i];
                break;
            }
            h[x]=(h[i]*h[prime[j]])%MOD;
        }
        mod(g[i]=g[i-1]+((LL)i%MOD*(LL)i%MOD*(LL)i%MOD*h[i]%MOD)%MOD);
    }
    int size=(int)ceil(sqrt(n));
    rep(i,1,n) {
        bel[i]=(i-1)/size+1;
        if (!st[bel[i]]) st[bel[i]]=i;
        ed[bel[i]]=std:: max(ed[bel[i]],i);
    }
    rep(b,1,bel[n]) {
        rep(i,st[b],ed[b]) {
            t[b][i-st[b]+1]=(t[b][i-st[b]]+f[i]*i%MOD*i%MOD)%MOD;
        }
    }
    drp(b,bel[n],1) {
        sum[b]=(sum[b+1]+t[b][ed[b]-st[b]+1])%MOD;
    }
}

inline void modify(int x) {
    int bx=bel[x];
    rep(i,x,ed[bx]) t[bx][i-st[bx]+1]=(t[bx][i-st[bx]]+f[i]*(LL)i%MOD*(LL)i%MOD)%MOD;
    drp(i,bx,1) sum[i]=(sum[i+1]+t[i][ed[i]-st[i]+1])%MOD;
}

inline LL query(int x,int y) {
    int bx=bel[x],by=bel[y];
    if (bx==by) return (t[bx][y-st[bx]+1]%MOD-t[bx][x-st[bx]]%MOD+MOD)%MOD;
    LL ret=0;
    mod(ret+=t[bx][ed[bx]-st[bx]+1]-t[bx][x-1-st[bx]+1]+MOD);
    mod(ret+=t[by][y-st[by]+1]);
    mod(ret+=sum[bx+1]-sum[by]+MOD);
    return ret;
}

inline void solve(int n) {
    LL ans=0;
    for (register int i=1,j;i<=n;i=j+1) {
        j=n/(n/i);
        mod(ans+=query(i,j)*g[n/i]%MOD);
    }
    printf("%lld\n", ans);
}

int main(void) {
    m=read(); n=read();
    pre_work();
    int a,b,k,d; LL x;
    while (m--) {
        a=read(); b=read(); x=read()%MOD; k=read();
        f[d=gcd(a,b)]=x*inv[a]%MOD*inv[b]%MOD;
        modify(d); solve(k);
    }
    return 0;
}
阅读更多

没有更多推荐了,返回首页