数论算法总结
数论算法
快速乘
ll mul(ll a,ll b,ll mod)
{
ll tmp=a*b-(ll)((long double)a/mod*b+0.5)*mod;
return tmp<0?tmp+mod:tmp;
}
同余系相关
Exgcd
存在 $ax+by=\gcd(a,b)因为 $所以 \(a x_{1}+b y_{1}=b x_{2}+\left(a-\left(\left\lfloor\frac{a}{b}\right\rfloor \times b\right)\right) y_{2}\)
\(a y_{2}+b x_{2}-\left\lfloor\frac{a}{b}\right\rfloor \times b y_{2}=a y_{2}+b\left(x_{2}-\left\lfloor\frac{a}{b}\right\rfloor y_{2}\right)\)
所以 \(x_{1}=y_{2}, y_{1}=x_{2}-\left\lfloor\frac{a}{b}\right\rfloor y_{2}\)
递归求解，终止时必然有 \(x=1,y=0\) 。
解集为 \(\{(x,y)|x=x_1+k\frac{b}{\gcd(a,b)},y=y_1-k\frac{a}{\gcd(a,b)}\}\)
void exgcd(ll a,ll b,ll &x,ll &y)
{
if(b==0){x=1,y=0;return;}
exgcd(b,a%b,x,y);
ll t=y;y=x-y*(a/b),x=t;
}
ExCRT
每次合并两个方程 \(x \equiv b_{1}\left(\bmod a_{1}\right)\) 和 \(x \equiv b_2\left(\bmod a_2\right)\)
存在 \(x=b_1+a_1k_1,x=b_2+a_2k_2\to b_1+a_1k_1=b_2+a_2k_2\)
求解二元一次方程 \(a_1k_1-a_2k_2=b_2-b_1\)
设 \(c\equiv k_1(\bmod \frac{a_2}{\gcd(a_1,a_2)})\)
可得新的方程为
\(x \equiv ca_{1}+b_{1}(\bmod \frac{a_1a_2}{\gcd(a_1,a_2)})\)
ll ExCRT(ll *a,ll *p,const int &n)
{
ll A=a[1],M=p[1],k1,k2;
for(int i=2;i<=n;i++)
{
ll c=((a[i]-A)%p[i]+p[i])%p[i],g=gcd(M,p[i]);
if(c%g)return -1;
exgcd(M,p[i],k1,k2);
k1=mul(k1,c/g,p[i]);
ll now=M/g*p[i];
A=(mul(k1,M,now)+A)%now,M=now;
}
return (A%M+M)%M;
}
ExBSGS
根据欧拉定理 \(a^{\phi(p)}\equiv 1(\bmod p)\) ，\(p\) 与 \(a\) 互质。
可以在 \(O(p)\) 的时间复杂度解决这个问题。
但是对于一般情况 \(p\) 不一定与 \(a\) 互质。
设 \(g=\gcd(a,p)\)
若 \(g|b\) ，则存在
否则，若 \(b=1\) 解为 \(x=0\)，\(b≠1\) 无解
由 $\frac{a}{g} $ 与 $ \frac{p}{g}$ 互质可得
继续分解直到 \(g=1\) 为止
这时候需要特判 \(a^1,a^2,...,a^t\) 是否存在与 \(b\) 同余的数。
因为 \(g=1\) ，所以原问题就被转换为了 \(a^x\equiv b(\bmod{p})\)，且 \(a,p\) 互质。
考虑分块求解这个问题，设 \(n=\lceil\sqrt{q}\rceil\)
显然可以使用 \(xn-y,x,y\leq n\) 来表示一个小于 \(p\) 的数
\(a^{xn}\equiv ba^{y}(\bmod{p})\)
只需要将所有的 \(ba^{y}\) 存入Hash表中，再对所有 \(y\) 查询合法的值。
int ExBSGS(int a,int b,int p)
{
if(a==0&&b==0)return 1;
if(a==0)return -1;
if(b==1)return 0;
unordered_map<int,int>vis;
int g=gcd(a,p),mag=1,t=0;
while(g!=1)
{
if(b%g)return -1;
++t,p/=g,b/=g,mag=mag*1LL*(a/g)%p;
if(mag==b)return t;
g=gcd(a,p);
}
vis.clear();
int n=sqrt(p)+1;
for(int i=0,m=b;i<n;i++,m=m*1LL*a%p)vis[m]=i;
a=Pow(a,n,p);
for(int i=1,m=mag*1LL*a%p;i<=n;i++,m=m*1LL*a%p)if(vis.count(m))return i*n-vis[m]+t;
return -1;
}
Lucas
首先可以观察得出 \(C_p^i=\frac{p}{i}C_{p-1}^{i-1}\) ，所以显然有 \(C_p^i\equiv0(\bmod p)\)
设 \(n=a_0+a_1p+a_2p^2+a_3^p+...+a_kp^k,m=b_0+b_1p+b_2p^2+...+b_kp^k\)
lucas定理也可以表示为一下形式
\((1+x)^n\) 中 \(x^m\) 的系数显然为 \(C_m^n\)
可得 \(x^m\) 的系数为 $\prod_{i=0}kC_{a_i} $
因为 \(m=b_0+b_1p+b_2p^2+...+b_kp^k\) 若 \(\exists a_i<b_i\) 则 \(m\) 无法被表示系数为 \(0\)
得证。
int lucas(int n,int m,int mod){if(n==0&&m==0)return 1;return lucas(n/mod,m/mod,mod)*1LL*C(n%mod,m%mod,mod)%mod;}
Exlucas
这里的 \(p\) 可能不是质数
根据唯一分解定理 \(p=p_1^{k_1}p_2^{k_2}...p_t^{k_t}\)
若 \(\forall k_i,k_i=1\) 可以直接对于每个质数使用lucas求解，再使用CRT合并。
否则需要求解子问题 \(C_n^m \bmod p_i^{k_i}\) 再使用CRT合并。
只需要解决 \(\frac{n!}{p^a}\bmod p^k\) ，其中 \(a\) 是 \(n!\) 包含的 \(p\) 的次数。
考虑将 \(n!\) 分组，\(n!=(p\times2p\times 3p\times 4p\times...\times\lfloor\frac{n}{p}\rfloor p)\times1\times2\times...\times(p-1)\times(p+1)\times...\times(2p-1)\times...\)
可以发现后面是模 \(p^k\) 是循环的，可以通过预处理求出
前面的 \(\lfloor\frac{n}{p}\rfloor!\) 是一个子问题递归求解即可。
复杂度 \(O(p\log n)\)
ll mul(ll a,ll b,ll mod)
{
ll tmp=a*b-(ll)((long double)a/mod*b+0.5)*mod;
return tmp<0?tmp+mod:tmp;
}
ll Pow(ll a,ll k,ll mod)
{
ll ret=1;
while(k){if(k&1)ret=ret*1LL*a%mod;k>>=1,a=a*1LL*a%mod;}
return ret;
}
ll fact(ll n,ll p,ll pk)
{
if(!n)return 1;
ll ans=1;
for(ll i=1;i<pk;i++)if(i%p)ans=ans*1LL*i%pk;
ans=Pow(ans,n/pk,pk);
for(ll i=1;i<=n%pk;i++)if(i%p)ans=ans*1LL*i%pk;
return ans*1LL*fact(n/p,p,pk)%pk;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
void exgcd(ll a,ll b,ll &x,ll &y)
{
if(b==0)return x=1,y=0,void();
exgcd(b,a%b,x,y);
ll t=x;x=y,y=t-y*(a/b);
}
ll inv(ll x,ll p)
{
ll t1,t2;
exgcd(x,p,t1,t2);
return (t1%p+p)%p;
}
ll getpower(ll x,ll p)
{
ll c=0;
while(x)c+=x/p,x/=p;
return c;
}
ll C(ll n,ll m,ll p,ll pk)
{
ll pa=getpower(n,p)-getpower(m,p)-getpower(n-m,p);
return Pow(p,pa,pk)*1LL*fact(n,p,pk)%pk*inv(fact(m,p,pk),pk)%pk*inv(fact(n-m,p,pk),pk)%pk;
}
ll ExCRT(ll *a,ll *p,const ll&n)
{
ll A=a[1],M=p[1],k1,k2;
for(ll i=2;i<=n;i++)
{
ll c=((a[i]-A)%p[i]+p[i])%p[i],g=gcd(M,p[i]);
if(c%g)return -1;
exgcd(M,p[i],k1,k2);
k1=mul(k1,c/g,p[i]);
ll now=M/g*p[i];
A=(mul(k1,M,now)+A)%now,M=now;
}
return (A%M+M)%M;
}
ll a[30],b[30];
ll cnt;
ll exlucas(ll n,ll m,ll p)
{
cnt=0;
ll lim=sqrt(p);
for(ll i=2;i<lim;i++)
{
if(p%i)continue;
ll t=1;while(p%i==0)p/=i,t*=i;
a[++cnt]=t,b[cnt]=C(n,m,i,t);
}
if(p>1)a[++cnt]=p,b[cnt]=C(n,m,p,p);
return ExCRT(b,a,cnt);
}
二次剩余
无特殊说明 \(p\) 均为奇素数。
对于 \(p,n\) ，若存在 \(x\) ，满足
则称 \(n\) 为模 \(p\) 意义下的二次剩余
勒让德符号
欧拉判别准则
证明：
若 \(n\) 为模 \(p\) 意义下的二次剩余，即存在 \(x^2\equiv n(\bmod p)\) ，由费马小定理可得 \(x^{p-1}\equiv 1(\bmod p)\)，显然 \(n^{\frac{p-1} {2}}\equiv 1(\bmod p)\) 。
若 \(n\) 为模 \(p\) 意义下的非二次剩余，假设存在 \(x^2\equiv n(\bmod p)\) ，此时 \(x^{p-1}\equiv -1(\bmod p)\) ，由费马小定理可知 \(x\) 不存在
若 $n\equiv 0(\bmod p) $ ，显然成立。
Cipolla
定理1
定理2
\(p\) 的二次剩余与非二次剩余的个数均为 \(\frac{p-1}{2}\) （不考虑 \(0\) 的情况下）。
证明：
我们只考虑所有的 \(n^2\)，假设有 \(x≠y\) 且 \(x^2≡y^2(\bmod p)\)，则 \(p∣(x^2−y^2)\)，即 \(p∣(x−y)(x+y)\)，若\(p∤(x−y)\)，则 \(p∣(x+y)\)，故 \(x+y≡0(\bmod p)\)，也就是不同的 \(x^2\) 共有 $\frac{p-1}{2} $ 个，即二次剩余有个 \(\frac{p-1}{2}\) 个。
算法流程
- 判断给定的数 \(x\) 是否是二次剩余。
- 随机一个 \(a\)，使其满足 \((a^2−x)\) 是二次非剩余。
- 令 \(\omega^2\equiv (a^2-x)\pmod{p}\) ，取 \(y\equiv \left(a+{\omega}\right)^{p+1\over 2}\) 为其中一个可行解。
证明：
int mod,w;
struct Complex
{
int a,b;
Complex(int A=0,int B=0){a=A,b=B;}
Complex operator * (const Complex &t)const
{
return Complex((a*1LL*t.a+b*1LL*t.b%mod*w)%mod,(a*1LL*t.b+b*1LL*t.a)%mod);
}
};
int Pow(int a,long long k)
{
int ret=1;
while(k){if(k&1)ret=ret*1LL*a%mod;k>>=1,a=a*1LL*a%mod;}
return ret;
}
int Pow(Complex a,long long k)
{
Complex ret=Complex(1,0);
while(k){if(k&1)ret=ret*a;k>>=1,a=a*a;}
return ret.a;
}
int Sqrt(int x)
{
if(!x)return 0;
if(Pow(x,(mod-1)>>1)==mod-1)return -1;
while(true)
{
int a=rand()*1LL*rand()%mod;
w=(a*1LL*a%mod+mod-x)%mod;
if(Pow(w,(mod-1)>>1)==mod-1)return Pow(Complex(a,1),(mod+1)>>1);
}
}
素数分解相关
Miller-Rabin
判断 \(q\) 是否为素数（\(q\leq 10^{18}\)）
根据费马小定理，当 \(a \bmod q ≠0\) 且 \(q\) 为质数时，\(a^{q-1}≡1 (\bmod q)\)
而当 \(q\) 不是质数时，不一定成立
因此就有了一个想法：随机一些 \(a\) ，对每一个 \(a\) 验证 ，如果有一个错误那么它一定是合数
但是，存在这样一类合数 \(q\) ，对于 \(1≤a<q\) 上面的式子都成立，比如 \(561\)
二次探测定理
如果 $x^2≡1(\bmod p) $，那么 \(x≡1(\bmod p)\) 或者 \(x≡p-1(\bmod p)\) ，这里 \(p\) 是质数。
可以发现对于合数，有一定概率不成立。
在检验时，如果 \(a^{q-1}≡1(\bmod q)\) 且 \(2|(q-1)\) ，则下一步计算 \(a ^{\frac{q-1}{2}} \bmod q\) ,如果算出来不是 \(1\) 或 \(q-1\) ，那么 \(q\) 是合数，如果算出来是 \(1\) 并且 \(2|\frac{q-1}{2}\) ，则继续计算 \(\frac{q-1}{4}\) ，直到出现奇数或者算出 \(q-1\) 。
ll mul(ll a,ll b,ll mod)
{
ll tmp=a*b-(ll)((long double)a/mod*b+0.5)*mod;
return tmp<0?tmp+mod:tmp;
}
ll Pow(ll a,ll k,ll mod)
{
ll ret=1;
while(k){if(k&1)ret=mul(ret,a,mod);k>>=1,a=mul(a,a,mod);}
return ret;
}
const int tprim[10]={0,2,3,5,7,11,13,17,19,23};
bool Miller_Rabin(ll x)
{
if(x<2)return 0;
ll k=x-1,cnt=0;
while(!(k&1))k>>=1,++cnt;
for(int i=1,j;i<=9;i++)
{
if(tprim[i]==x)return 1;
if(Pow(tprim[i],x-1,x)!=1)return 0;
ll now=Pow(tprim[i],k,x);
if(now==1||now==x-1)continue;
now=mul(now,now,x);
for(j=1;j<=cnt;++j,now=mul(now,now,x))if(now==x-1)break;
if(j>cnt)return 0;
}
return 1;
}
递归改成递推，这样复杂度是 \(O(k\log n)\)
莫比乌斯反演
莫比乌斯反演函数
性质
莫比乌斯反演定理
若函数 \(F(n)\) 满足
则存在
证明：
考虑 \(f(i)\) 的系数可以很简单地得到证明。
推论：
若函数 \(F(d)\) 满足
则存在
证明：
考虑每个 \(f(k\cdot d)\) 的系数被贡献时满足 \(i\cdot j=k\)，显然 \(i\) 是 \(k\) 的约数，被贡献的系数值为 \(\mu(i)\)。
这个推论可以很简单地用于求解类似 \(\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]\) 问题
一些例题
基础的套路不在过多阐述。
Problem 1.
设 \(f(d)\) 满足 （满足 \(\gcd(i,j)=d\) 的 \(\frac{ij}{d}\) 的和）
显然
由莫比乌斯反演的性质可得
枚举 \(k\)
这里可以直接数论分块 \(O(n)\) 解决。
继续将这个式子变形， 考虑枚举 \(T=kd\)
可以发现后面的 \(g(T)=\sum_{d|T}\mu(d)d\) 是一个积性函数，可以线性筛出来。
这样对于单次询问复杂度 \(O(\sqrt{n})\)
Problem 2.
\(d(i)\) 为 \(i\) 的约数个数
有一个结论
显然每个质因子 \(p\) 的贡献是独立的，\(d(x)\) 为每种质因子贡献的乘积，设 \(i=i'p^{k_1},j=j'p^{k_2}\) 可以发现左边贡献为 \(k_1+k_2+1\) ，右边贡献为 \(k_1+1+k_2+1-1\) （\([\gcd(p^{t_1},p^{t_2})=1],t_1\leq k_1,t_2\leq k_2\)），所以得证。
按照套路考虑枚举 \(d\) ，显然此时 \(x,y\) 对 \(d\) 有贡献，仅当 \(\gcd(x,y)\) 是 \(d\) 的倍数。
所以 \(x\) 为 \(d\) 的倍数，\(y\) 为 \(d\) 的倍数。
考虑设 \(f(n)\) 为
显然 $f(n) $ 不难发现 \(f(n)=\sum_{i=1}^n d(i)\)
考虑线性筛筛出 \(d(i)\)
这个显然可以数论分块，总复杂度 $O(n) $
总结
一般来说莫比乌斯反演的题目，主要考虑化为 \([\gcd(x,y)=1]\) ，通过交换枚举顺序构造积性函数，或者使用数论分块求解，最重要的是抓住运算的性质和熟练掌握反演公式。
Min_25 筛
规定 \(P\) 是所有质数组成的集合，若无特殊说明 \(p\) 为质数。
其中 \(f(n)\) 为积性函数，并且 \(f(p^k)\) 能够快速计算，\(f(p)\) 能够表示为一个简单多项式。
筛质数的答案
首先需要对每个 \(n=\lfloor\frac{N}{i}\rfloor\)求出 \(\sum_{i=1}^n[i\in P]f(i)\)
首先线性筛出 \(\sqrt{n}\) 以内的所有质数，\(P_i\) 表示第 \(i\) 个质数。
因为 \(f(p)\) 能够写成一个简单多项式
所以本质上来说需要求解的是 \(\sum_{i=1}^n [i\in P] i^k\)
考虑构造一个函数 \(g(n,j)\) 满足
$\min (i) $ 表示 \(i\) 的最小质因子。
其实不难发现 \(g(n,j)\) 相当于运行了 \(j\) 次埃氏筛后，没有被筛掉的数和质数的 \(i^k\) 的和。
考虑如何求 \(g(n,j)\) ，发现这个转移可以分类讨论
考虑第 \(j\) 次筛质数，显然筛掉的是最小质因子大于 \(P_j\) 的数，而满足这个条件数必然大于 \(P_j^2\)
所以如果 \(n<P_j^2\) ，这次筛是不会筛掉任何数的，显然 \(g(n,j)=g(n,j-1)\) 。
否则 \(n\geq P_j^2\) ，显然会筛掉 \(P_j\) 的所有倍数（最小质因子为 $P_j $）。
不难发现 \(i^k\) 是一个完全积性函数，可以将需要删掉的 \(P_j\) 的倍数表示为 \(P_j\cdot s\) (\(2\leq s\leq \lfloor\frac{n}{P_j}\rfloor\)，且 \(s\) 的最小质因子大于等于 \(P_j\))，不难发现 \(g(\lfloor\frac{n}{P_j}\rfloor,j-1)\) 为所有合法的 \(s^k\) 的和加上 \(\sum_{i=1}^{j-1}P_i^k\) ，而 \(g(p_{j-1},j-1)\) 恰好等于 \(\sum_{i=1}^{j-1}P_i^k\) ，所以 \(g(\lfloor\frac{n}{P_j}\rfloor,j-1)-g(p_{j-1},j-1)\) 就是所有合法的 \(s^k\) 的和，而 \(i^k\) 是完全积性函数，只需要乘上 \(P_j^k\) 就是需要删掉的数。
边界条件 \(g(n,0)=\sum_{i=1}^n i^k\) 显然 \(k\) 较小的情况可以直接手推公式，\(k\) 较大可以考虑插值（\(k\) 大了时间复杂度爆炸）。
筛所有数的答案
如果分质数与合数讨论显然可以得到
前面一部分即为所有质数的贡献，而后面求合数的贡献，非常好理解就是暴力枚举 \(P_k\) 及其次幂然后求值，要注意的是 \(S(\lfloor\frac{n}{P_k^e}\rfloor,k+1)f(P_k^e)\) 并不包含出 \(P_k^{e+1}\) 的函数值，所以还需要加上 \(f(P_k^{e+1})\) 。
总复杂度 \(O(\frac{n^{\frac{3}{4}}}{\log n})\)
关于实现方面的细节
有一个性质 \(\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor=\lfloor\frac{n}{ab}\rfloor\) 所以我们只需要预处理所有的 \(g(\lfloor\frac{n}{x}\rfloor,j)\) 就可以了。
在这里如果使用 map 储存复杂度显然会变大，考虑到 \(\lfloor\frac{n}{\lfloor\frac{n}{x}\rfloor}\rfloor\) 与 \(\lfloor\frac{n}{x}\rfloor\) 中必然有一个数 \(\leq \sqrt{n}\) ，所以可以直接用两个数组储存编号。
\(S(n,j)\) 运算的时候不需要记忆化，因为调用的所有 \(S(n,j)\) 中没有相同的整数对 \((n,j)\) 。
\(\sum_{i=1}^ni=\frac{n(n+1)}{2}\)
\(\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}{6}\)
\(\sum_{i=1}^ni^3=\frac{n^2(n+1)^2}{4}\)
\(\sum_{i=1}^{n}i^4=\frac{n(n+1)(2n+1)(3n^2+3n-1)}{30}\)
举个简单的例子
#include<bits/stdc++.h>
using namespace std;
namespace Min_25
{
const int mod=1e9+7,inv2=(mod+1)/2,inv3=(mod+1)/3;
const int maxn=1e6+10;
int md(int x){return x>=mod?x-mod:x;}
bool flag[maxn];
int g1[maxn],g2[maxn];
long long w[maxn];
int id1[maxn],id2[maxn],m;
int sum1[maxn],sum2[maxn],lim;
int prime[maxn],cnt=0;
void make_prime(int n)
{
for(int i=2;i<=n;i++)
{
if(!flag[i])
{
prime[++cnt]=i;
sum1[cnt]=(sum1[cnt-1]+i);
sum2[cnt]=(sum2[cnt-1]+i*1LL*i)%mod;
}
for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
{
flag[prime[j]*i]=1;
if(i%prime[j]==0)break;
}
}
}
long long n;
void G()
{
m=0;
for(long long i=1,j;i<=n;i=j+1)
{
w[++m]=n/i,j=n/w[m];
if(w[m]<=lim)id1[w[m]]=m;
else id2[n/w[m]]=m;
int t=w[m]%mod;
g1[m]=(t+1)*1LL*t/2%mod-1;
g2[m]=(t+1)*1LL*t/2%mod*(t*2+1)%mod*inv3%mod-1;
}
for(int i=1;i<=cnt;i++)
for(int j=1;j<=m&&prime[i]*1LL*prime[i]<=w[j];j++)
{
int lst=w[j]/prime[i]<=lim?id1[w[j]/prime[i]]:id2[n/(w[j]/prime[i])];
g1[j]=md(g1[j]+mod-(g1[lst]-sum1[i-1]+mod)*1LL*prime[i]%mod);
g2[j]=md(g2[j]+mod-(g2[lst]-sum2[i-1]+mod)*1LL*prime[i]%mod*prime[i]%mod);
}
}
int S(long long a,int j)
{
if(prime[j]>=a)return 0;
int id=a<=lim?id1[a]:id2[n/a];
int ans=md(md(g2[id]+mod-g1[id])+mod-md(sum2[j]+mod-sum1[j]));
for(int k=j+1;prime[k]*1LL*prime[k]<=a&&k<=cnt;k++)
{
long long pe=prime[k];
for(int e=1;pe<=a;e++,pe*=prime[k])
{
long long tmp=pe%mod;
ans=(ans+tmp*(tmp-1)%mod*(S(a/pe,k)+(e!=1)))%mod;
}
}
return ans;
}
int Sum(long long N)
{
n=N;
lim=sqrt(n);
make_prime(lim);
G();
return (S(n,0)+1)%mod;
}
}
long long n;
int main()
{
scanf("%lld",&n);
printf("%d\n",Min_25::Sum(n));
}
杜教筛
常见积性函数
- \(\mu(x)\) 莫比乌斯反演函数。
- \(\varphi(x)\) 欧拉函数 \(\varphi(x)=\sum\limits_{i=1}^x[\gcd(x,i)=1]\) 。
- \(d(x)\) 约数个数函数 \(d(x)=\sum\limits_{i=1}^n[i|n]\)
- \(\sigma(x)\) 约数个数和函数 \(\sigma(x)=\sum\limits_{i=1}^n[i|n]i\)
- \(I(x)\) 恒等函数函数值恒为 \(1\)
- \(\epsilon(x)\) 元函数 \(\epsilon(x)=[x=1]\)
- \(id(x)\) 单位函数 \(id(x)=x\)
狄利克雷卷积
满足交换律，结合律，分配律。
不难发现存在 \(f*\epsilon=f\) 。
杜教筛
需要计算积性函数 \(f(x)\) 的前缀和。
记
考虑寻找到两个积性函数使得 \(h=g*f\) ，那么有
考虑将 \(d=1\) 对应的项提出来
也就是说只要 \(h(x)\) 前缀和很好求，那么这个式子就可以使用整除分块完成计算。
例子
-
\(\mu*I=\epsilon\)
-
\(\varphi*I=id\)
int smu[maxn];
long long sphi[maxn];
int mu[maxn],phi[maxn];
int prime[maxn],cnt=0,flag[maxn];
void make_prime(int n)
{
mu[1]=1,phi[1]=1;
for(int i=2;i<=n;i++)
{
if(!flag[i])
{
prime[++cnt]=i;
mu[i]=-1,phi[i]=i-1;
}
for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
{
flag[prime[j]*i]=1;
if(i%prime[j]==0)
{
mu[prime[j]*i]=0;
phi[prime[j]*i]=phi[i]*prime[j];
break;
}
mu[prime[j]*i]=-mu[i];
phi[prime[j]*i]=phi[i]*phi[prime[j]];
}
}
for(int i=1;i<=n;i++)smu[i]=smu[i-1]+mu[i],sphi[i]=sphi[i-1]+phi[i];
}
unordered_map<long long,long long>phimp;
unordered_map<long long,int>mump;
long long Sphi(long long n)
{
if(n<=5e6)return sphi[n];
if(phimp.count(n))return phimp[n];
long long ans=n*(n+1)/2;
for(long long i=2,j;i<=n;i=j+1)
{
j=n/(n/i);
ans-=(j-i+1)*Sphi(n/i);
}
return phimp[n]=ans;
}
int Smu(long long n)
{
if(n<=5e6)return smu[n];
if(mump.count(n))return mump[n];
int ans=1;
for(long long i=2,j;i<=n;i=j+1)
{
j=n/(n/i);
ans-=(j-i+1)*Smu(n/i);
}
return mump[n]=ans;
}
