b×f(a,a+b)=(a+b)∗f(a,b)
很像辗转相减法。。
那么每次修改点
(a,b)
的值,会修改所有满足
gcd(i,j)==gcd(a,b)
的点
(i,j)
的值。
记
d=gcd(a,b)
,那么
fi,j=x∗i/d∗j/d
(gcd(i,j)==d
那么可以转化为,每次修改对角线上的值。记
Numi
为
fi,i
令 G(n)=∑nI=1∑nj=1[gcd(i,j)==1]i∗j
根据 BZOJ2154: Crash的数字表格 的做法
nn√ 求 G(n) 无法接受,复杂度还是高了一些。
考虑递推。当且仅当 n|d 时, ⌊nd⌋ 比 ⌊n−1d⌋ 大1
那么 G(n)=G(n−1)+∑i|nμ(i)i2(S2(⌊ni⌋)−S2(⌊n−1i⌋)
令 H(n)=∑i|nμ(i)i2(S2(⌊ni⌋)−S2(⌊n−1i⌋)=n3∑i|nμ(i)1i
令 h(n)=∑i|nμ(i)1i
发现 h(n) 是一个积性函数。
证明一下:
对于两个互质的数 x,y
记 x 的因子集为
记 y 的因子集为
那么 x∗y 的因子集 C 为
举个例子 x=10,y=21
A=1,2,5,10 s1=4
B=1,3,7,21 s2=4
那么
显然
μ(a)1a∗μ(b)1b=μ(a∗b)1a∗b
那么相当于
h(xy)
为
A
中每个
然后就可以线性筛出
h(n)
然后就可以下底函数分块来求了。但是有修改操作,如果用树状数组,修改复杂度为
O(logn)
,但是一次查询的复杂度为
O(lognn√)
,会
T
掉。
那么可以考虑牺牲一下修改的复杂度,用分块来维护前缀和,修改变成
然后这个题终!于!做!完!了!
19s
算是给卡过去了。。好像
S2(n)=∑ni=1φ(i)∗i2
但是不理解为什么啊
QAQ
,如有神犇愿意解答感激不尽!
【代码】
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#define N 4000005
#define M 600005
#define Mod 1000000007
#define INF 0x7fffffff
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const ull base=31;
ll read()
{
ll x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return x*f;
}
int m,n,block,Mx;
int Prime[N],bl[N],R[N];bool Not_Prime[N]={1,1};
ll f[N],h[N],g[N],Inv[N],sum[2005][2005],Sum[2005];
ll SSqr(int x){
return 1LL*x*x%Mod*x%Mod;
}
ll gcd(ll a,ll b){
return b==0?a:gcd(b,a%b);
}
void Pre()
{
Inv[1]=h[1]=g[1]=1;block=sqrt(m);
for(register int i=2;i<=m;i++)
{
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;j<=Prime[0]&&i*Prime[j]<=m;j++)
{
Not_Prime[i*Prime[j]]=1;
if(i%Prime[j]==0) {
h[i*Prime[j]]=h[i];
break;
}
h[i*Prime[j]]=h[i]*h[Prime[j]]%Mod;
}
g[i]=(g[i-1]+SSqr(i)*h[i]%Mod)%Mod;
}
for(register int i=1;i<=m;i++) {
f[i]=1LL*i*i;
bl[i]=(i-1)/block+1;
Sum[bl[i]]=(Sum[bl[i]]+f[i])%Mod;
sum[bl[i]][i-R[bl[i]-1]]=(sum[bl[i]][i-R[bl[i]-1]-1]+f[i])%Mod;
if(i==n||(i%block==0)) R[bl[i]]=i;
}
Mx=bl[m];
for(int i=1;i<=Mx;i++)
Sum[i]=(Sum[i]+Sum[i-1])%Mod;
}
void Update(ll d)
{
int Bl=bl[d];
for(int i=d;i<=R[Bl];i++)
sum[Bl][i-R[Bl-1]]=(sum[Bl][i-1-R[Bl-1]]+f[i])%Mod;
Sum[Bl]=(Sum[Bl-1]+sum[Bl][R[Bl]-R[Bl-1]])%Mod;
for(int i=Bl+1;i<=Mx;i++)
Sum[i]=Sum[i-1]+sum[i][R[i]-R[i-1]];
}
ll Get_Sum(int x,int y)
{
int sx=bl[x],sy=bl[y];
if(sx==sy) return (sum[sx][y-R[sx-1]]-sum[sx][x-1-R[sx-1]]+Mod)%Mod;
return ((sum[sx][R[sx]-R[sx-1]]-sum[sx][x-1-R[sx-1]]+sum[sy][y-R[sy-1]]+Sum[sy-1]-Sum[sx])%Mod+Mod)%Mod;
}
void Solve(int x)
{
int pos;
ll ans=0;
for(int i=1;i<=x;i=pos+1)
{
pos=(x/(x/i));
ans=(ans+Get_Sum(i,pos)*g[x/i]%Mod)%Mod;
}
printf("%lld\n",ans);
}
int main()
{
n=read(),m=read();
Pre();
for(int i=1;i<=n;i++)
{
static int a,b,k;ll x;
a=read(),b=read(),x=read(),k=read();x%=Mod;
int GCD=gcd(a,b);
f[GCD]=x*Inv[a/GCD]%Mod*Inv[b/GCD]%Mod;
Update(GCD);
Solve(k);
}
return 0;
}