数论初步
参考资料
《算法竞赛进阶指南》
基本数论概念
整除、因数、倍数什么的就不用说了吧。
数论函数
指定义域为正整数的函数。下面讨论的『函数』如非特殊说明，都指数论函数。
积性函数
定义
若有 \(\forall \gcd(a,b)=1,f(ab)=f(a)f(b)\)，则称 \(f\) 为积性函数。
若有 \(\forall a,b,f(ab)=f(a)f(b)\)，则称 \(f\) 为完全积性函数。
性质
若 \(f(x),g(x)\) 都是积性函数，则以下函数也是积性函数：
前三个证明显然，第四个其实是这两个函数的狄利克雷卷积，后面会讲。
欧几里得算法
内容及其证明
求 gcd 的常用方法。其内容为：
证明如下：
当 \(a<b\) 时，定理显然成立。
当 \(a\ge b\) 时，设 \(a=kb+r\)，其中 \(0\le r<b\)，即 \(r=a\bmod b\)。因为 \(\forall d\mid a,d\mid b\) 都有 \(d\mid kb\)，所以 \(d\mid (a-kb)\) 即 \(d\mid r\)。于是我们得到结论：\(\forall d\mid a,d\mid b,\quad d\mid (a\bmod b)\)。因此，\(b\) 与 \(a\bmod b\) 的公约数集相对于 \(a\) 与 \(b\) 的公约数集是没有改变的，那么两个集合中最大的数自然不变。
证毕。
算法应用及其时间复杂度证明
开始前，先来证明一个结论：当 \(a\ge b\) 时，\(a\bmod b<\frac a2\)。
注：证明过程中的 \(x\) 为任意不小于 1 的正实数。
\(\because \lfloor x \rfloor>\frac x2\)（一个大于等于 1 的数的一半一定小于其整数部分）
\(\therefore \lfloor\frac{a}{b}\rfloor>\frac{a}{2b}\)
\(\therefore \lfloor\frac{a}{b}\rfloor b>\frac{a}{2}\)
\(\therefore a-\lfloor\frac{a}{b}\rfloor b<\frac{a}{2}\)
\(\therefore a\bmod b<\frac{a}{2}\)
根据欧几里得算法，代码很容易得出：
int gcd(int a,int b){
return b?gcd(b,a%b):a;
}
其时间复杂度为 \(\mathcal O(\log a+\log b)\)。证明如下：
在求 \(\gcd(a,b)\) 时，会出现下面两种情况：
- \(a<b\)，此时会由
gcd(a,b)递归到gcd(b,a)。- \(a\ge b\)，此时会由
gcd(a,b)递归到gcd(b,a%b)。这样就使得其中一个数至少折半。因此这一情况最多发生 \(\mathcal O(\log a+\log b)\) 次。由于前者发生后一定会发生后者，所以前者的发生次数一定不多于后者，所以总的时间复杂度仍为 \(\mathcal O(\log a+\log b)\)。
欧拉函数
内容
欧拉函数为 \(\varphi(n)\)，表示小于等于 \(n\) 的与 \(n\) 互质的数的个数。即：
性质
性质一
若在算数基本定理中，\(n=\prod_{i=1}^m p_i^{c_i}\)，则：
证明：容斥原理。设每个 \(n\) 的质因子 \(p_i\) 在 \([1,n]\) 中的倍数集合为 \(A_i\)。则有：
又因为我们知道 \(|A_i|=\lfloor \frac n{p_i}\rfloor\)，于是就可以用容斥的公式：
我们还知道，\(|\bigcap_{i=1}^kA_{i}|=\lfloor\frac{n}{\prod_{i=1}^kp_i}\rfloor\)。于是：
性质二
欧拉函数是积性函数。
证明：（其中 \(\gcd(a,b)=1\)，在唯一分解定理中，\(a=\prod_{i=1}^np_i^{c_i}\)，\(b=\prod_{i=1}^mq_i^{k_i}\)）。
https://www.zhihu.com/question/410737433
性质三
证明（运用了莫比乌斯反演，可以先把整篇文章看完再回过头来看）：
upd: 我太 naive 了。现在告诉大家一种更简单的证明方法。
这个式子其实在算不超过 \(n\) 的与 \(n\) 互质的数的和。我们知道 \(\gcd(a,b)=\gcd(b,a-b)\)，所以对于每个 \(i\perp n\)，都有 \((n-i)\perp n\)，它们的和等于 \(n\)，所以事实上就是 \(\frac{\varphi(n)}{2}\) 个 \(n\) 相加。
埃氏筛与欧拉筛
埃氏筛
这种筛法虽然时间复杂度逊于欧拉筛，但是对于题目的一些特殊要求能够更好地应付。例如 筛一段区间内的质数 就只能用埃氏筛。
埃氏筛的代码如下：
for(ll i=2;i<=N;i++){
if(!vis[i]){
p.push_back(i);
if(i*i<=N){
for(ll j=i*i;j<=N;j+=i)vis[j]=1;
}
}
}
可以看出，埃氏筛是对于每个质数，把能整除它的所有和数全部标记一遍。但是一个小于 \(i^2\) 的合数必定有比 \(i\) 更小的质因子，所以它肯定已经被筛过，没有必要再筛。
如果想要筛一段区间 \([L,R]\) 内的质数，只需要用 \([1,\sqrt R]\) 内的质数去筛就可以了。
设需要筛出的质数范围为 \([1,N]\)，则其时间复杂度为 \(\mathcal O(N\log\log N)\)，证明需要用微积分，然而我不会，会微积分的可以取 OI Wiki 看证明。 现在会了
欧拉筛
现在我们想筛出 \(2\) 到 \(N\) 的所有质数。
欧拉筛的核心思想是从大到小累计质数，以达到每个合数只筛一次的目的。算法用一个数组来记录每个数的最小质因子。设这个数组为 \(v\)，我们可以用以下方式维护它：
- 从小到大用 \(i\) 把 \(2\) 到 \(N\) 的每个数扫一遍。
- 如果 \(v_i\) 没有被修改过，说明 \(i\) 没有被筛，所以 \(i\) 是质数，其最小质因子为其本身，即执行 \(v_i\gets i\)。
- 无论 \(i\) 是不是质数，都需要在 \(v_i\) 的基础上累计质因子。即扫描不大于 \(v_i\) 的每个质数 \(p\)，\(v_{pi}\gets p\)。
注意第三点，因为质因子是从大到小累积的，所以 \(p\) 一定是 \(pi\) 的最小质因子。
该算法从 \(2\) 到 \(N\) 的每个数都只会被遍历一次，筛一次，所以时间复杂度为 \(\mathcal O(N)\)。
P3383 线性筛质数模板的核心代码：
void init(int N){
for(int i=2;i<=N;i++){
if(!v[i]){v[i]=i;p.push_back(i);}
for(int pr:p){
if(pr>v[i]||pr>N/i)break;
v[pr*i]=pr;
}
}
}
线性筛欧拉函数
线性筛过程中，假设我们正在用一个数 \(n\) 和一个质数 \(p\) 筛 \(pn\)，且已知 \(\varphi(n)\)，需要求出 \(\varphi(pn)\)。设在唯一分解定理中，\(n=\prod_{i=1}^mq_i^{c_i}\)
- 当 \(p\mid n\) 时，\(\varphi(pn)=pn\prod_{i=1}^m(\frac{q_i-1}{q_i})=p\cdot\varphi(n)\)
- 当 \(p\nmid n\) 时，根据积性函数的定义，\(\varphi(pn)=\varphi(p)\cdot\varphi(n)\)
然后就能求出 \([2,N]\) 中的每一个整数的欧拉函数值了。只需在上面的程序中稍加修改：
inline void init(int n){
for(int i=2;i<=n;i++){
if(!v[i]){v[i]=i;phi[i]=i-1;p.push_back(i);}
for(int pr:p){
if(pr>v[i]||pr>n/i)break;
v[pr*i]=pr;
phi[pr*i]=phi[i]*(pr==v[i]?pr:phi[pr]);
}
}
}
线性筛任意积性函数
已知 \(f\) 是一个积性函数。设在唯一分解定理中，\(n=\prod_{i=1}^mp_i^{c_i}\)。则有：
因此，只要在欧拉筛时，统计每个数最小质因数 \(p\) 及其次数 \(c\)，就可以递推计算积性函数：
以后，在讨论积性函数时，只讨论其在 \(p^k\) 的取值然后相乘的方法，是比较常见的。
例题
P5481 [BJOI2015] 糖果
你可能会问这不是一道组合计数问题吗，为啥放到数论这里。
因为这道题的难点不在于计数，而在于任意模数。
计数式子插板法随便搞搞就得出来了，是
所以任意模数怎么办呢？如果直接 exgcd 算很可能出现逆元不存在的情况，而且还不一定是模数的倍数，所以答案不一定是 \(0\)。
然后我就忍不住看了题解，发现只要把 \(m!\) 分解质因数，然后把上面的式子除掉这些质因数就行了。
这里有一个分解阶乘质因数的技巧：枚举所有小于等于 \(m\) 的质数 \(p\)，然后计算 \(p\) 在 \(m!\) 中出现的次数。这个次数就等于能整除 \(p\) 的数的数量加上能整除 \(p^2\) 的数的数量加上能整除 \(p^3\) 的数的数量，等等。这样一来，如果一个数含有 \(k\) 个 \(p\)，其贡献就为 \(k\)，所以这个算法是正确的。故 \(p\) 在 \(m!\) 中出现的的次数为：
代码如下：
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1e5+5;
int pr[N],c[N],d[N],tot;int P,v[N];
inline void init(int m){
const int n=N-5;
for(int i=2;i<=n;i++){
if(!v[i])pr[++tot]=v[i]=i;
for(int j=1;j<=tot;j++){
if(pr[j]>v[i]||pr[j]>n/i)break;
v[pr[j]*i]=pr[j];
}
}
for(int i=1;i<=tot;i++){
if(m<i)break;
ll x=pr[i];
while(m>=x){
c[i]+=m/x;
x*=pr[i];
}
}
}
inline ll dow(ll a,int b){
ll ret=1;
while(b--){
ret=(ret*a)%P;
a--;
}
return ret;
}
int n,m,K;
int main(){
scanf("%d%d%d%d",&n,&m,&K,&P);
init(m);ll ans=1;
for(int i=K;i<K+m;i++)d[i-K]=i;
for(int i=1,p;i<=tot;i++){
p=pr[i];if(!c[i])break;
for(int j=(K%p?(K/p+1)*p:K);j<K+m;j+=p){
while(c[i]&&d[j-K]%p==0)d[j-K]/=p,c[i]--;
}
}
for(int i=K;i<K+m;i++)ans=(ans*d[i-K])%P;
printf("%lld\n",dow(ans,n));
return 0;
}
P2158 仪仗队
一个人 \((x,y)\) 能被看到，当且仅当 \(x+y=1\) 或 \(\gcd(x,y)=1\)（这里下标从 0 开始）。然后因为正方形是沿对称轴对称的，所以我们可以先算出答案的不严格一半，即 \(x\ge y\) 的情况。
于是，问题转化为：求 \(\text{ans}=\sum_{x=1}^N\sum_{y=1}^x[\gcd(x,y)=1]=\sum_{x=1}^n\varphi(x)\)。
线性筛一下就好了，答案为 \(2(\text{ans}-1)+3\)，其中 \(-1\) 是减去 \((1,1)\) 的情况（因为这个点在对称轴上）\(+3\) 是加上 \((0,1),(1,0),(1,1)\) 的情况。
细节：注意 \(\varphi(1)=1\)，需要提前赋值。注意输入要减一，而且要特判输入等于一的情况。
该算法的时间复杂度为 \(\mathcal O(n)\)。
我的代码：
#include<cstdio>
#include<vector>
using namespace std;
const int N=4e4+5;
int v[N],phi[N];
long long ans;
vector<int>p;
inline void init(int n){
ans=phi[1]=1;
for(int i=2;i<=n;i++){
if(!v[i]){v[i]=i;phi[i]=i-1;p.push_back(i);}
for(int pr:p){
if(pr>v[i]||pr>n/i)break;
v[pr*i]=pr;
phi[pr*i]=phi[i]*(pr==v[i]?pr:phi[pr]);
}
ans+=phi[i];
}
}
int n;
int main(){
scanf("%d",&n);
init(n-1);
printf("%lld",(n-1)?((ans-1)<<1)+3:0ll);
return 0;
}
同余
若 \(a\bmod m=b\bmod m\)，则称 \(a\equiv b\pmod m\)。
同余类和剩余系
定义
对于 \(a\in[0,m-1]\)，将集合 \(\{x\in\mathbb{N}\mid x\equiv a\pmod m\}=\{a+km\mid k\in\mathbb{N}\}\) 称为模 \(m\) 的一个同余类，简记为 \(\overline a\)。
\(m\) 个同余类构成 \(m\) 的完全剩余系（即 \(\{\overline 0,\overline 1,\cdots,\overline{m-1}\}\)）。
\(1\) 到 \(m\) 中与 \(m\) 互质的数所代表的剩余系有 \(\varphi(m)\) 个，它们构成 \(m\) 的简化剩余系。
性质
简化剩余系关于模 \(m\) 乘法封闭。
证明
设 \(a,b\) 是 \(m\) 简化剩余系中的两个数，\(a=k_1m+x,b=k_2m+y\ (1\le x,y<m)\)。
则 \(ab=(k_1m+x)(k_2m+y)=k_1k_2m^2+k_2mx+k_1my+xy\)，所以 \(ab\equiv xy\pmod m\)
根据定义，\(\gcd(x,m)=\gcd(y,m)=1\)，所以 \(\gcd(xy,m)=1\)。
根据欧几里得算法公式，\(\gcd(xy,m)=\gcd(m,xy\bmod m)=1\)，则简化剩余系关于模 \(m\) 乘法封闭。
证毕。
欧拉定理和费马小定理
内容
欧拉定理：若正整数 \(a,n\) 互质，则 \(a^{\varphi(n)}\equiv 1\pmod n\)。
费马小定理：若 \(p\) 是质数，则对于任意整数 \(a\)，有 \(a^p\equiv a\pmod n\)。
证明
设 \(n\) 的简化剩余系为 \(\{\overline{a_1},\overline{a_2},\cdots,\overline{a_{\varphi(n)}}\}\)。根据定义，\(1\le a_i<n\)。
\(\because\gcd(a,n)=1,\gcd(a_i,n)=1\)
\(\therefore\gcd(aa_i,n)=1\) 即 \(\gcd(n,aa_i\bmod n)=1\)
\(\therefore\forall a_i\)，\(\overline{aa_i\bmod n}\) 都在 \(n\) 的简化剩余系中。
设 \(\overline{aa_i\bmod n}\) 和 \(\overline{aa_j\bmod n}\) 是相同的同余类，则 \(aa_i\equiv aa_j\pmod n\)，即 \(a(a_i-a_j)\equiv 0\pmod n\)
\(\therefore n\mid a(a_i-a_j)\) 即 \(n\mid (a_i-a_j)\)
\(\therefore a_i\equiv a_j\pmod n\)
\(\therefore a_i=a_j\)
于是我们得到一个成立的命题：『如果 \(\overline{aa_i\bmod n}\) 和 \(\overline{aa_j\bmod n}\) 是相同的同余类，那么\(a_i=a_j\)』。因此它的逆否命题『如果 \(a_i\ne a_j\)，那么 \(\overline{aa_i\bmod n}\) 和 \(\overline{aa_j\bmod n}\) 不是相同的同余类』也成立。
\(\therefore\{\overline{aa_1\bmod n},\overline{aa_2\bmod n},\cdots,\overline{aa_{\varphi(n)}\bmod n}\}=\{\overline{a_1},\overline{a_2},\cdots,\overline{a_{\varphi(n)}}\}\)
于是：
\(\therefore a^{\varphi(n)}\equiv 1\pmod n\)
至于费马小定理，其实就是欧拉定理的一个特例（除了 \(p\mid n\) 的情况，不过这种情况下费马小定理显然成立）。
欧拉定理的推论
内容
证明
扩展欧拉定理
例题
P5091 【模板】扩展欧拉定理
套公式即可。代码是很久以前写的，所以码风可能和现在不一样。
#include<cctype>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
inline ll read(){
char c=getchar();ll x=0;bool f=0;
for(;!isdigit(c);c=getchar())f^=!(c^45);
for(;isdigit(c);c=getchar())x=(x<<1)+(x<<3)+(c^48);
return f?-x:x;
}
inline ll getB(ll mod){
char c=getchar();ll b=0;bool f=0;
for(;!isdigit(c);c=getchar());
for(;isdigit(c);c=getchar()){
b=(b<<1)+(b<<3)+(c^48);
if(b>mod)b%=mod,f=1;
}
return f?b+mod:b;
}
inline ll phi(ll x){
ll res=x;
for(register ll i=2;i*i<=x;i++){
if(x%i==0){
res=res*(i-1)/i;
while(x%i==0)x/=i;
}
}
if(x>1)res=res*(x-1)/x;
return res;
}
inline ll quickPower(ll a,ll b,ll n){
ll res=1;
while(b){
if(b&1){
res=(res*a)%n;
}
a=(a*a)%n;
b>>=1;
}
return res;
}
ll a,m,phim,b;
int main(){
a=read();m=read();
phim=phi(m);
b=getB(phim);
printf("%lld",quickPower(a,b,m));
return 0;
}
数论分块
过程
对于形如 \(\sum_{i=1}^nf(i)g(\lfloor\frac ni\rfloor)\) 的式子，由于 \(\left\lfloor\dfrac ni\right\rfloor\) 的取值至多只有 \(2\sqrt n\) 种，因此可以把和式分块计算，而 \(f(i)\) 使用前缀和处理。
结论：值 \(\lfloor\frac ni\rfloor\) 所在块的右端点为
我们也可以类似地搞一个二维（或多维）数论分块，用于求解类似 \(\sum_{i=1}^nf(i)g(\lfloor\frac ni\rfloor)h(\lfloor\frac mi\rfloor)\) 的式子。这时我们应该取右端点为所有维度右端点的最小值。例如二维数论分块应为 r=min(n/(n/i),m/(m/i))。
显然，二维数论分块的时间复杂度仍然是 \(\mathcal O(\sqrt n)\) 的（假设 \(n\) 和 \(m\) 同级）。
证明
令 \(k=\lfloor\dfrac ni\rfloor\)。
例题
UVA11526
题意
求 \(\sum_{i=1}^n\lfloor\dfrac ni\rfloor\)。\(T\) 组数据。
分析
模板题。按照上面的流程即可。注意边界情况和运算顺序。
时间复杂度为 \(\mathcal O(T\sqrt n)\)。
代码
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
ll H(ll n){
ll ret=0;
for(ll l=1,r=0;l<=n;l=r+1){
r=n/(n/l);ret+=(r-l+1)*(n/l);
}
return ret;
}
int main(){
ll T,n;
scanf("%lld",&T);
while(T--){
scanf("%lld",&n);
printf("%lld\n",H(n));
}
return 0;
}
狄利克雷卷积
从现在开始，数论函数一律使用粗黑体（如 \(\mathbf f,\mathbf g\)）或希腊字母（如 \(\varepsilon\)）。
定义两个数论函数的加法为 \((\mathbf{f+g})(n)=\mathbf f(n)+\mathbf g(n)\)，点乘为 \((\mathbf{f\cdot g})(n)=\mathbf{f}(n)\mathbf{g}(n)\)。下面将省略点乘符号。
当 \(\mathbf f\) 是完全积性函数时，\(\mathbf{(fg)*(fh)=f(g*h)}\)。
定义
运算
为狄利克雷卷积。也可以写成：
性质
有单位元 \(\varepsilon(n)=[n=1]\)，即 \(\mathbf f*\varepsilon=\mathbf f\)。
当 \(\mathbf{f,g}\) 都是积性函数时，\(\mathbf{f*g}\) 也是积性函数。
一个积性函数的逆也是积性函数。
对于任意 \(\mathbf f(1)\ne 0\) 的函数 \(\mathbf f\)，存在逆元 \(\mathbf g\) 使得 \(\mathbf{f*g}=\varepsilon\)。
若 \(\mathbf g\) 是 \(\mathbf f\) 的逆元，只需要满足以下式子即可：
证明
挑一些不那么明显的。
结合律：
这就与三个函数的顺序无关了。
求 \(\mathbf f\) 的逆元 \(\mathbf g\)：
其他先鸽了。
一些常见的数论函数及其性质
定义
定义 \(\mu\) 为 \(\mathbf 1\) 的逆元。即 \(\mu*\mathbf 1=\varepsilon\)。
性质
定义中提到的所有函数均为积性函数。其中前两个是完全积性函数。
其实只要记住前三个，后面的都可以自己推出来。
这个式子自己代入定义就很容易得到了。
证明
前两个都很显然。这里证明第三个。设 \(p\) 为任意质数。
由于 \(\mathbf 1*\varphi\) 是积性函数，所以有：
倒数第二个可以用前面提到的欧拉函数的性质证明。以后做题推式子要注意不要忘记前面提到的公式。
最后一个可以通过分类讨论质数在 \(x,y\) 中的分布情况证明。
莫比乌斯反演
定理
若 \(\mathbf g=\mathbf{f*1}\)，则有 \(\mathbf f=\mathbf{g}*\mu\)。
当然，一般的题目肯定不会直接写出卷积的形式。
所以我们这么写：
魔力筛
虽然在网上没看见过这种叫法，倒是搜到了魔力筛让天价药物降价，但机房都这么叫，我也暂且这么叫吧。
这是一个可以在 \(\mathcal O(n\log\log n)\) 的时间复杂度求出 任意数论函数 与 \(\mu\) 的狄利克雷卷积的算法。当然，必须保证这个数论函数能被提前计算出来。
设 \(\mathbf g_{i,n}=\sum_{d\mid n}\mathbf f(d)\mu(\dfrac nd)\)，其中 \(d\) 只含前 \(i\) 种质因子。则有：
需要滚动数组优化。
虽然可以被解释为高维差分，但我不会。
其实这个 DP 方程也比较好理解。第一种情况显然正确。我们来看看第二种：
其中 \(d\) 只包含前 \(i\) 种质因子，\(x,y\) 只包含前 \(i-1\) 种质因子。我们知道，每多一个质因子，\(\mu\) 的取值就会乘以 \(-1\)，因此上式的意义是：强制不选 \(p_i\)，然后乘以 \(-1\) 累加在答案上。如果含有平方因子，则值为 \(0\)，不影响答案。
复杂度和埃氏筛一样。
参考代码：
inline void calc(int n){
for(int i=1;i<=n;i++)g[i]=f[i];
for(int p:pr){
for(int i=n/p;i>=1;i--){
g[i*p]=g[i*p]-g[i];
}
}
}
应用
练习题
题目描述
给定 \(n,m\le 10^6\)，求
题目分析
我们知道 \(d\mid\gcd(a,b)\iff d\mid a\land d\mid b\)。
有一个常见的套路：若要求 \(\mathbf f(\gcd(x,y))\)，可以先找到一个数论函数 \(\mathbf g=\mathbf{1*f}\)，这样子就有：
于是解决这道题就变得很简单了。我们知道 \(\mathbf{id}=\mathbf 1*\varphi\)。所以好像并不需要莫比乌斯反演：
然后就可以用二维数论分块搞了。不过现在应该不会出这么套路的题了吧。
[国家集训队]Crash的数字表格 / JZPTAB
题目描述
给定 \(n,m\)，求
对 \(20101009\) 取模的结果。
数据范围：\(1\le n,m\le 10^7\)。
题目分析
令 \(\mathbf{id}^{-1}=\mathbf{g*1}\)，则 \(\mathbf g=\mathbf{id}^{-1}*\mu\)。
\(\mathbf g\) 显然是积性函数，并且在质数幂处很容易求。因此我们可以线性计算之。
时间复杂度为 \(\mathcal O(n)\)。数论分块好像没有必要，如果不用数论分块甚至可能会快一点
代码
#include<cstdio>
#include<vector>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1e7+5,P=20101009;
int v[N],t[N];
vector<int>pr;
ll inv[N],g[N],s[N];
inline void init(int n){
inv[1]=1;g[1]=1;s[1]=1;
for(int i=2;i<=n;i++){
inv[i]=(P-P/i)*inv[P%i]%P;
if(!v[i]){
g[i]=(inv[i]+P-1)%P;
pr.push_back(t[i]=v[i]=i);
}else{
if(i==t[i])g[i]=(inv[i]-inv[i/v[i]]+P)%P;
else g[i]=g[t[i]]*g[i/t[i]]%P;
}
for(int p:pr){
if(p>v[i]||p>n/i)break;
v[i*p]=p;t[i*p]=(v[i]==p?t[i]*p:p);
}
s[i]=(s[i-1]+g[i]*i%P*i%P)%P;
}
}
inline ll calc(int n,int m){
ll ret=0;
for(int l=1,r;l<=n;l=r+1){
r=min({n,n/(n/l),m/(m/l)});
ret=(ret+(s[r]-s[l-1]+P)%P*(n/l)%P*(n/l+1)%P*(m/l)%P*(m/l+1)%P)%P;
}
return ret*inv[4]%P;
}
int main(){
int n,m;scanf("%d%d",&n,&m);
if(n>m)swap(n,m);init(n);
printf("%lld\n",calc(n,m));
return 0;
}
[SDOI2017]数字表格
题目描述
给定 \(n,m\)，求
对 \(10^9+7\) 取模的结果。其中 \(F_i\) 表示斐波那契数列的第 \(i\) 项。
有 \(T\) 组数据。数据范围：\(1\le T\leq 10^3\)，\(1\le n,m\le 10^6\)。
题目分析
我们定义 \(\mathbf f(n)=\ln F_n\)，\(\mathbf{g*1}=\mathbf f\)。则：
这个式子可以外层数论分块，内层暴力枚举计算。
单次询问的时间复杂度为 \(\mathcal O(\sqrt n)\)。注意对 \(\mu\) 的处理。
另一种方法是使用类似魔力筛的方法计算内层，从而使预处理的时间复杂度降到 \(\mathcal O(n\log\log n)\)，请自行参阅题解。
代码
#include<cstdio>
#include<vector>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1e6+5,P=1e9+7;
inline ll _pow(ll a,ll b){
ll ret=1;a%=P;
while(b){if(b&1)ret=(ret*a)%P;a=a*a%P;b>>=1;}
return ret;
}
int mu[N],v[N];
ll F[N],G[N];
vector<int>pr;
inline void init(){
mu[1]=1;F[0]=0;F[1]=1;G[0]=G[1]=1;
const int n=N-5;
for(int i=2;i<=n;i++){
G[i]=1;F[i]=(F[i-1]+F[i-2])%P;
if(!v[i])pr.push_back(v[i]=i),mu[i]=-1;
for(int p:pr){
if(p>v[i]||p>n/i)break;
v[i*p]=p;
if(v[i]==p)mu[i*p]=0;
else mu[i*p]=mu[i]*mu[p];
}
}
ll mF;
for(int i=2;i<=n;i++){
mF=_pow(F[i],P-2);
for(int j=i;j<=n;j+=i){
if(mu[j/i]==-1)G[j]=G[j]*mF%P;
else if(mu[j/i]==1)G[j]=G[j]*F[i]%P;
}
}
for(int i=2;i<=n;i++){
G[i]=G[i-1]*G[i]%P;
}
}
int main(){
int T,n,m;scanf("%d",&T);
ll ans;init();
while(T--){
scanf("%d%d",&n,&m);
if(n>m)swap(n,m);
ans=1;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=ans*_pow(G[r]*_pow(G[l-1],P-2)%P,(ll)(n/l)*(m/l))%P;
}
printf("%lld\n",ans);
}
return 0;
}
YY 的 GCD
题目大意
设 \(P\) 为所有质数构成的集合。求
\(T\) 组数据。数据范围：\(T=10^4\)，\(n,m\le 10^7\)。TL：4s。
题目分析
设 \(\mathbf f(n)=[n\in P]\)，\(\mathbf{g*1}=\mathbf f\)。
然后就和上面一样了...... 直接写结果吧。
由于 \(\mathbf f\) 不是积性函数，我们需要用魔力筛计算 \(\mathbf g\)。然后就可以数论分块搞了。
时间复杂度为 \(\mathcal O(n\log\log n+T\sqrt n)\)。
代码
#include<cstdio>
#include<vector>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1e7+5;
int v[N],g[N];
vector<int>pr;
inline void init(){
const int n=N-5;
for(int i=2;i<=n;i++){
if(!v[i])pr.push_back(v[i]=i),g[i]=1;
for(int p:pr){
if(p>v[i]||p>n/i)break;
v[p*i]=p;
}
}
for(int p:pr){
for(int i=n/p;i>=1;i--)g[i*p]-=g[i];
}
for(int i=1;i<=n;i++)g[i]+=g[i-1];
}
int main(){
init();
int n,m,T;scanf("%d",&T);
ll ans;
while(T--){
scanf("%d%d",&n,&m);
if(n>m)swap(n,m);ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min({n,n/(n/l),m/(m/l)});
ans+=(ll)(g[r]-g[l-1])*(n/l)*(m/l);
}
printf("%lld\n",ans);
}
return 0;
}
浙公网安备 33010602011771号