数论函数学习笔记

笔记

积性函数

对于互质的两个数a，b，f(ab)=f(a)*f(b)

完全积性函数

对于任意两个数，f(ab)=f(a)f(b)

欧拉函数

欧拉函数表示[1,n]中与n互质的数，$\varphi (1) =1$
$\varphi (n) =n*(1-\frac{1}{p_1})(1-\frac{1}{p_2})...(1-\frac{1}{p_k})$
欧拉函数是积性函数，但不是完全积性函数
有一些特殊的小性质
[1,n]中与n互质的数的和是$\frac{n*\varphi (n)}{2}$
因为互质的数总是成对出现，例如k和n互质，n-k必然也和n互质，否则k不与n互质

莫比乌斯函数

如果n有平方数因子的话$\mu (n)=0$
如果n有k个质因子的话$\mu (n)=(-1)^k$
如果n=1的话$\mu(1)=1$
$\mu(n) \left\{\begin{matrix} 1& n=1\\ (-1)^k& n=p_1p_2p_3...p_k \\ 0& else \end{matrix}\right.$

除数函数

$\sigma _k(n)$表示n的所有正因子k次幂之和
$\sigma (n)$表示n的所有正因子之和

幂函数

$Id_k(n)=n^k$
$Id(n)=n$

单位函数

$\varepsilon (1)=1$其余时候全是0

狄利克雷卷积

$d(n)=\sum_{d|n}1$
$\sigma (n)=\sum_{d|n}d$
$\varphi (n)=\sum_{d|n}\mu (d)\frac{n}{d}$
这个式子可以用容斥来考虑
1的时候式子是n，然后对于$n$的一个质因子$p$，我们筛去所有$\frac{n}{p}$个，也就是p的倍数，但是会重复筛掉一些数，如果有$p_1p_2|d$那么这个d就被删掉了两次，我们就加回来，如果有$p_1p_2p_3|d$它被减去了3次，又被恢复了3次，然后再减去它……所以就得到了这个式子
$\varepsilon (n)=\sum_{d|n}\mu(d)$
设n中有k个不同的质因子
$\sum_{d|n}\mu(d) = \sum_{i=0}^{k} (-1)^i*C_{i}^{k}$
考虑杨辉三角，组合数之间分别填上+号和-号，最后的和一定是0
只有当n=1的时候，这个式子的值才是1
$n=\sum_{d|n}\varphi(d)$

莫比乌斯反演

若有两个函数f,g
$f(n)=\sum_{d|n}g(d)$
则他们也满足
$g(n)=\sum_{d|n}\mu(\frac{n}{d})f(d)$
$f(n)=\sum_{d|n}g(d)$
$g(n)=\sum_{m|n}[\frac{m}{n}=1]g(m)$
$g(n)=\sum_{m|n}\sum_{d|\frac{n}{m}}\mu(d)g(m)$
$g(n)=\sum_{d|n}\sum_{m|\frac{n}{d}}\mu(d)g(m)$
$g(n)=\sum_{d|n}\mu(d)\sum_{m|\frac{n}{d}}g(m)$
$g(n)=\sum_{d|n}\mu(d)f(\frac{n}{d})$
它还有另一种形式

若两个函数满足

$f(n)=\sum_{n|d}g(d)$
则$g(n)=\sum_{n|d}\mu(\frac{d}{n})f(d)$
$g(n)=\sum_{n|d}\mu(\frac{d}{n})\sum_{d|t}g(t)$
$g(n)=\sum_{n|t}g(t)\sum_{d|\frac{t}{n}}\mu(d)$
$g(n)=\sum_{n|t}g(t)[n==t]$
$g(n)=g(n)$

BZOJ 1257: [CQOI2007]余数之和sum

Description

给出正整数n和k，计算j(n, k)=k mod 1 + k mod 2 + k mod 3 + … + k mod n的值，其中k mod i表示k除以i的余数。例如j(5, 3)=3 mod 1 + 3 mod 2 + 3 mod 3 + 3 mod 4 + 3 mod 5=0+1+0+3+3=7

Input

输入仅一行，包含两个整数n, k。

Output

输出仅一行，即j(n, k)。

Sample Input

5 3

Sample Output

7

HINT

50%的数据满足：1<=n, k<=1000 100%的数据满足：1<=n ,k<=10^9

题解

可以拆式子

$x\%i=x- \lfloor \frac{x}{i} \rfloor*i$
而$\lfloor \frac{x}{i} \rfloor$取值只有$2\sqrt x$个
因为如果$i<\sqrt n$，那么有$\sqrt n$种不同取值，$i>=\sqrt n$时，取值只能在$\sqrt n$以内，那么值只有$2\sqrt x$

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
//#define ivorysi
using namespace std;
typedef long long ll;
int n,k;
ll ans;
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    scanf("%d%d",&n,&k);
    if(n>k) {
    ans=(ll)(n-k)*k;
    n=k;
    }
    int t,r;
    for(int i=1;i<=n;i=r+1) {
    t=k/i;r=k/t;
    if(r>n) r=n;
    ans+=(ll)(r-i+1)*k-(ll)(r-i+1)*(i+r)/2*t;
    }
    printf("%lld\n",ans);
    return 0;
}

