BZOJ3529: [Sdoi2014]数表

BZOJ3529

看着这个a的限制看着就很讨厌，就先拿掉。
弱化版的题目就是求

$$\sum_{i=1}^n\sum_{j=1}^m \sum_{d=1}^{min(n,m)}F(d)[gcd(i,j)==d]$$
$F(i)$表示i的约数和。
稍作变形：
$$\sum_{d=1}^{min(n,m)}F(d)\sum_{i=1}^{n/d}\sum_{j=1}^{m/d} [gcd(i,j)==1]$$
$$\sum_{d=1}^{min(n,m)}F(d)\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{t|gcd(i,j)}μ(t)$$
$$\sum_{d=1}^{min(n,m)}F(d)\sum_{t=1}^{min(n/d,m/d)}μ(t)\lfloor\frac{n}{dt}\rfloor\lfloor\frac{m}{dt}\rfloor$$
令 $T=dt$
$$\sum_{T=1}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d|T}F(d) μ(\frac{T}{d})$$
预处理出 $\sum_{d|T}F(d) μ(\frac{T}{d})$就可以解决了。
再考虑加入限制 $a$
即只有$F(d)<=a$时才加入答案。
就将 $F(d)$排序，将询问离线下来按 $a$排序。每次在树状数组中暴力插入$\sum_{d|T}F(d) μ(\frac{T}{d})$
然后查询分块统计答案。
对于取模运算可以让 $int$自然溢出，最后再&2147483647就可以了。

【代码】

#include <cstdio>
#include <iostream>
#include <algorithm>
#define N 100005
#define INF 0x7fffffff
using namespace std;
typedef long long ll;
typedef pair<int,int> pa;
int read()
{
    int 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 T,mx;
int Miu[N],p[N],ans[N],szsz[N];
pa f[N];
bool Not_Prime[N];

class Query{
    public:
    int id,x,y,z;
    Query(){}
    Query(int ii,int xx,int yy,int zz){id=ii,x=xx,y=yy,z=zz;}
}Q[N];

bool operator <(Query a,Query b){
    return a.z<b.z;
}

int lowbit(int x){return x&-x;}
void Sum_Up(int x,int y){
    for(int i=x;i<=mx;i+=lowbit(i)) szsz[i]+=y;
}
int Query_Ans(int x){
    int rtn=0;
    for(int i=x;i;i-=lowbit(i)) rtn+=szsz[i];
    return rtn;
}

void Pre_Work()
{
    Miu[1]=1;
    for(int i=2;i<=mx;i++)
    {
        if(!Not_Prime[i]) p[++p[0]]=i,Miu[i]=-1;
        for(int j=1;j<=p[0]&&i*p[j]<=mx;j++)
        {
            Not_Prime[i*p[j]]=1;
            if(i%p[j]==0) {Miu[i*p[j]]=0;break;}
            Miu[i*p[j]]=-Miu[i];
        }
    }
    for(int i=1;i<=mx;i++) 
        for(int j=i;j<=mx;j+=i)
            f[j].first+=i;
    for(int i=1;i<=mx;i++) f[i].second=i;
    sort(f+1,f+1+mx);
}

void Input_Init()
{
    T=read();
    for(int i=1;i<=T;i++)
    {
        static int x,y,z;
        x=read(),y=read(),z=read();
        if(x>y) swap(x,y);
        Q[i]=Query(i,x,y,z);mx=max(mx,x);
    }
    sort(Q+1,Q+1+T);
}

int Get_Ans(int x,int y)
{
    int pos,rtn=0;
    for(int i=1;i<=x;i=pos+1)
    {
        pos=min(x/(x/i),y/(y/i));
        rtn+=(x/i)*(y/i)*(Query_Ans(pos)-Query_Ans(i-1));
    }
    return rtn;
}

void Solve()
{
    int now=0;
    for(int i=1;i<=T;i++)
    {
        while(now+1<=mx&&f[now+1].first<=Q[i].z)
        {
            now++;
            for(int j=f[now].second;j<=mx;j+=f[now].second)
                Sum_Up(j,f[now].first*Miu[j/f[now].second]);
        }
        ans[Q[i].id]=Get_Ans(Q[i].x,Q[i].y);
    }
    for(int i=1;i<=T;i++) printf("%d\n",ans[i]&INF);
}

int main()
{
    Input_Init();
    Pre_Work();
    Solve();
    return 0;
}