1. 程式人生 > >BZOJ3529: [Sdoi2014]數表(莫比烏斯反演,離線)

BZOJ3529: [Sdoi2014]數表(莫比烏斯反演,離線)

Description

有一張 n×m 的數表,其第 i 行第 j 列(1 <= i <= n, 1 <= j <= m)的數值為 能同時整除 i 和 j 的所有自然數之和。給定 a , 計算數表中不大於 a 的數之和。

Input

輸入包含多組資料。 輸入的第一行一個整數Q表示測試點內的資料組數 接下來Q行,每行三個整數n,m,a(|a| < =10^9)描述一組資料。 1 < =N.m < =10^5  , 1 < =Q < =2×10^4

Output

對每組資料,輸出一行一個整數,表示答案模2^31的值。

Sample Input

2
4 4 3
10 10 5

Sample Output

20
148

解題思路:

這道題就是讓我們求${\sum_{i=1}^{N}}{\sum_{j=1}^{M}}{\sigma(gcd(i,j))}({\sigma(gcd(i,j))}<=a)$ a比較讓人噁心,考慮將${\sigma(gcd(i,j))}$按次序加入答案,直接統計,就是一種離線的做法。 那麼就不需要考慮a了,答案就變成了(設n<=m) ${\sum_{i=1}^{n}}{\sum_{j=1}^{m}}{\sigma(gcd(i,j))}$
$={\sum_{d=1}^{n}}{\sigma(d)}{\sum_{d|i}}{n}{\sum_{d|j}^{m}}[gcd(i,j)==d]$ $={\sum_{d=1}^{n}}{\sigma(d)}{\sum_{i=1}^{\left \lfloor {\frac{n}{d}} \right \rfloor}}{\sum_{j=1}^{\left \lfloor {\frac{n}{d}} \right \rfloor}}{[gcd(i,j)==1]}$ $={\sum_{d=1}^{n}}{\sigma(d)}{\sum_{k=1}^{\left \lfloor {\frac{n}{d}} \right \rfloor}}{\mu(k)}{\left \lfloor {\frac{n}{dk}} \right \rfloor}{\left \lfloor {\frac{m}{dk}} \right \rfloor}$
設T=dk $={\sum_{T=1}^{n}}{\sum_{d|T}}{\sigma(d)}{\mu(\frac{T}{d})}{\left \lfloor {\frac{n}{T}} \right \rfloor}{\left \lfloor \frac{m}{T} \right \rfloor}$ 按順序加入${\sigma(d)}$就好了 取模卡我好長時間,沒想到是這種方法。 程式碼:
  1 #include<cstdio>
  2 #include<cstring>
  3 #include<algorithm>
  4 typedef long long lnt;
  5 const int N=300010;
  6 struct int_2{int a,b;bool friend operator < (int_2 x,int_2 y){if(x.a!=y.a)return x.a<y.a;return x.b<y.b;}}F[N];
  7 int prime[N];
  8 int miu[N];
  9 bool vis[N];
 10 int cnt;
 11 int line[N];
 12 struct qust{
 13     int n,m,a;
 14     int no;
 15     int ans;
 16 }q[N];
 17 int lowbit(int x)
 18 {
 19     return x&(-x);
 20 }
 21 void update(int pos,int v)
 22 {
 23     while(pos<N)
 24     {
 25         line[pos]+=v;
 26         pos+=lowbit(pos);
 27     }
 28     return ;
 29 }
 30 int query(int pos)
 31 {
 32     int ans=0;
 33     while(pos)
 34     {
 35         ans+=line[pos];
 36         pos-=lowbit(pos);
 37     }
 38     return ans;
 39 }
 40 void gtp(void)
 41 {
 42     miu[1]=1;
 43     for(int i=2;i<N;i++)
 44     {
 45         if(!vis[i])
 46         {
 47             prime[++cnt]=i;
 48             miu[i]=-1;
 49         }
 50         for(int j=1;j<=cnt&&prime[j]*i<N;j++)
 51         {
 52             int x=i*prime[j];
 53             vis[x]=true;
 54             if(i%prime[j]==0)
 55             {
 56                 miu[x]=0;
 57                 break;
 58             }
 59             miu[x]=-miu[i];
 60         }
 61     }
 62     for(int i=1;i<N;i++)
 63     {
 64         for(int j=i;j<N;j+=i)
 65             F[j].a+=i;
 66         F[i].b=i;
 67     }
 68     return ;
 69 }
 70 bool cmp(qust x,qust y)
 71 {
 72     return x.a<y.a;
 73 }
 74 bool cmq(qust x,qust y)
 75 {
 76     return x.no<y.no;
 77 }
 78 int main()
 79 {
 80     gtp();
 81     int T;
 82     scanf("%d",&T);
 83     for(int i=1;i<=T;i++)
 84     {
 85         scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a);
 86         q[i].no=i;
 87     }
 88     std::sort(q+1,q+T+1,cmp);
 89     std::sort(F+1,F+N);
 90     for(int i=1,j=1;i<=T;i++)
 91     {
 92         for(;j<N&&F[j].a<=q[i].a;j++)
 93         {
 94             for(int k=F[j].b;k<N;k+=F[j].b)
 95                 update(k,F[j].a*miu[k/F[j].b]);
 96         }
 97         int n=q[i].n,m=q[i].m;
 98         if(n>m)
 99             std::swap(n,m);
100         for(int u=1,v;u<=n;u=v+1)
101         {
102             v=std::min(n/(n/u),m/(m/u));
103             q[i].ans+=(query(v)-query(u-1))*(n/u)*(m/u);
104         }
105     }
106     std::sort(q+1,q+T+1,cmq);
107     for(int i=1;i<=T;i++)
108         printf("%d\n",q[i].ans&(0x7fffffff));
109     return 0;
110 }