BZOJ 2705: [SDOI2012]Longge的问题

Description

Longge的数学成绩非常好，并且他非常乐于挑战高难度的数学问题。现在问题来了：给定一个整数N，你需要求出∑gcd(i, N)(1<=i <=N)。

Input

一个整数，为N。

Output

一个整数，为所求的答案。

Sample Input

6

Sample Output

15

HINT

【数据范围】
对于60%的数据，0 对于100%的数据，0

题解

$g(n)=\sum_{d|n}d*\varphi(\frac{n}{d})$
如何快速求出这$\sqrt n$个$\varphi (d)$
我们枚举每个质因子有几个就好了，可以dfs

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 1000005
//#define ivorysi
using namespace std;
typedef long long ll;
ll n;
ll ans;
ll prime[MAXN];
int cnt[MAXN],tot;
void dfs(int dep,ll now,ll eu) {
    if(dep>tot) {
    ans+=(n/now)*eu;
    return;
    }
    dfs(dep+1,now,eu);
    ll t=1;
    siji(i,1,cnt[dep]) {
    t=t*prime[dep];
    dfs(dep+1,now*t,eu*t/prime[dep]*(prime[dep]-1));
    }
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    scanf("%lld",&n);
    int t=n;
    for(int i=2;i<=t/i;++i) {
    if(t%i==0) {
        prime[++tot]=i;
        while(t%i==0) {
        ++cnt[tot];
        t/=i;
        }
    }
    }
    if(t!=1) {
    prime[++tot]=t;
    cnt[tot]=1;
    }
    dfs(1,1,1);
    printf("%lld\n",ans);
    return 0;
}

BZOJ 2226: [Spoj 5971] LCMSum

Description

Given n, calculate the sum LCM(1,n) + LCM(2,n) + .. + LCM(n,n), where LCM(i,n) denotes the Least Common Multiple of the integers i and n.

Input

The first line contains T the number of test cases. Each of the next T lines contain an integer n.

Output

Output T lines, one for each test case, containing the required sum.

Sample Input

3
1
2
5

Sample Output

1
4
55

HINT

Constraints
1 <= T <= 300000
1 <= n <= 1000000

题解

除法比乘法慢了那么多！！！
这道题卡常！！！
$g(n)=\sum_{d|n}n\frac{\frac{n}{d}\varphi(\frac{n}{d})}{2}$
我们预处理每个数的与它互质且比它小的数的和，然后对于每个n枚举约数
拼命卡常啊QAQ

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 1000005
//#define ivorysi
using namespace std;
typedef long long ll;
ll phi[MAXN*2];
bool nonprime[MAXN*2];
int prime[MAXN*2],tot;
int T,n;
ll ans;
int read() {
    char c=getchar();
    int res=0;
    while(c<'0' || c>'9') c=getchar();
    while(c>='0' && c<='9') {
    res=res*10+c-'0';
    c=getchar();
    }
    return res;
}
void out(ll x) {
    if(x>=10) {
    out(x/10);
    }
    putchar('0'+x%10);
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    phi[1]=1;
    siji(i,2,1000000) {
    if(!nonprime[i]) {
        phi[i]=i-1;
        prime[++tot]=i;
    }
    siji(j,1,tot) {
        if((ll)prime[j]*i>1000000) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) {
        phi[i*prime[j]]=phi[i]*prime[j];
        break;
        }
        phi[i*prime[j]]=phi[i]*(prime[j]-1);
    }
    }
    siji(i,2,1000000) phi[i]=phi[i]*i/2;
    T=read();
    while(T--) {
    n=read();
    ans=0;
    for(int i=1;i*i<=n;++i) {
        if(n%i==0) {
        ans+=phi[i];
        if(i*i<n) {ans+=phi[n/i];}
        }
    }
    ans=ans*n;
    out(ans);
    putchar('\n');
    }
    return 0;
}

