BZOJ 3529: [Sdoi2014]数表 数学 + 莫比乌斯反演 + 取模优化 + 线性筛 + 树状数组

http://www.lydsy.com/JudgeOnline/problem.php?id=3529

Description

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

Input

    输入包含多组数据。
    输入的第一行一个整数Q表示测试点内的数据组数，接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。

Output

    对每组数据，输出一行一个整数，表示答案模2^31的值。

Sample Input

2
4 4 3
10 10 5

Sample Output

20
148

HINT

1 < =N．m < =10^5  ， 1 < =Q < =2×10^4

同时整除i和j的所有自然数，必然是gcd(i,j)的约数咯，之和，就是所有约数和咯，而F（x）=x的约数和，这个东西是积性函数，可以o(n)线性筛得到，当然也可以暴力nlogn预处理。

接下来膜拜PoPoQQQ 的ppt就好了

推公式过程类似 求 ∑∑gcd(i,j)类似，只不过现在把gcd(i,j)换成F（gcd(i,j)） 。那么,最后得到的公式为

对于前面部分显然可以分段sqrt(n)做，只要预处理了后面的部分，那么后面的部分我们正着搞不好搞（直接求每一个d的所有 约数x的f(x)*u(d/x)），可以反着搞，类似暴力求f(i)的思路

枚举每一个i，对于这个i，暴力更新所有i的倍数，最后得到的表就是我们所要的。（复杂度nlogn）

然后由于这个限制F（i）<=a,我们可以离线所有的查询，  每次询问将所有F(i)<=a的i插入 用树状数组维护前缀和即可。

   预处理mu 复杂度o(n)，预处理f(i)复杂度nlogn

求表 f(i)*u的是 nlogn

离线查询插入i到树状数组 复杂度 nlognlogn（第一个logn是枚举i的倍数，第二个是树状数组的）

solve查询答案复杂度是q*sqrt(n)*logn（q为查询数，sqrt为分块的，logn为树状数组的）

这个取模啊，可以优化:

对2的幂取模（2^k），实际等同于和(2^k - 1)进行“位与”操作。

注意，如果不是31次的话，得到的结果应该是同余的（即会出现负的）

参考代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long  ll;
const ll p =1000000007;
const int  N=100000;
bool is_prime[N+500];
int prime[N+50];
int mu[N+50];
int tot;
/****************************/
int tree[N+50];
int lowbit(int x)
{
    return x&-x;
}
void add(int x,int v)
{
    for (int i=x; i<=N; i+=lowbit(i))
        tree[i]+=v;
}
int get(int x)
{
    int sum=0;
    for (int i=x; i; i-=lowbit(i))
        sum+=tree[i];
    return sum;
}

void Moblus()       //o(n)
{
    tot = 0;
    mu[1] = 1;
    for(ll i = 2; i < N; i++)
    {
        if(!is_prime[i])
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(ll j = 0; j < tot && i*prime[j] < N; j++)
        {
            is_prime[i*prime[j]] = 1;
            if(i % prime[j])
            {
                mu[i*prime[j]] = -mu[i];
            }
            else
            {
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
}
struct node
{
    int x;
    int id;
};
node ud[N+50];
bool cmp(node a,node b )
{
    return a.x<b.x;
}
void pre()  // nlogn+n+nlogn
{
    for (int i=1; i<=N; i++)
        for (ll j=i ; j<=N; j+=i)
            ud[j].x+=i;
    for (int i=1; i<=N; i++)
        ud[i].id=i;
    sort(ud+1,ud+1+N,cmp);
}
//找[1,n],[1,m]内gcd(i,j)为质数的对数
int solve(int n,int m )
{
    if (n>m)swap(n,m);
    int ret=0;
    for (int i=1,last; i<=n; i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        ret+=( get(last)-get(i-1))*(n/i)*(m/i) ;
    }
    return ret&( 2147483647);
}
struct query
{
    int n,m,a;
    int id;
};
query qq[21324];
bool cmp2(query a,query b)
{
    return a.a<b.a;
}
int ans[21234];
void update(int id,int x)
{
    for (int j=id,cun=1; j<=N; j+=id,cun++)
        add(j,x*mu[cun]);
}
int main()
{
     Moblus();
    pre();      //预处理f(i)约数和
    int t;
    cin>>t;
    for (int i=1; i<=t; i++)
    {
        int n,m,a;
        scanf("%d%d%d",&n,&m,&a);
        qq[i].id=i, qq[i].n=n, qq[i].m=m,  qq[i].a=a;
    }
    sort(qq+1,qq+1+t,cmp2);
    int j=1;
    for (int i=1; i<=t; i++)
    {
        int n=qq[i].n;
        int m=qq[i].m;
        int a=qq[i].a;
        while(ud[j].x<=a)  update(ud[j].id,ud[j].x)  ,j++;
        ans[qq[i].id]=solve(n,m);
    }
    for (int i=1; i<=t; i++)
        printf("%d\n",ans[i]);

    return 0;
}