【BZOJ2154】Crash的数字表格

2154: Crash的数字表格

Time Limit: 20 Sec Memory Limit: 259 MB

Submit: 1302 Solved: 517

[Submit][Status][Discuss]

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。

HINT

Source

数论

写吐了反演之后颓了颓别的水题然后继续回来反演╮(╯▽╰)╭

首先∑i=1n∑j=1mlcm(i,j)=∑i=1n∑j=1mi∗jgcd(i,j)\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)=\sum_{i=1}^n\sum_{j=1}^m\frac{i*j} {gcd(i,j)}

显然可以枚举gcd(i,j)的值乱搞

但是上面的i*j很烦人

所以我们来换一种反演姿势不用那种水题反演了= =

通常情况下我们的反演函数都是单参数的,那为什么不开个脑洞试试双参数呢?

令f(x,y)=∑i=1x∑j=1yi∗j(gcd(i,j)=1)f(x,y)=\sum_{i=1}^x\sum_{j=1}^y i*j(gcd(i,j)=1)

对应F(x,y)=∑i=1x∑j=1yi∗jF(x,y)=\sum_{i=1}^x\sum_{j=1}^y i*j

少了一个限制条件.考虑怎样使用反演公式.

很显然F(x,y)=x(x+1)2⋅y(y+1)2F(x,y)=\frac{x(x+1)}2\cdot\frac{y(y+1)}2这样可以很快得出F(x,y)的值.当然这对我们找他们两个的关系并没有任何帮助233只是到时候程序跑起来快些虽然很关键但是目前为止并不必要

还是来从原始公式找一找关系:

如果我们想要求f,那肯定要先想枚举gcd(i,j)=1之类的东西

既然f里gcd(i,j)=1,那么i*j就可以变成一大坨互异素数的乘积,那下剩下那些i*j很显然是含有平方因子的数.这些数的μ\mu当然是0,也就是说我们可以考虑一个μ(...)F(...,...)\mu(...)F(...,...)这样形式的式子.

为了让那些μ\mu是0的数显而易见一些,我们假设f是这样的东西:

f(x,y)=∑i=1min(x,y)μ(i)F(...,...)f(x,y)=\sum_{i=1}^{min(x,y)}\mu(i)F(...,...)

考虑我在反演笔记里的反演公式形式1,也就是最普遍原始的反演公式:f(n)=∑d|nμ(d)F(nd)f(n)=\sum_{d|n}\mu(d)F(\frac n d)

当然想到搞成

f(x,y)=∑i=1min(x,y)μ(i)F(⌊xi⌋,⌊yi⌋)f(x,y)=\sum_{i=1}^{min(x,y)}\mu(i)F(\lfloor \frac x i\rfloor,\lfloor \frac y i\rfloor)

至于证明方式,跟当时普通单参数的证明方式类同.

做出这个式子之后观察一下,发现:不对!等式左右两边并不相等!(就算一眼看不出来写个程序验证一下也知道了= =)

似乎还差了一点什么…让我们调用我们敏锐的观(nao)察(bu)力,是不是这样就好了↙:

f(x,y)=∑i=1min(x,y)i2⋅μ(i)F(⌊xi⌋,⌊yi⌋)f(x,y)=\sum_{i=1}^{min(x,y)}i^2\cdot\mu(i)F(\lfloor \frac x i\rfloor,\lfloor \frac y i\rfloor)

至此这两个双参函数的反演完成.我们开始着手应对原题要求的答案.

已有:ans=∑i=1n∑j=1mi∗jgcd(i,j)ans=\sum_{i=1}^n\sum_{j=1}^m\frac{i*j} {gcd(i,j)}

我们令k=gcd(i,j)k=gcd(i,j)枚举gcd之后发现可以和上面f和F互推一样的方式发展成:

ans=∑k=1min(m,n)k2⋅f(⌊nk⌋⌊mk⌋)kans=\sum_{k=1}^{min(m,n)}\frac {k^2\cdot f(\lfloor \frac n k\rfloor \lfloor \frac m k \rfloor)}{k}

=∑k=1min(m,n)k⋅f(⌊nk⌋⌊mk⌋)=\sum_{k=1}^{min(m,n)} k\cdot f(\lfloor \frac n k\rfloor \lfloor \frac m k \rfloor)

因此用枚举除法取值的方法处理⌊nk⌋⌊mk⌋\lfloor \frac n k\rfloor \lfloor \frac m k \rfloor,这一步复杂度O(n√)O(\sqrt n),而f函数里面还有一个F(⌊xi⌋,⌊yi⌋)F(\lfloor \frac x i\rfloor,\lfloor \frac y i\rfloor),所以再对外层枚举的每一步内枚举这个取值,又是一个O(n√)O(\sqrt n),最后O(n)O(n)愉快解决(虽然常数会比较大)

在上面的公式里一直开了一天多脑洞才搞出来的我终于得到了解脱QωQ

注意:这题最蛋疼的并不是公式,而是很多细节的处理

由于答案需要用到LongLong而且有取模的过程,会出现很多蛋疼的情况

如果你像我一样手残打错个括号什么的,请享受这份调代码的愉♂悦

[code]#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define MAXN 10000010
#define P 20101009
using namespace std;
bool not_prime[MAXN];
int prime[MAXN],top;
int mu[MAXN]={0,1};
long long prev[MAXN];
long long ans;
int n,m;
void check_prime()
{
    for (int i=2;i<=max(n,m);i++)
    {
        if (!not_prime[i])
            prime[++top]=i,mu[i]=-1;
        for (int j=1;j<=top&&i*prime[j]<=max(n,m);j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
}
long long Calc(long long x)
{
    return (x%P)*((x+1)%P)%P*10050505%P;
}
long long get_f(int a,int b)
{
    int last=0;
    long long ret=0;
    for (int i=1;i<=min(a,b);i=last+1)
    {
        last=min(a/(a/i),b/(b/i));
        ret+=(Calc(a/i)*Calc(b/i)%P)*(prev[last]-prev[i-1])%P;
        ret%=P;
    }
    return ret;
}
long long get_ans(int x,int y)
{
    long long ret=0;
    int last=0;
    for (int i=1;i<=min(x,y);i=last+1)
    {
        last=min(x/(x/i),y/(y/i));
        ret+=get_f(x/i,y/i)*(Calc(last)-Calc(i-1))%P;
        ret%=P;
    }
    if (ret<0) ret=(ret+P)%P;
    return ret;
}
int main()
{
    freopen("nt2011_table.in","r",stdin);
    freopen("nt2011_table.out","w",stdout);
    scanf("%d%d",&n,&m);
    check_prime();
    for (int i=1;i<=max(n,m);i++) prev[i]=(prev[i-1]+(long long)i*i*mu[i])%P;
    long long ans=get_ans(n,m);
    cout<<ans<<endl;
}