BZOJ 1101: [POI2007]Zap

Description

FGD正在破解一段密码，他需要回答很多类似的问题：对于给定的整数a,b和d，有多少正整数对x,y，满足x<=a
，y<=b，并且gcd(x,y)=d。作为FGD的同学，FGD希望得到你的帮助。

Input

第一行包含一个正整数n，表示一共有n组询问。（1<=n<= 50000）接下来n行，每行表示一个询问，每行三个
正整数，分别为a,b,d。（1<=d<=a,b<=50000）

Output

对于每组询问，输出到输出文件zap.out一个正整数，表示满足条件的整数对数。

Sample Input

2
4 5 2
6 4 3

Sample Output

3
2
//对于第一组询问，满足条件的整数对有(2,2)，（2,4），（4，2）。对于第二组询问，满足条件的整数对有（
6,3），（3,3）。

题解

我们转化一下
求$[\frac{a}{d},\frac{b}{d}]$之间有多少个互质的数
以下还用$a,b$表示区间左右端点
然后列出式子
$\sum_{i=1}^{a}\sum_{j=1}^{b}[gcd(i,j)==1]$
$\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{d|gcd(a,b)}\mu (d)$
$\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{d|i , d|j}\mu (d)$
$\sum_{d}\sum_{d|i}\sum_{d|j}\mu (d)$
$\sum_{d}\mu (d) \cdot \lfloor \frac{a}{d}\rfloor \cdot \lfloor\frac{b}{d}\rfloor$

然后因为$\lfloor \frac{n}{d}\rfloor$只有$2\sqrt n$种取值，且每段都对应一段连续的取值，然后这个东西的取值不超过$2\sqrt n+2\sqrt m$

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 1000005
//#define ivorysi
using namespace std;
typedef long long ll;
const int N = 50000;
int miu[N+10],prime[N+10],tot,sum[N+10],T;
bool nonprime[N+10];
ll ans=0;
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    miu[1]=1;
    siji(i,2,N) {
    if(!nonprime[i]) {
        prime[++tot]=i;
        miu[i]=-1;
    }
    siji(j,1,tot) {
        if(i*prime[j]>N) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) break;
        miu[i*prime[j]]=-miu[i];
    }
    }
    siji(i,1,N) {
    sum[i]=sum[i-1]+miu[i];
    }
    scanf("%d",&T);
    int a,b,d;
    while(T--) {
    scanf("%d%d%d",&a,&b,&d);
    a/=d;b/=d;
    ans=0;
    int n=min(a,b);
    int r;
    for(int i=1;i<=n;i=r+1) {
        r=min(a/(a/i),b/(b/i));
        if(r>n) r=n;
        ans+=(ll)(a/i)*(b/i)*(sum[r]-sum[i-1]);
    }
    printf("%lld\n",ans);
    }
    return 0;
}

BZOJ 2440: [中山市选2011]完全平方数

Description

小 X 自幼就很喜欢数。但奇怪的是，他十分讨厌完全平方数。他觉得这些
数看起来很令人难受。由此，他也讨厌所有是完全平方数的正整数倍的数。然而
这丝毫不影响他对其他数的热爱。
这天是小X的生日，小 W 想送一个数给他作为生日礼物。当然他不能送一
个小X讨厌的数。他列出了所有小X不讨厌的数，然后选取了第 K个数送给了
小X。小X很开心地收下了。
然而现在小 W 却记不起送给小X的是哪个数了。你能帮他一下吗？

Input

包含多组测试数据。文件第一行有一个整数 T，表示测试
数据的组数。
第2 至第T+1 行每行有一个整数Ki，描述一组数据，含义如题目中所描述。

Output

含T 行，分别对每组数据作出回答。第 i 行输出相应的
第Ki 个不是完全平方数的正整数倍的数。

Sample Input

4
1
13
100
1234567

Sample Output

1
19
163
2030745

HINT

