[SDOI2014]数表
题目描述
有一张N*m的数表，其第i行第j列（1 < =i < =礼，1 < =j < =m）的数值为能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。
输入输出格式
输入格式：
输入包含多组数据。 输入的第一行一个整数Q表示测试点内的数据组数，接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。
输出格式：
对每组数据，输出一行一个整数，表示答案模2^31的值。
输入输出样例
输入样例#1：
2 4 4 3 10 10 5
输出样例#1：
20 148
说明
1 < =N．m < =10^5 ， 1 < =Q < =2*10^4
题解：
令F[i]为i的因数和,g[i]为gcd=i的数量
则ans=∑∑F[gcd(i,j)]
=∑dF[d]g[d]
g[d]=∑i∑j[gcd(i,j)==d] (i<=n,j<=m)
=∑∑[gcd(i,j)==1] (i<=n/d,j<=m/d)
=∑∑∑pμ(p) (p|i,p|j)
=∑pμ(p)∑∑1 (p<=n/d,i<=n/dp,j<=m/dp)
=∑pμ(p)[n/dp][m/dp]
ans=∑dF[d]∑pμ(p)[n/dp][m/dp]
令k=dp
g[d]=∑k[n/k][m/k]μ(k/d) (k<=n，d|k)
ans=∑k[n/k][m/k]∑dF[d]μ(k/d) (d|k)
令f[k]=∑dF[d]μ(k/d) (d|k)
处理处f的前缀和就行了，对于每个F[i],把它的f[x*i]全部加上F[i]*mu[x]
至于小于a的条件，离线按a排序，将F按从小到大顺序加入维护
分块就不说了
1 #include<cstdio> 2 #include<iostream> 3 #include<algorithm> 4 #include<cstring> 5 using namespace std; 6 struct data 7 { 8 int id,n,m,a; 9 }a[100001]; 10 int Mod=(1<<31); 11 int f[100001],mu[100001],ms[100001],prime[100001],tot; 12 int c[400001],p[100001],now=1,ans[100001]; 13 bool vis[100001]; 14 bool cmp2(data a,data b) 15 { 16 return a.a<b.a; 17 } 18 bool cmp(int a,int b) 19 { 20 return f[a]<f[b]; 21 } 22 void mobius() 23 {int i,j; 24 mu[1]=1; 25 ms[1]=1;f[1]=1; 26 for (i=2;i<=100000;i++) 27 { 28 if (vis[i]==0) 29 { 30 mu[i]=-1; 31 f[i]=i+1; 32 ms[i]=i; 33 prime[++tot]=i; 34 } 35 for (j=1;j<=tot,prime[j]*i<=100000;j++) 36 { 37 vis[i*prime[j]]=1; 38 if (i%prime[j]) 39 { 40 ms[i*prime[j]]=prime[j]; 41 f[prime[j]*i]=f[i]*f[prime[j]]; 42 mu[i*prime[j]]=-mu[i]; 43 } 44 else 45 { 46 ms[prime[j]*i]=ms[i]*prime[j]; 47 if (ms[i]==i) f[prime[j]*i]=(ms[prime[j]*i]*prime[j]-1)/(prime[j]-1); 48 else f[prime[j]*i]=f[i/ms[i]]*f[prime[j]*ms[i]]; 49 mu[i*prime[j]]=0; 50 break; 51 } 52 } 53 } 54 } 55 void add(int x,int d) 56 { 57 while (x<=100000) 58 { 59 c[x]+=d; 60 x+=(x&(-x)); 61 } 62 } 63 int query(int x) 64 { 65 int s=0; 66 while (x) 67 { 68 s+=c[x]; 69 x-=(x&(-x)); 70 } 71 return s; 72 } 73 void solve(int x) 74 {int j,pos,i,lasts,ss; 75 for (;f[p[now]]<=a[x].a;now++) 76 { 77 for (j=1;p[now]*j<=100000;j++) 78 add(p[now]*j,f[p[now]]*mu[j]); 79 } 80 pos=1; 81 lasts=0; 82 for (i=1;i<=a[x].n;i=pos+1) 83 { 84 pos=min(a[x].n/(a[x].n/i),a[x].m/(a[x].m/i)); 85 ss=query(pos); 86 ans[a[x].id]=(ans[a[x].id]+(a[x].n/i)*(a[x].m/i)*(ss-lasts)); 87 //cout<<ans[a[x].id]<<' '; 88 lasts=ss; 89 } 90 while (ans[a[x].id]<0) ans[a[x].id]+=Mod; 91 //cout<<endl; 92 } 93 int main() 94 {int T,i,j; 95 cin>>T; 96 for (i=1;i<=T;i++) 97 { 98 scanf("%d%d%d",&a[i].n,&a[i].m,&a[i].a); 99 if (a[i].n>a[i].m) swap(a[i].n,a[i].m); 100 a[i].id=i; 101 } 102 mobius(); 103 for (i=1;i<=100000;i++) p[i]=i; 104 sort(p+1,p+100001,cmp); 105 sort(a+1,a+T+1,cmp2); 106 for (i=1;i<=T;i++) 107 solve(i); 108 for (i=1;i<=T;i++) 109 printf("%d\n",ans[i]); 110 }
