3529: [Sdoi2014]数表 - BZOJ
Description
有一张N×m的数表，其第i行第j列（1 < =i < =n，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
这个太卡时间了，搞了我好久，不过最后跑了8s，利用了他取模的数，C++正好可以直接爆int，然后小于0就加2的31次方就行了，我就用longint强行截取了后半部分相当于爆int
首先a[i,j]上的数我们可以看成是F[gcd(i,j)]，F[i]我们都预处理出来
然后我们要求的就是 ΣF[i]*g[i] (F[i]<=a,g[i]是gcd=i的个数)
所以我们要求的就是 ΣF[i]*Σtrunc(n/d)*trunc(m/d)*μ(d/i) (F[i]<=a,i|d)
所以我们要求的就是 Σtrunc(n/d)*trunc(m/d)*ΣF[i]*μ(d/i) (F[i]<=a,i|d)
因为trunc(n/d)*trunc(m/d)只有根号n级别的个数，所以我们要处理ΣF[i]*μ(d/i)的前缀和，我用的是树状数组
然后每次询问都可以根号n*logn回答了
但是还有一个条件，就是F[i]<=a
前缀和还是动态的，怎么办
首先把F[i]排序是肯定的，对于每一个i，影响到的是他的倍数，这个很麻烦啊
于是就离线做了，把a排序，然后暴力修改前缀和，即枚举倍数
因为排序了，所以我们最多每个i都枚举倍数一次，应该是nlognlogn的
然后就写完了
1 const 2 maxn=100100; 3 h=1<<31; 4 type 5 node=record 6 n,m,a,id:longint; 7 end; 8 var 9 q:array[0..maxn]of node; 10 flag:array[0..maxn]of boolean; 11 p,u,k,f,ans,c:array[0..maxn]of longint; 12 t,tot,max:longint; 13 14 function min(x,y:longint):longint; 15 begin 16 if x<y then exit(x); 17 exit(y); 18 end; 19 20 procedure swap(var x,y:longint); 21 var 22 t:longint; 23 begin 24 t:=x;x:=y;y:=t; 25 end; 26 27 procedure swap(var x,y:node); 28 var 29 t:node; 30 begin 31 t:=x;x:=y;y:=t; 32 end; 33 34 procedure sort(l,r:longint); 35 var 36 i,j:longint; 37 y:int64; 38 begin 39 i:=l; 40 j:=r; 41 y:=f[k[(l+r)>>1]]; 42 repeat 43 while f[k[i]]<y do 44 inc(i); 45 while f[k[j]]>y do 46 dec(j); 47 if i<=j then 48 begin 49 swap(k[i],k[j]); 50 inc(i); 51 dec(j); 52 end; 53 until i>j; 54 if i<r then sort(i,r); 55 if j>l then sort(l,j); 56 end; 57 58 procedure pre; 59 var 60 i,j,s:longint; 61 begin 62 f[1]:=1; 63 u[1]:=1; 64 for i:=2 to max do 65 begin 66 if flag[i]=false then 67 begin 68 inc(tot); 69 p[tot]:=i; 70 f[i]:=i+1; 71 u[i]:=-1; 72 end; 73 for j:=1 to tot do 74 begin 75 if int64(i)*p[j]>max then break; 76 flag[i*p[j]]:=true; 77 if i mod p[j]=0 then 78 begin 79 s:=p[j]; 80 while i mod (int64(s)*p[j])=0 do 81 s:=s*p[j]; 82 if s=i then f[i*p[j]]:=(s*p[j]*p[j]-1)div(p[j]-1) 83 else f[i*p[j]]:=f[i div s]*f[s*p[j]]; 84 break; 85 end 86 else 87 begin 88 f[i*p[j]]:=f[i]*(p[j]+1); 89 u[i*p[j]]:=-u[i]; 90 end; 91 end; 92 end; 93 for i:=1 to max do 94 k[i]:=i; 95 sort(1,max); 96 end; 97 98 procedure sort2(l,r:longint); 99 var 100 i,j,y:longint; 101 begin 102 i:=l; 103 j:=r; 104 y:=q[(l+r)>>1].a; 105 repeat 106 while q[i].a<y do 107 inc(i); 108 while q[j].a>y do 109 dec(j); 110 if i<=j then 111 begin 112 swap(q[i],q[j]); 113 inc(i); 114 dec(j); 115 end; 116 until i>j; 117 if i<r then sort2(i,r); 118 if j>l then sort2(l,j); 119 end; 120 121 procedure init; 122 var 123 i:longint; 124 begin 125 read(t); 126 for i:=1 to t do 127 begin 128 read(q[i].n,q[i].m,q[i].a); 129 q[i].id:=i; 130 if q[i].n>q[i].m then swap(q[i].n,q[i].m); 131 if max<q[i].n then max:=q[i].n; 132 end; 133 sort2(1,t); 134 end; 135 136 procedure add(x,y:longint); 137 begin 138 while x<=100000 do 139 begin 140 c[x]:=longint(int64(c[x])+y); 141 x:=x+(x and (-x)); 142 end; 143 end; 144 145 function sum(x:longint):longint; 146 begin 147 sum:=0; 148 while x>0 do 149 begin 150 sum:=longint(int64(sum)+c[x]); 151 x:=x-(x and (-x)); 152 end; 153 end; 154 155 procedure main; 156 var 157 i,j,kk,s,lasta,s1,s2:longint; 158 begin 159 lasta:=1; 160 for i:=1 to t do 161 begin 162 while (lasta<=max) and (f[k[lasta]]<=q[i].a) do 163 begin 164 s:=k[lasta]; 165 j:=1; 166 while s<=max do 167 begin 168 if u[j]<>0 then add(s,f[k[lasta]]*u[j]); 169 inc(s,k[lasta]); 170 inc(j); 171 end; 172 inc(lasta); 173 end; 174 kk:=1; 175 while kk<=q[i].n do 176 begin 177 s1:=q[i].n div kk; 178 s2:=q[i].m div kk; 179 s:=min(trunc(q[i].n/s1),trunc(q[i].m/s2)); 180 ans[q[i].id]:=longint(int64(ans[q[i].id])+int64(longint(int64(s1)*s2))*longint(int64(sum(s))-sum(kk-1))); 181 kk:=s+1; 182 end; 183 end; 184 for i:=1 to t do 185 if ans[i]<0 then writeln(ans[i]+h) 186 else writeln(ans[i]); 187 end; 188 189 begin 190 init; 191 pre; 192 main; 193 end.