对于 100%的数据有 1 ≤ Ki ≤ 10^9,T ≤ 50

题解

这道题考虑容斥原理，和$\mu (n)$的性质
二分一个x，考虑快速$[1,x]$之间有多少个非完全平方数的个数
这个就相当于，1的倍数，减去1个质数平方的倍数，加上2个质数平方的倍数……
发现这个东西符号刚好和$\mu (n)$的正负一样

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 1000005
//#define ivorysi
using namespace std;
typedef long long ll;
const int N = 100000;
int miu[N+10],prime[N+10],tot,sum[N+10],T;
bool nonprime[N+10];
ll ans=0;
ll check(ll mid) {
    ll res=0;
    for(ll i=1;i*i<=mid;++i) {
    if(miu[i]==0) continue;
    res+=(ll)miu[i]*(mid/(i*i));
    }
    return res;
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    miu[1]=1;
    siji(i,2,N) {
    if(!nonprime[i]) {
        prime[++tot]=i;
        miu[i]=-1;
    }
    siji(j,1,tot) {
        if(i*prime[j]>N) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) break;
        miu[i*prime[j]]=-miu[i];
    }
    }
    scanf("%d",&T);
    ll k;
    while(T--) {
    scanf("%lld",&k);
    ll l=1,r=10000000000LL;
    while(l<r) {
        ll mid=(l+r)>>1;
        if(check(mid)>=k) r=mid;
        else l=mid+1;
    }
    printf("%lld\n",r);
    }
    return 0;
}

BZOJ 2820: YY的GCD

Description

神犇YY虐完数论后给傻×kAc出了一题给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对kAc这种
傻×必然不会了，于是向你来请教……多组输入

Input

第一行一个整数T 表述数据组数接下来T行，每行两个正整数，表示N, M

Output

T行，每行一个整数表示第i组数据的结果

Sample Input

2
10 10
100 100

Sample Output

30
2791

HINT

T = 10000
N, M <= 10000000

题解

这道题按照之前某道1101的思路，可以很容易想到一个超时思路
枚举一个质数
$\sum_{p}^{min(n,m)}\sum_{d}^{min(n,m)}\mu(d)\cdot\lfloor \frac{n}{pd}\rfloor\cdot\lfloor\frac{m}{pd}\rfloor$
然后我们可以设$T=pd$
$\sum_{T}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\cdot\lfloor \frac{m}{T}\rfloor \sum_{p|T}\mu (\frac{T}{p})$

我们只要能求出对于$sum[T]=\sum_{p|T}\mu (\frac{T}{p})$然后处理成前缀和就很妙了，就变成1101了
然后我们枚举每个质数，去更新它的倍数的话，根据调和级数，均摊是$O(\log n)$的
然后质数一共有$O(\frac{n}{\log n})$个所以复杂度是$O(n)$可以暴力预处理

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 10000005
//#define ivorysi
using namespace std;
typedef long long ll;
ll sum[MAXN];
bool nonprime[MAXN];
int prime[MAXN],tot,n,m,miu[MAXN];
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    miu[1]=1;
    siji(i,2,10000000) {
    if(!nonprime[i]) {
        prime[++tot]=i;
        miu[i]=-1;
    }
    siji(j,1,tot) {
        if(i>10000000/prime[j]) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) break;
        miu[i*prime[j]]=-miu[i];
    }
    }
    siji(i,1,tot) {
    int t=prime[i];
    while(t<=10000000) sum[t]+=miu[t/prime[i]],t+=prime[i];
    }
    siji(i,1,10000000) sum[i]+=sum[i-1];
    int T;
    scanf("%d",&T);
    while(T--) {
    scanf("%d%d",&n,&m);
    int t=min(n,m),r;
    ll ans=0;
    for(int i=1;i<=t;i=r+1) {
        r=min(n/(n/i),m/(m/i));
        if(r>t) r=t;
        ans+=(ll)(sum[r]-sum[i-1])*(n/i)*(m/i);
    }
    printf("%lld\n",ans);
    }
    return 0;
}

BZOJ 2005: [Noi2010]能量采集

Description

