数表（bzoj 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

/*
    在PoPoQQQ的ppt上看的，感觉各种化简很神。 
*/
#include<cstdio>
#include<iostream>
#include<algorithm>
#define N 1000010
using namespace std;
int prime[N],num,mul[N],t[N],mark[N],ans[N],mx;
struct node{int n,m,a,id;}q[N];
pair<int,int> F[N];
bool operator<(node a,node b){
    return a.a<b.a;
}
void add(int x,int val){
    for(int i=x;i<=mx;i+=i&-i) t[i]+=val;
}
int query(int x){
    int tmp=0;
    for(int i=x;i;i-=i&-i) tmp+=t[i];
    return tmp;
}
void get_prime(){
    mul[1]=1;
    for(int i=2;i<=mx;i++){
        if(!mark[i]) prime[++num]=i,mul[i]=-1;
        for(int j=1;j<=num&&i*prime[j]<=mx;j++){
            mark[i*prime[j]]=1;
            mul[i*prime[j]]=-mul[i];
            if(i%prime[j]==0){mul[i*prime[j]]=0;break;}
        }
    }
    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;
}
void solve(int x){
    int id=q[x].id,n=q[x].n,m=q[x].m;
    for(int i=1,j;i<=q[x].n;i=j+1){
        j=min(n/(n/i),m/(m/i));
        ans[id]+=(n/i)*(m/i)*(query(j)-query(i-1));
    }
}
int main(){
    int Q;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);
        mx=max(mx,q[i].n);q[i].id=i;
    }
    get_prime();
    sort(q+1,q+Q+1);
    sort(F+1,F+mx+1);
    int now=0;
    for(int i=1;i<=Q;i++){
        while(now+1<=mx&&F[now+1].first<=q[i].a){
            now++;
            for(int j=F[now].second;j<=mx;j+=F[now].second)
                add(j,F[now].first*mul[j/F[now].second]);
        }
        solve(i);
    }
    for(int i=1;i<=Q;i++)
        printf("%d\n",ans[i]&0x7fffffff);
    return 0;
}