数论算法
快速乘
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(leftlfloorfrac{a}{b} ight floor imes b ight) ight) y_{2})
(a y_{2}+b x_{2}-leftlfloorfrac{a}{b} ight floor imes b y_{2}=a y_{2}+bleft(x_{2}-leftlfloorfrac{a}{b} ight floor y_{2} ight))
所以 (x_{1}=y_{2}, y_{1}=x_{2}-leftlfloorfrac{a}{b} ight floor y_{2})
递归求解，终止时必然有 (x=1,y=0) 。
解集为 ({(x,y)|x=x_1+kfrac{b}{gcd(a,b)},y=y_1-kfrac{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(mod a_{1} ight)) 和 (x equiv b_2left(mod a_2 ight))
存在 (x=b_1+a_1k_1,x=b_2+a_2k_2 o b_1+a_1k_1=b_2+a_2k_2)
求解二元一次方程 (a_1k_1-a_2k_2=b_2-b_1)
设 (cequiv k_1(mod frac{a_2}{gcd(a_1,a_2)}))
可得新的方程为
(x equiv ca_{1}+b_{1}(mod 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(mod 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^xequiv b(mod{p}))，且 (a,p) 互质。
考虑分块求解这个问题，设 (n=lceilsqrt{q} ceil)
显然可以使用 (xn-y,x,yleq n) 来表示一个小于 (p) 的数
(a^{xn}equiv ba^{y}(mod{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^iequiv0(mod 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}{b_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 mod p_i^{k_i}) 再使用CRT合并。
只需要解决 (frac{n!}{p^a}mod p^k) ，其中 (a) 是 (n!) 包含的 (p) 的次数。
考虑将 (n!) 分组，(n!=(p imes2p imes 3p imes 4p imes... imeslfloorfrac{n}{p} floor p) imes1 imes2 imes... imes(p-1) imes(p+1) imes... imes(2p-1) imes...)
可以发现后面是模 (p^k) 是循环的，可以通过预处理求出
前面的 (lfloorfrac{n}{p} floor!) 是一个子问题递归求解即可。
复杂度 (O(plog 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^2equiv n(mod p)) ，由费马小定理可得 (x^{p-1}equiv 1(mod p))，显然 (n^{frac{p-1} {2}}equiv 1(mod p)) 。
若 (n) 为模 (p) 意义下的非二次剩余，假设存在 (x^2equiv n(mod p)) ，此时 (x^{p-1}equiv -1(mod p)) ，由费马小定理可知 (x) 不存在
若 $nequiv 0(mod p) $ ，显然成立。
Cipolla
定理1
定理2
(p) 的二次剩余与非二次剩余的个数均为 (frac{p-1}{2}) （不考虑 (0) 的情况下）。
证明：
我们只考虑所有的 (n^2)，假设有 (x≠y) 且 (x^2≡y^2(mod p))，则 (p∣(x^2−y^2))，即 (p∣(x−y)(x+y))，若(p∤(x−y))，则 (p∣(x+y))，故 (x+y≡0(mod p))，也就是不同的 (x^2) 共有 $frac{p-1}{2} $ 个，即二次剩余有个 (frac{p-1}{2}) 个。
算法流程
- 判断给定的数 (x) 是否是二次剩余。
- 随机一个 (a)，使其满足 ((a^2−x)) 是二次非剩余。
- 令 (omega^2equiv (a^2-x)pmod{p}) ，取 (yequiv left(a+{omega} ight)^{p+1over 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) 是否为素数（(qleq 10^{18})）
根据费马小定理，当 (a mod q ≠0) 且 (q) 为质数时，(a^{q-1}≡1 (mod q))
而当 (q) 不是质数时，不一定成立
因此就有了一个想法：随机一些 (a) ，对每一个 (a) 验证 ，如果有一个错误那么它一定是合数
但是，存在这样一类合数 (q) ，对于 (1≤a<q) 上面的式子都成立，比如 (561)
二次探测定理
如果 $x^2≡1(mod p) $，那么 (x≡1(mod p)) 或者 (x≡p-1(mod p)) ，这里 (p) 是质数。
可以发现对于合数，有一定概率不成立。
在检验时，如果 (a^{q-1}≡1(mod q)) 且 (2|(q-1)) ，则下一步计算 (a ^{frac{q-1}{2}} mod 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(klog n))
莫比乌斯反演
莫比乌斯反演函数
性质
莫比乌斯反演定理
若函数 (F(n)) 满足
则存在
证明：
考虑 (f(i)) 的系数可以很简单地得到证明。
推论：
若函数 (F(d)) 满足
则存在
证明：
考虑每个 (f(kcdot d)) 的系数被贡献时满足 (icdot j=k)，显然 (i) 是 (k) 的约数，被贡献的系数值为 (mu(i))。
这个推论可以很简单地用于求解类似 (sum_{i=1}^nsum_{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_1leq k_1,t_2leq 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=lfloorfrac{N}{i} floor)求出 (sum_{i=1}^n[iin P]f(i))
首先线性筛出 (sqrt{n}) 以内的所有质数，(P_i) 表示第 (i) 个质数。
因为 (f(p)) 能够写成一个简单多项式
所以本质上来说需要求解的是 (sum_{i=1}^n [iin 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)) 。
否则 (ngeq P_j^2) ，显然会筛掉 (P_j) 的所有倍数（最小质因子为 $P_j $）。
不难发现 (i^k) 是一个完全积性函数，可以将需要删掉的 (P_j) 的倍数表示为 (P_jcdot s) ((2leq sleq lfloorfrac{n}{P_j} floor)，且 (s) 的最小质因子大于等于 (P_j))，不难发现 (g(lfloorfrac{n}{P_j} floor,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(lfloorfrac{n}{P_j} floor,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(lfloorfrac{n}{P_k^e} floor,k+1)f(P_k^e)) 并不包含出 (P_k^{e+1}) 的函数值，所以还需要加上 (f(P_k^{e+1})) 。
总复杂度 (O(frac{n^{frac{3}{4}}}{log n}))
关于实现方面的细节
有一个性质 (lfloorfrac{lfloorfrac{n}{a} floor}{b} floor=lfloorfrac{n}{ab} floor) 所以我们只需要预处理所有的 (g(lfloorfrac{n}{x} floor,j)) 就可以了。
在这里如果使用 map 储存复杂度显然会变大，考虑到 (lfloorfrac{n}{lfloorfrac{n}{x}
floor}
floor) 与 (lfloorfrac{n}{x}
floor) 中必然有一个数 (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
",Min_25::Sum(n));
}
杜教筛
常见积性函数
- (mu(x)) 莫比乌斯反演函数。
- (varphi(x)) 欧拉函数 (varphi(x)=sumlimits_{i=1}^x[gcd(x,i)=1]) 。
- (d(x)) 约数个数函数 (d(x)=sumlimits_{i=1}^n[i|n])
- (sigma(x)) 约数个数和函数 (sigma(x)=sumlimits_{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;
}