栋栋有一块长方形的地，他在地上种了一种能量植物，这种植物可以采集太阳光的能量。在这些植物采集能量后，
栋栋再使用一个能量汇集机器把这些植物采集到的能量汇集到一起。 栋栋的植物种得非常整齐，一共有n列，每列
有m棵，植物的横竖间距都一样，因此对于每一棵植物，栋栋可以用一个坐标(x, y)来表示，其中x的范围是1至n，
表示是在第x列，y的范围是1至m，表示是在第x列的第y棵。 由于能量汇集机器较大，不便移动，栋栋将它放在了
一个角上，坐标正好是(0, 0)。 能量汇集机器在汇集的过程中有一定的能量损失。如果一棵植物与能量汇集机器
连接而成的线段上有k棵植物，则能量的损失为2k + 1。例如，当能量汇集机器收集坐标为(2, 4)的植物时，由于
连接线段上存在一棵植物(1, 2)，会产生3的能量损失。注意，如果一棵植物与能量汇集机器连接的线段上没有植
物，则能量损失为1。现在要计算总的能量损失。 下面给出了一个能量采集的例子，其中n = 5，m = 4，一共有20
棵植物，在每棵植物上标明了能量汇集机器收集它的能量时产生的能量损失。 在这个例子中，总共产生了36的能
量损失。

Input

仅包含一行，为两个整数n和m。

Output

仅包含一个整数，表示总共产生的能量损失。

Sample Input

【样例输入1】

5 4

【样例输入2】

3 4

Sample Output

【样例输出1】

36

【样例输出2】

20

对于100%的数据：1 ≤ n, m ≤ 100,000。

题解

这道题只要想起来某个公式，还是很简单的
这个公式就是
$n=\sum_{d|n}\varphi(d)$
这个东西可以通过数学归纳法证明，就是证明如果n是一个质数的指数幂，然后n是两个质数的指数幂……就证出来了
然后我们可以列式子
$ans=2\cdot \sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)-n\cdot m$
然后我们发现最难算的是这个（废话）
$\sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)$
$\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d|gcd(i,j)} \varphi(d)$
$\sum_{d}\sum_{d|i}\sum_{d|j} \varphi(d)$
$\sum_{d} \varphi(d)\cdot \lfloor \frac{n}{d} \rfloor \cdot \lfloor \frac{m}{d}\rfloor$
这道题就愉快的出解了

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <vector>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 100005
//#define ivorysi
using namespace std;
typedef long long ll;
ll sum[MAXN];
bool nonprime[MAXN];
int prime[MAXN],tot,n,m,phi[MAXN];
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    phi[1]=1;
    siji(i,2,100000) {
    if(!nonprime[i]) {
        prime[++tot]=i;
        phi[i]=i-1;
    }
    siji(j,1,tot) {
        if(i>100000/prime[j]) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) {
        phi[i*prime[j]]=phi[i]*prime[j];
        break;
        }
        phi[i*prime[j]]=phi[i]*(prime[j]-1);
    }
    }
    siji(i,1,100000) phi[i]+=phi[i-1];
    scanf("%d%d",&n,&m);
    int t=min(n,m),r;
    ll ans=0;
    for(int i=1;i<=t;i=r+1) {
    r=min(n/(n/i),m/(m/i));
    if(r>t) r=t;
    ans+=(ll)(phi[r]-phi[i-1])*(n/i)*(m/i);
    }
    ans=ans*2;
    ans-=(ll)n*m;
    printf("%lld\n",ans);
    return 0;
}

BZOJ 3529: [Sdoi2014]数表

Description

有一张N×m的数表，其第i行第j列（1 < =i < =礼，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

题解

$\sum F(gcd(i,j))$满足$F(gcd(i,j))<=a$然后$1<=i<=n,1<=j<=m$
我们可以这么算
设g(p)表示以p为最大公约数的对数有多少个
$g(p)=\sum_{i}^{\frac{n}{p}}\sum_{j}^{\frac{m}{p}}gcd(i,j)==1$
$g(p)=\sum_{i}^{\frac{n}{p}}\sum_{j}^{\frac{m}{p}}\sum_{d|i,d|j}\mu (d)$
$g(p)=\sum_{d}\sum_{d|i}^{\frac{n}{p}}\sum_{d|j}^{\frac{m}{p}}\mu (d)$
$g(p)=\sum_{d}\mu (d) \lfloor \frac{n}{pd}\rfloor \cdot \lfloor \frac{m}{pd} \rfloor$
$g(p)=\sum_{p|T}\mu (\frac{T}{p}) \lfloor \frac{n}{T}\rfloor \cdot \lfloor \frac{m}{T} \rfloor$

