Crash的LCM矩阵问题-优快云博客

Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Sample Input

4 5

Sample Output

122
【数据规模和约定】
100%的数据满足N, M ≤ 107。

jzp神题

似乎跟gcd问题有点像，实际上区别较大

求sigma(lcm(x,y))(1<=x<=n,1<=y<=m)

明显的，枚举最大公约数d之后就变为sigma(x*y)(1<=x<=n/d,1<=y<=m/d,gcd(x,y)=1)=f[d]

那么gcd=d的lcm和为f[d]*d

虽然d的取值在[1,n]但是n/d和 m/d的取值是sqrt(n)级别的

就可以用解决gcd问题的类似的方法来解决这个假设l到r的n/d和m/d是相等的，我们可以将其一起统计

对答案的贡献为sigma(x*y)(1<=x<=n/d,1<=y<=m/d,gcd(x,y)=1)*(r-l+1)*(l+r)/2

接下来问题是如何快速求sigma(x*y)(1<=x<=n/d,1<=y<=m/d,gcd(x,y)=1)

设这个为g[a,b]好了

那么类似的我们可以求f[a,b,k]=sigma(xy)(1<=x<=a,1<=y<=b,gcd[x,y]=wk)

明显f[a,b,k]可以直接算

f[a,b,k]=(1+..+a/k)*k*(1+..+b/k)*k=(a/k+1)*(a/k)/2*k*(b/k+1)*(b/k)/2*k

显然g[a,b]=sigma(mob(k)*f[a,b,k])

问题来了

每个f[a,b,k]都含有k^2这一项,这个变换不连续……

但是显然的f[a,b,k]只会乘mob(k)

那我们定义o(i)=mob(i)*i*i,w(a,b,k)=(a/k+1)*(a/k)/2*(b/k+1)*(b/k)/2

g[a,b]=sigma(o(k)*w(a,b,k))

发现w(a,b,k)是分sqrt段的

o可以预处理出来

就可以sqrt级别时间求g函数了

那么复杂度为O(sqrt(n)*sqrt(n))=O(n)

program bzoj2154;
const
    maxn=20101009;
var
    n,m,i,j,k:longint;
    ans:int64;
    f:array [0..10000001] of longint;
    p:array [0..10000001] of int64;

function calc (a,b:longint):int64;inline;
var
    i,k:longint;
    ans:int64;
begin
    ans:=0;
    i:=1;
    while i<=a do
        begin
            if a div (a div i)<b div (b div i) then
                k:=a div (a div i)
                                                        else
                k:=b div (b div i);
            ans:=(ans+(p[k]-p[i-1]) mod maxn*
                        (int64(a div i+1)*(a div i) div 2 mod maxn) mod maxn*
                        (int64(b div i+1)*(b div i) div 2 mod maxn) mod maxn) mod maxn;
            i:=k+1;
        end;
    exit(ans);
end;

begin
    read(n,m);
    if n>m then
        begin
            n:=n xor m;
            m:=n xor m;
            n:=n xor m;
        end;
    for i:=2 to n do
        if f[i]=0 then
            for j:=i to n div i do
                f[i*j]:=i;
    p[1]:=1;
    for i:=2 to n do
        if f[i]=0 then p[i]:=-1
                    else
        begin
            k:=i div f[i];
            if k mod f[i] = 0 then p[i]:=0
                                    else p[i]:=-p[k];
        end;
    for i:=1 to n do
        p[i]:=p[i]*i*i mod maxn;
    for i:=2 to n do
        p[i]:=(p[i]+p[i-1]) mod maxn;
    i:=1;
    while i<=n do
        begin
            if n div (n div i)<m div (m div i) then
                k:=n div (n div i)
                                                        else
                k:=m div (m div i);
            ans:=(ans+int64(k-i+1)*(k+i) div 2 mod maxn * calc(n div i,m div i)) mod maxn;
            i:=k+1;
        end;
    ans:=(ans + maxn) mod maxn;
    writeln(ans);
end.