bzoj 3529: [Sdoi2014]数表 莫比乌斯反演+树状数组

题意

有一张N×m的数表，其第i行第j列（1 < =i < =礼，1 < =j < =m）的数值为能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。
1 < =N．m < =10^5 ， 1 < =Q < =2×10^4

分析

多做了几道反演题终于找到了一点所谓的套路了，也就是把式子化简之后求前缀和然后用分块处理。

直接上PoPoQQQ大佬%%%的题解

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define N 100005
#define ll long long
#define inf 0x7fffffff
using namespace std;

const int MAXN=100000;

int Q,prime[N],not_prime[N],tot,mu[N];
int c[N];
struct F{int id;ll val;}f[N];
struct query{int id,n,m,a;ll ans;}q[N];

bool cmpa(query a,query b)
{
    return a.a<b.a;
}

bool cmpval(F a,F b)
{
    return a.val<b.val;
}

bool cmpid(query a,query b)
{
    return a.id<b.id;
}

int check(int x,int y)
{
    int s=0;
    while (x%y==0)
    {
        x/=y;
        s++;
    }
    return s;
}

ll ksm(int x,int y)
{
    if (!y) return 1;
    if (y==1) return x;
    ll w=ksm(x,y/2);
    w=(ll)w*w;
    if (y%2==1) w=(ll)(w*x)&inf;
    return w;
}

void getmu(int n)
{
    mu[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i])
        {
            prime[++tot]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=1;i<=n;i++)
    {
        f[i].val=1;
        f[i].id=i;
    }
    for (int i=1;i<=tot;i++)
        for (int j=prime[i];j<=n;j+=prime[i])
        {
            int w=check(j,prime[i]);
            f[j].val=(ll)f[j].val*(ksm(prime[i],w+1)-1)/(prime[i]-1);
        }
    /*for (int i=1;i<=n;i++)
        for (int j=i;j<=n;j+=i)
            s[j]+=mu[j/i]*f[i];*/
}

void ins(int x,int y)
{
    while (x<=MAXN)
    {
        c[x]+=y;
        x+=x&(-x);
    }
}

ll find(int x)
{
    ll ans=0;
    while (x)
    {
        ans+=c[x];
        x-=x&(-x);
    }
    return ans;
}

void add(int x,int y)
{
    for (int i=x;i<=MAXN;i+=x)
        ins(i,mu[i/x]*y);
}


ll solve(int n,int m)
{
    ll ans=0,last;
    for (int i=1;i<=n;i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        ans+=(ll)(n/i)*(m/i)*(find(last)-find(i-1));
    }
    return ans;
}

int main()
{
    getmu(MAXN);
    scanf("%d",&Q);
    for (int i=1;i<=Q;i++)
    {
        scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);
        if (q[i].n>q[i].m) swap(q[i].n,q[i].m);
        q[i].id=i;
    }
    sort(q+1,q+Q+1,cmpa);
    sort(f+1,f+MAXN+1,cmpval);
    int now=1,nowq=1;
    while (nowq<=Q)
    {
        while (now<=MAXN&&f[now].val<=q[nowq].a)
        {
            add(f[now].id,f[now].val);
            now++;
        }
        while (nowq<=Q&&(now>MAXN||q[nowq].a<f[now].val))
        {
            q[nowq].ans=solve(q[nowq].n,q[nowq].m);
            nowq++;
        }
    }
    sort(q+1,q+Q+1,cmpid);
    for (int i=1;i<=Q;i++)
        printf("%lld\n",q[i].ans&inf);
    return 0;
}