然后这个式子就可以改成
$ans=\sum_{i=1}^{min(n,m)}F(i)\cdot \sum_{i|D} \mu(\frac{D}{i})\lfloor \frac{n}{D}\rfloor \cdot \lfloor \frac{m}{D} \rfloor$
$ans=\sum_{i=D}^{min(n,m)}\lfloor \frac{n}{D}\rfloor \cdot \lfloor \frac{m}{D} \rfloor \sum_{i|D}\mu(\frac{D}{i})F(i)$
我们可以线性递推出$F(i)$
我们发现如果$i\%prime[j]!=0$ 那么约数和是$F(i)*(prime[j]+1)$
如果$i\%prime[j]==0$ 那么约数和是$(F(i)-F(\frac{i}{prime[j]}))*prime[j]+F(i)$
我们需要把$\sum_{i|D}\mu(\frac{D}{i})F(i)$维护成一个前缀和，就可以快速处理了
考虑到$F(i)<=a$我们离线处理，把询问按a排序，把$F(i)$排序，然后不断加入$F(i)$并更新到i的倍数上，用树状数组维护前缀和
然后这个东西是$O(n \log^2 n)$的
每次询问是$O(\sqrt n \log n)$
总复杂度是
$O(n \log^2 n +q\sqrt n \log n)$

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 100005
//#define ivorysi
using namespace std;
typedef long long ll;
const int mod=(1<<30)-1+(1<<30);
const int N=100000;
int Q;
int g[N+5],prime[N+5],tot,miu[N+5],pos[N+5],ans[N+5],tr[N+5];
int f[N+5];
bool nonprime[N*2];
struct node {
    int n,m,id,a;
    bool operator < (const node &rhs) const {
    return a<rhs.a;
    }
}qry[N];
bool cmp(int a,int b) {
    return f[a]<f[b];
}
int lowbit(int x) {return x&(-x);}
void Insert(int pos,int x) {
    if(x==0) return;
    while(pos<=N) {
    tr[pos]+=x;
    pos+=lowbit(pos);
    }
}
int query(int pos) {
    int res=0;
    while(pos>0) {
    res+=tr[pos];
    pos-=lowbit(pos);
    }
    return res;
}
int _query(int l,int r) {
    return query(r)-query(l-1);
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    miu[1]=1;
    f[1]=1;
    siji(i,2,N) {
        if(!nonprime[i]) {
            prime[++tot]=i;
            miu[i]=-1;
            f[i]=i+1;
        }
        siji(j,1,tot) {
            if(prime[j]>N/i) break;
            nonprime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                f[i*prime[j]]=((ll)(f[i]-f[i/prime[j]])*prime[j]+f[i]);
                break;
            }
            miu[i*prime[j]]=-miu[i];
            f[i*prime[j]]=f[i]*prime[j]+f[i];
        }
    }
    siji(i,1,N) pos[i]=i;
    sort(pos+1,pos+N+1,cmp);
    scanf("%d",&Q);
    siji(i,1,Q) {
        scanf("%d%d%d",&qry[i].n,&qry[i].m,&qry[i].a);
        qry[i].id=i;
    }
    sort(qry+1,qry+Q+1);
    int now=1;
    siji(i,1,Q) {
        while(now<=N && f[pos[now]]<=qry[i].a) {
        int t=pos[now];
        while(t<=N) {
            Insert(t,(f[pos[now]]&mod)*miu[t/pos[now]]);
            t+=pos[now];
        }
        ++now;
        }
        int t=min(qry[i].n,qry[i].m),r;
        for(int j=1;j<=t;j=r+1) {
        r=min(qry[i].n/(qry[i].n/j),qry[i].m/(qry[i].m/j));
        if(r>t) r=t;
        ans[qry[i].id]+=_query(j,r)*(qry[i].n/j)*(qry[i].m/j);
        }
        ans[qry[i].id]&=mod;
    }
    siji(i,1,Q) {
        printf("%d\n",ans[i]);
    }
    return 0;
}

BZOJ 2154: Crash的数字表格

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 ≤ 10^7。

题解

