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
4 4 3
10 10 5
Sample Output
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;
}