3529: [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)描述一组数据。

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的

然后就写完了

const
    maxn=100100;
    h=1<<31;
type
    node=record
      n,m,a,id:longint;
    end;
var
    q:array[0..maxn]of node;
    flag:array[0..maxn]of boolean;
    p,u,k,f,ans,c:array[0..maxn]of longint;
    t,tot,max:longint;
 
function min(x,y:longint):longint;
begin
    if x<y then exit(x);
    exit(y);
end;
 
procedure swap(var x,y:longint);
var
    t:longint;
begin
    t:=x;x:=y;y:=t;
end;
 
procedure swap(var x,y:node);
var
    t:node;
begin
    t:=x;x:=y;y:=t;
end;
 
procedure sort(l,r:longint);
var
    i,j:longint;
    y:int64;
begin
    i:=l;
    j:=r;
    y:=f[k[(l+r)>>1]];
    repeat
      while f[k[i]]<y do
        inc(i);
      while f[k[j]]>y do
        dec(j);
      if i<=j then
      begin
        swap(k[i],k[j]);
        inc(i);
        dec(j);
      end;
    until i>j;
    if i<r then sort(i,r);
    if j>l then sort(l,j);
end;
 
procedure pre;
var
    i,j,s:longint;
begin
    f[1]:=1;
    u[1]:=1;
    for i:=2 to max do
      begin
        if flag[i]=false then
        begin
          inc(tot);
          p[tot]:=i;
          f[i]:=i+1;
          u[i]:=-1;
        end;
        for j:=1 to tot do
          begin
            if int64(i)*p[j]>max then break;
            flag[i*p[j]]:=true;
            if i mod p[j]=0 then
              begin
                s:=p[j];
                while i mod (int64(s)*p[j])=0 do
                  s:=s*p[j];
                if s=i then f[i*p[j]]:=(s*p[j]*p[j]-1)div(p[j]-1)
                else f[i*p[j]]:=f[i div s]*f[s*p[j]];
                break;
              end
            else
              begin
                f[i*p[j]]:=f[i]*(p[j]+1);
                u[i*p[j]]:=-u[i];
              end;
          end;
      end;
    for i:=1 to max do
      k[i]:=i;
    sort(1,max);
end;
 
procedure sort2(l,r:longint);
var
    i,j,y:longint;
begin
    i:=l;
    j:=r;
    y:=q[(l+r)>>1].a;
    repeat
      while q[i].a<y do
        inc(i);
      while q[j].a>y do
        dec(j);
      if i<=j then
      begin
        swap(q[i],q[j]);
        inc(i);
        dec(j);
      end;
    until i>j;
    if i<r then sort2(i,r);
    if j>l then sort2(l,j);
end;
 
procedure init;
var
    i:longint;
begin
    read(t);
    for i:=1 to t do
      begin
        read(q[i].n,q[i].m,q[i].a);
        q[i].id:=i;
        if q[i].n>q[i].m then swap(q[i].n,q[i].m);
        if max<q[i].n then max:=q[i].n;
      end;
    sort2(1,t);
end;
 
procedure add(x,y:longint);
begin
    while x<=100000 do
      begin
        c[x]:=longint(int64(c[x])+y);
        x:=x+(x and (-x));
      end;
end;
 
function sum(x:longint):longint;
begin
    sum:=0;
    while x>0 do
      begin
        sum:=longint(int64(sum)+c[x]);
        x:=x-(x and (-x));
      end;
end;
 
procedure main;
var
    i,j,kk,s,lasta,s1,s2:longint;
begin
    lasta:=1;
    for i:=1 to t do
      begin
        while (lasta<=max) and (f[k[lasta]]<=q[i].a) do
          begin
            s:=k[lasta];
            j:=1;
            while s<=max do
              begin
                if u[j]<>0 then add(s,f[k[lasta]]*u[j]);
                inc(s,k[lasta]);
                inc(j);
              end;
            inc(lasta);
          end;
        kk:=1;
        while kk<=q[i].n do
          begin
            s1:=q[i].n div kk;
            s2:=q[i].m div kk;
            s:=min(trunc(q[i].n/s1),trunc(q[i].m/s2));
            ans[q[i].id]:=longint(int64(ans[q[i].id])+int64(longint(int64(s1)*s2))*longint(int64(sum(s))-sum(kk-1)));
            kk:=s+1;
          end;
      end;
    for i:=1 to t do
      if ans[i]<0 then writeln(ans[i]+h)
      else writeln(ans[i]);
end;
 
begin
    init;
    pre;
    main;
end.