严谨认真的推一下式子
$\sum_{d}^{n}\frac{d^2\cdot F(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}{d}$
我们枚举约数，然后$F(x,y)$表示$1<=i<=x,1<=j<=y$中$i,j$互质$i\cdot j$的和
$\sum_{d}^{n}d\cdot F(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)$
然后考虑怎么求$F(x,y)$
$F(x,y)=\sum_{i=1}^{min(x,y)}\sum_{j=1}^{min(x,y)}[gcd(i,j)==1]i\cdot j$
$F(x,y)=\sum_{i=1}^{min(x,y)}\sum_{j=1}^{min(x,y)}\sum_{d|i,d|j}\mu(d) i\cdot j$
$F(x,y)=\sum_{d}\sum_{d|i}\sum_{d|j}\mu(d) i\cdot j$
$F(x,y)=\sum_{d}d^2\cdot \mu(d) Sum(\lfloor\frac{x}{d}\rfloor,\lfloor\frac{y}{d}\rfloor)$
$Sum(x,y)=\sum_{i=1}^{x}\sum_{j=1}^{y} i\cdot j=\frac{x(x+1)}{2}\cdot \frac{y(y+1)}{2}$
枚举$d$是$O(\sqrt n)$的，计算$F(x,y)$是$O(\sqrt n)$的，总体是$O(n)$的，已经可以解决本题

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 100005
#define mod 20101009
//#define ivorysi
using namespace std;
typedef long long ll;
const int N=10000000;
bool nonprime[N+5];
int prime[N+5],tot,miu[N+5],n,m,ans;
int sum(int x,int y) {
    return (1LL*x*(x+1)/2%mod*(1LL*y*(y+1)/2%mod))%mod;
}
int calc(int x,int y) {
    int t=min(x,y),r;
    int res=0;
    for(int i=1;i<=t;i=r+1) {
    r=min(x/(x/i),y/(y/i));
    if(r>t) r=t;
    res+=1LL*sum(x/i,y/i)*(miu[r]-miu[i-1])%mod;
    res=(res%mod+mod)%mod;
    }
    return res;
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    miu[1]=1;
    siji(i,2,N) {
        if(!nonprime[i]) {
        prime[++tot]=i;
        miu[i]=-1;
        }
        siji(j,1,tot) {
        if(prime[j]>N/i) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) break;
        miu[i*prime[j]]=-miu[i];
        }
    }
    siji(i,1,N) miu[i]=1LL*i*i%mod*miu[i]%mod;
    siji(i,1,N) miu[i]+=miu[i-1],miu[i]%=mod;
    scanf("%d%d",&n,&m);
    int r,t=min(n,m);
    for(int i=1;i<=t;i=r+1) {
        r=min(n/(n/i),m/(m/i));
        if(r>t) r=t;
        ans+=calc(n/i,m/i)*(1LL*(i+r)*(r-i+1)/2%mod)%mod;
        ans%=mod;
    }
    ans=(ans%mod+mod)%mod;
    printf("%d\n",ans);
    return 0;
}

BZOJ 2693: jzptab

Description

$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j) \% 100000009$($10^8+9$)

Input

一个正整数T表示数据组数
接下来T行 每行两个正整数 表示N、M

Output

T行 每行一个整数 表示第i组数据的结果

Sample Input

1
4 5

Sample Output

122

HINT

T <= 10000
N, M<=10000000

题解

我们可以优化一下上一个式子，让它可以满足多组询问
$\sum_{d=1}^{min(n,m)}d\cdot F(\lfloor \frac{n}{d}\rfloor,\lfloor \frac{m}{d}\rfloor)$
$\sum_{d=1}^{min(n,m)}d\cdot \sum_{i=1}^{min(n*m)} \mu(i)\cdot i^2 \cdot Sum(\lfloor \frac{n}{di}\rfloor,\lfloor \frac{m}{di}\rfloor)$
$\sum_{D=1}^{min(n,m)}Sum(\lfloor \frac{n}{D}\rfloor,\lfloor \frac{m}{D}\rfloor) \sum_{i|D} \mu(i)\cdot i^2 \cdot \frac{D}{i}$
我们发现只要求出
$\sum_{i|D} \mu(i)\cdot i^2 \cdot \frac{D}{i}$的前缀和我们就可以在$O(\sqrt n)$中完成一次询问
考虑到，积性函数的约数和还是积性函数
$f(ab)=f(a)\cdot f(b)$
$F(a)=\sum_{d|a}f(a)$
$F(ab)=\sum_{d|a}\sum_{h|b}f(d)\cdot f(h)$

然后我们可以在线性筛的同时求出这个函数值
当$i\%prime[j]==0$的时候我们直接让这个函数值乘上$prime[j]$即可，因为多出来的约数的贡献$\mu(i*prime[j])$全是0，所以直接考虑$\frac{D}{i}$中$D*prime[j]$

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 100005
#define mod 100000009
//#define ivorysi
using namespace std;
typedef long long ll;
const int N=10000000;
bool nonprime[N+5];
int prime[N+5],tot,func[N+5],n,m,ans,T;
int sum(int x,int y) {
    return (1LL*x*(x+1)/2%mod*(1LL*y*(y+1)/2%mod))%mod;
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    func[1]=1;
    siji(i,2,N) {
        if(!nonprime[i]) {
        prime[++tot]=i;
        func[i]=((i-1LL*i*i%mod)%mod+mod)%mod;
        }
        siji(j,1,tot) {
        if(prime[j]>N/i) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) {
            func[i*prime[j]]=1LL*prime[j]*func[i]%mod;
            break;
        }
        func[i*prime[j]]=1LL*func[prime[j]]*func[i]%mod;
        }
    }
    siji(i,1,N) func[i]+=func[i-1],func[i]%=mod;
    scanf("%d",&T);
    while(T--) {
        ans=0;
        scanf("%d%d",&n,&m);
        int t=min(n,m),r;
        for(int i=1;i<=t;i=r+1) {
        r=min(n/(n/i),m/(m/i));
        if(r>t) r=t;
        ans+=1LL*sum(n/i,m/i)*(func[r]-func[i-1])%mod;
        ans=(ans%mod+mod)%mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}

BZOJ 2694: Lcm

Description

对于任意的>1的n gcd(a, b)不是n^2的倍数
也就是说gcd(a, b)没有一个因子的次数>=2
然后求lcm(a,b)

Input

一个正整数T表示数据组数
接下来T行 每行两个正整数 表示N、M

Output

T行 每行一个整数 表示第i组数据的结果

Sample Input

4
2 4
3 3
6 5
8 3

Sample Output

24
28
233
178

HINT

HINT
T <= 10000
N, M<=4000000

题解

枚举每个合法的d
$\sum_{d}d \sum_{i=1}^{n} i^2\cdot \mu(i) \cdot Sum(\lfloor \frac{n}{di}\rfloor,\lfloor \frac{m}{di}\rfloor)$
$\sum_{T}^{min(n,m)} Sum(\lfloor \frac{n}{T}\rfloor,\lfloor \frac{m}{T}\rfloor)\sum_{d|T} (\frac{T}{d})^2\cdot d \cdot \mu(\frac{T}{d})$
然后用合法的$d$更新它的倍数，根据调和级数，那么这是$n \log n$的

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define siji(i,x,y) for(int i=(x);i<=(y);++i)
#define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
#define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
#define MAXN 100005
//#define ivorysi
using namespace std;
typedef long long ll;
const int N=4000000;
bool nonprime[N+5];
int prime[N+5],tot,miu[N+5],n,m,ans,T,func[N+5];
int sum(int x,int y) {
    return (x*(x+1)/2)*(y*(y+1)/2);
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    miu[1]=1;
    siji(i,2,N) {
        if(!nonprime[i]) {
        prime[++tot]=i;
        miu[i]=-1;
        }
        siji(j,1,tot) {
        if(prime[j]>N/i) break;
        nonprime[i*prime[j]]=1;
        if(i%prime[j]==0) break;
        miu[i*prime[j]]=-miu[i];
        }
    }
    siji(i,1,N) {
        if(miu[i]==0) continue;
        int t=i;
        while(t<=N) {
        func[t]+=t/i*t*miu[t/i];
        t+=i;
        }
    }
    siji(i,1,N) {func[i]+=func[i-1];}
    scanf("%d",&T);
    while(T--) {
        ans=0;
        scanf("%d%d",&n,&m);
        int t=min(n,m),r;
        for(int i=1;i<=t;i=r+1) {
        r=min(n/(n/i),m/(m/i));
        if(r>t) r=t;
        ans+=sum(n/i,m/i)*(func[r]-func[i-1]);
        }
        printf("%d\n",ans&((1<<30)-1));
    }
    return 0;
}