数论与数论函数

数论与数论函数

全文参考oi-wiki.org

基础知识

快速幂

ll qpow(ll a,ll b,ll mod){
	ll res=1;
	a=a%mod;
	while(b){
		if(b&1)res=res*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return res;
}

可以分块预处理任意进制的快速幂,或者“光速幂”。

快速乘法

用于 a ∗ b a*b ab l o n g l o n g long long longlong的情况

ll qmul(ll a,ll b,ll mod){
	ll res=0;
	while(b){
		if(b&1)res=(res+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return res;
}

还有更快的方法,传说中的 O ( 1 ) O(1) O(1)快速乘:

ll qmul(ll x,ll y,ll mod){
    return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;     
}

整除及其性质

素数

1、素数的个数:
素数计数函数:小于等于 x x x的素数的个数,用 π ( x ) \pi(x) π(x)表示,随着 x x x的增大,有近似结论: π ( x ) ∼ x ln ⁡ ( x ) \pi(x)\sim\frac{x}{\ln(x)} π(x)ln(x)x
2、素数判定:
暴力判定: O ( n ) O(\sqrt {n}) O(n )
素性测试: O ( k log ⁡ 3 n ) O(k \log^3n) O(klog3n)(详见匡斌模板)

最大公约数

欧几里得算法
ll gcd(ll a,ll b){
	return b==0?a:gcd(b,a%b);
}
多个数的最大公约数

g c d ( a , b , c ) = g c d ( g c d ( a , b ) , c ) gcd(a,b,c)=gcd(gcd(a,b),c) gcd(a,b,c)=gcd(gcd(a,b),c)

最小公倍数

l c m ( a , b ) = a ∗ b / g c d ( a , b ) lcm(a,b)=a*b/gcd(a,b) lcm(a,b)=ab/gcd(a,b)

多个数的最小公倍数

l c m ( a , b , c ) = l c m ( l c m ( a , b ) , c ) lcm(a,b,c)=lcm(lcm(a,b),c) lcm(a,b,c)=lcm(lcm(a,b),c)
谨记 l c m ( a , b , c ) ≠ a ∗ b ∗ c / g c d ( a , b , c ) lcm(a,b,c)\neq a*b*c/gcd(a,b,c) lcm(a,b,c)=abc/gcd(a,b,c)

扩展欧几里得定理

常用于求解 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b)

int Exgcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  }
  int d = Exgcd(b, a % b, x, y);
  int t = x;
  x = y;
  y = t - (a / b) * y;
  return d;
}

函数返回值为 g c d ( a , b ) gcd(a,b) gcd(a,b),如果 a x + b y = c ax+by=c ax+by=c c c c不是 g c d ( a , b ) gcd(a,b) gcd(a,b)的倍数,则方程无解。

欧拉函数

φ ( n ) \varphi(n) φ(n)表示小于等于 n n n n n n互质的数的个数。
特别地,令 φ ( 1 ) = 1 \varphi(1)=1 φ(1)=1
n = p 1 k 1 p 2 k 2 … p s k s n=p_1^{k_1}p_2^{k_2}…p_s^{k_s} n=p1k1p2k2psks,那么 φ ( n ) = n ∗ ∏ i = 1 s p i − 1 p i \varphi(n)=n*\prod_{i=1}^s{\frac{p_i-1}{p_i}} φ(n)=ni=1spipi1
n = p k n=p^k n=pk, p p p是质数,那么 φ ( n ) = p k − 1 ∗ ( p − 1 ) \varphi(n)=p^{k-1}*(p-1) φ(n)=pk1(p1)

欧拉定理

g c d ( a , m ) = 1 gcd(a,m)=1 gcd(a,m)=1,则 a φ ( m ) ≡ 1 ( m o d m ) a^{\varphi(m)}\equiv 1\pmod{m} aφ(m)1(modm)
当m为素数时即为费马小定理。

扩展欧拉定理

a b ≡ { a b   m o d   φ ( p ) ,   gcd ⁡ ( a ,   p ) = 1 a b , gcd ⁡ ( a , p ) ≠ 1 , b < φ ( p ) a b   m o d   φ ( p ) + φ ( p ) , gcd ⁡ ( a , p ) ≠ 1 , b ≥ φ ( p ) ( m o d p ) a^b\equiv\begin{cases} a^{b\bmod{\varphi(p)}},\,&\gcd(a,\,p)=1\\ a^b,& \gcd(a,p)\ne1,b<\varphi(p)\\ a^{b\bmod{\varphi(p)+\varphi(p)}},&\gcd(a,p)\ne1,b\ge\varphi(p)\end{cases}\pmod{p} ababmodφ(p),ab,abmodφ(p)+φ(p),gcd(a,p)=1gcd(a,p)=1,b<φ(p)gcd(a,p)=1,bφ(p)(modp)
这个定理常用来欧拉降幂

同余与同余方程

剩余系

完全剩余系:对于 n n n,就是集合 { 0 , 1 , … , n − 1 } \{0,1,…,n-1\} {0,1,,n1}
简化剩余系:又叫既约剩余系,对于 n n n,是完全剩余系中与 n n n互质的元素构成的子集。

裴蜀定理

a , b a,b a,b是不全为零的整数,则存在整数 x , y x,y x,y,使得 a x + b y = gcd ⁡ ( a , b ) ax+by=\gcd(a,b) ax+by=gcd(a,b)
这个定理可以扩展到 n n n个数的情况

乘法逆元

如果一个线性同余方程 a x ≡ 1 ( m o d b ) ax\equiv 1\pmod b ax1(modb),则 x x x称为 a   m o d   b a \bmod b amodb 的逆元,记作 a − 1 a^{-1} a1.

扩展欧几里得法
void exgcd(int a, int b, int& x, int& y) {
  if (b == 0) {
    x = 1, y = 0;
    return;
  }
  exgcd(b, a % b, y, x);
  y -= a / b * x;
}
费马小定理

a − 1 = a p − 2   m o d   p a^{-1}=a^{p-2}\bmod p a1=ap2modp,可以用快速幂求 a p − 2 a^{p-2} ap2

线性求逆元

1 , 2 , … , n 1,2,…,n 1,2,,n中每个数关于 p p p的逆元
首先,有 1 − 1 ≡ 1 ( m o d p ) 1^{-1}\equiv 1\pmod p 111(modp)
p = k i + j , j < i , 1 < i < p p=ki+j,j<i,1<i<p p=ki+j,j<i,1<i<p
k i + j ≡ 0 ( m o d p ) ki+j\equiv 0\pmod p ki+j0(modp)
两边同时乘 i − 1 , j − 1 i^{-1},j^{-1} i1,j1:
k j − 1 + i − 1 ≡ 0 ( m o d p ) kj^{-1}+i^{-1}\equiv 0\pmod p kj1+i10(modp)
i − 1 ≡ − k j − 1 ( m o d p ) i^{-1}\equiv -kj^{-1}\pmod p i1kj1(modp)
i − 1 ≡ − ( p i ) ( p   m o d   i ) − 1 ( m o d p ) i^{-1}\equiv -(\frac{p}{i})(p\bmod i)^{-1}\pmod p i1(ip)(pmodi)1(modp)

inv[1] = 1;
for (int i = 2; i <= n; ++i) 
	inv[i] = (long long)(p - p / i) * inv[p % i] % p;
线性求任意 n n n个数的逆元

如果要求任意给定 n n n个数的逆元( 1 ≤ a i < p 1\le a_i<p 1ai<p)。
首先计算 n n n个数的前缀积,记为 s n s_n sn,求得其逆元,记为 s v n sv_n svn,当把它乘上 a n a_n an时,就会得到 a 1 a_1 a1 a n − 1 a_{n-1} an1的积的逆元,记为 s v n − 1 sv_{n-1} svn1,同理可得所有的 s v i sv_i svi,于是 a i − 1 a_i^{-1} ai1就可以用 s i − 1 ∗ s v i s_{i-1}*sv_i si1svi得到。

线性同余方程

形如 a x ≡ b ( m o d c ) ax\equiv b\pmod c axb(modc)的方程被称为线性同余方程,
其解和 a x + c y = b ax+cy=b ax+cy=b是等价的
我们可以先解出 a x + b y = gcd ⁡ ( a , b ) ax+by=\gcd(a,b) ax+by=gcd(a,b)的一组 x 0 , y 0 x_0,y_0 x0,y0,然后两边同时除以 gcd ⁡ ( a , b ) \gcd(a,b) gcd(a,b),在乘 c c c。就找到了方程的一个解。
gcd ⁡ ( a , b ) = 1 \gcd(a,b)=1 gcd(a,b)=1,且 x 0 , y 0 x_0,y_0 x0,y0为方程 a x + b y = c ax+by=c ax+by=c的一组解,则该方程的任意解可表示为: x = x 0 + b t , y = y 0 − a t x=x_0+bt,y=y_0-at x=x0+bt,y=y0at,对任意整数t都成立。

int ex_gcd(int a, int b, int& x, int& y) {
  if (b == 0) {
    x = 1;
    y = 0;
    return a;
  }
  int d = ex_gcd(b, a % b, x, y);
  int temp = x;
  x = y;
  y = temp - a / b * y;
  return d;
}
bool liEu(int a, int b, int c, int& x, int& y) {
  int d = ex_gcd(a, b, x, y);
  if (c % d != 0) return 0;
  int k = c / d;
  x *= k;
  y *= k;
  return 1;
}

中国剩余定理

一般这玩意都是用扩展中国剩余定理(模数不互质)

#include<iostream>
#include<vector>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long lt;

lt read()
{
    lt f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}

const int maxn=100010;
int n;
lt ai[maxn],bi[maxn];

lt mul(lt a,lt b,lt mod)
{
    lt res=0;
    while(b>0)
    {
        if(b&1) res=(res+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return res;
}

lt exgcd(lt a,lt b,lt &x,lt &y)
{
    if(b==0){x=1;y=0;return a;}
    lt gcd=exgcd(b,a%b,x,y);
    lt tp=x;
    x=y; y=tp-a/b*y;
    return gcd;
}

lt excrt()
{
    lt x,y,k;
    lt M=bi[1],ans=ai[1];//第一个方程的解特判
    for(int i=2;i<=n;i++)
    {
        lt a=M,b=bi[i],c=(ai[i]-ans%b+b)%b;//ax≡c(mod b)
        lt gcd=exgcd(a,b,x,y),bg=b/gcd;
        if(c%gcd!=0) return -1; //判断是否无解,然而这题其实不用

        x=mul(x,c/gcd,bg);
        ans+=x*M;//更新前k个方程组的答案
        M*=bg;//M为前k个m的lcm
        ans=(ans%M+M)%M;
    }
    return (ans%M+M)%M;
}

int main()
{
    n=read();
    for(int i=1;i<=n;++i)
    bi[i]=read(),ai[i]=read();//读入时看好先读哪个
    printf("%lld",excrt());
    return 0;
}

BSGS

常用来在 O ( p ) O(\sqrt{p}) O(p )的时间复杂度内求解 a x ≡ b   m o d   p a^x\equiv b\bmod p axbmodp,其中 a a a p p p互质,方程的解 x x x满足 0 ≤ x < p 0\le x< p 0x<p
x = A ⌈ p ⌉ − B x=A\lceil\sqrt{p}\rceil-B x=Ap B,其中 0 ≤ A , B ≤ ⌈ p ⌉ 0\le A,B\le \lceil\sqrt{p}\rceil 0A,Bp ,则有 a A ⌈ p ⌉ − B ≡ b ( m o d p ) a^{A\lceil\sqrt{p}\rceil-B}\equiv b\pmod p aAp Bb(modp),则有 a A ⌈ p ⌉ ≡ b a B a^{A\lceil\sqrt{p}\rceil}\equiv ba^B aAp baB.
首先,计算并记录所有 b a B ba^B baB的取值,再枚举 A A A,计算 a A ⌈ p ⌉ a^{A\lceil\sqrt{p}\rceil} aAp 是否有 b a B ba^B baB与之相等,从而可以得到所有的 x x x x = A ⌈ p ⌉ − B x=A\lceil\sqrt{p}\rceil-B x=Ap B
待补充……

原根

( a , m ) = 1 (a,m)=1 (a,m)=1,使得 a l ≡ 1 ( m o d m ) a^l\equiv 1\pmod m al1(modm)成立最小的 l l l,称为 a a a关于模 m m m的阶,记为 o r d m a ord_ma ordma
o r d m a = l ord_ma=l ordma=l,则 o r d m a t = l ( t , l ) ord_ma^t=\frac{l}{(t,l)} ordmat=(t,l)l
由欧拉定理,设 o r d m a = l ord_ma=l ordma=l,则 a n ≡ 1 ( m o d m ) a^n\equiv 1\pmod m an1(modm)当且仅当 l ∣ n l|n ln,特别地, l ∣ φ ( m ) l|\varphi(m) lφ(m).
性质:
1、设 p p p是素数, o r d p a = l ord_pa=l ordpa=l,那么有且仅有 φ ( l ) \varphi(l) φ(l)个关于模 p p p的阶为 l l l且两两互不同余的数
2、设 o r d m a = l ord_ma=l ordma=l,则 1 , a , a 2 , … , a l − 1 1,a,a^2,…,a^{l-1} 1,a,a2,,al1关于模 m m m两两互不同余
3、设 p p p是素数, l ∣ p − 1 l|p-1 lp1,则存在 φ ( l ) \varphi(l) φ(l)个关于模 p p p的阶为 l l l且两两互不同余的数。
4、若 m = p 1 a 1 p 2 a 2 … p k a k m=p_1^{a_1}p_2^{a_2}…p_k^{a_k} m=p1a1p2a2pkak,则 o r d m a = [ o r d p 1 a 1 , o r d p 2 a 2 , … , o r d p k a k ] ord_ma=[ord_{p_1}^{a_1},ord_{p_2}^{a_2},…,ord_{p_k}^{a_k}] ordma=[ordp1a1,ordp2a2,,ordpkak]

原根

( g , m ) = 1 , (g,m)=1, (g,m)=1, o r d m g = φ ( m ) ord_mg=\varphi(m) ordmg=φ(m),则称 g g g m m m的一个原根。
g g g m m m的一个原根当且仅当{ g , g 2 , … , g φ ( m ) g,g^2,…,g^{\varphi(m)} g,g2,,gφ(m)}构成模 m m m的一个既约剩余系。

是否有原根

m m m有原根,则 m m m一定是: 2 , 4 , p a , 2 p a 2,4,p^a,2p^a 2,4,pa,2pa,其中 p p p是奇素数, a a a为正整数。

求所有原根

g g g m m m的一个原根,则集合 S = { g s ∣ 1 ≤ s ≤ φ ( m ) , ( s , φ ( m ) ) = 1 } S=\{g^s|1\le s\le \varphi(m),(s,\varphi(m))=1\} S={gs1sφ(m),(s,φ(m))=1}给出 m m m的全部原根。

求一个原根

p − 1 p-1 p1进行质因数分解得到不同的质因子 d 1 , d 2 , … , d m d_1,d_2,…,d_m d1,d2,,dm,对于任意的 1 < a < p 1<a<p 1<a<p,要判定 a a a是否是模 p p p的原根,只需要检验 a p − 1 d 1 , a p − 1 d 2 , … , a p − 1 d m a^{\frac{p-1}{d_1}},a^{\frac{p-1}{d_2}},…,a^{\frac{p-1}{d_m}} ad1p1,ad2p1,,admp1 m m m个数中是否存在一个数在模 p p p意义下与1同余。若存在,则 a a a不是 p p p的原根;若不存在,则 a a a p p p的原根。
或者说: ( g , m ) = 1 (g,m)=1 (g,m)=1,设 p 1 , p 2 , … , p k p_1,p_2,…,p_k p1,p2,,pk φ ( m ) \varphi(m) φ(m)的所有不同的素因子,则 g g g m m m的原根,当且仅当对任意 1 ≤ i ≤ k 1\le i\le k 1ik,都有 g φ ( m ) p i ≢ 1 ( m o d m ) g^{\frac{\varphi(m)}{p_i}}\not\equiv 1\pmod m gpiφ(m)1(modm)

注:998244353 的原根是3;1e9+7的原根是5

下面给出暴力求原根的代码:(python版本)

def check(p):
    fac = []
    n = p - 1
    i = 2
    while i * i <= n:
        if n % i == 0:
            fac.append(i)
            while n % i == 0:
                n = n // i
        i += 1
    fac.append(n)
    i = 2
    while i < p:
        flag = True
        for j in fac:
            if pow(i, (p - 1) // j, p) == 1:
                flag = False
                break
        if flag:
            return i
        i += 1


n = int(input())
print(check(n))

卢卡斯定理(大组合数取模)

Lucas定理内容如下:
C n m   m o d   p = C ⌊ n p ⌋ ⌊ m p ⌋ ⋅ C n   m o d   p m   m o d   p   m o d   p C_n^m\bmod p=C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{m}{p}\rfloor}\cdot C_{n\bmod p}^{m\bmod p}\bmod p Cnmmodp=CpnpmCnmodpmmodpmodp

long long Lucas(long long n, long long m, long long p) {
  if (m == 0) return 1;
  return (C(n % p, m % p, p) * Lucas(n / p, m / p, p)) % p;
}
exLucas定理

当模数 p p p不是素数时,就需要用到exLucas定理,待补充。。
下面是代码,求 C n m   m o d   p C_n^m\bmod p Cnmmodp

#include<bits/stdc++.h>
#define ll long long
using namespace std;
#ifndef Fading
inline char gc(){
    static char now[1<<16],*S,*T;
    if (T==S){T=(S=now)+fread(now,1,1<<16,stdin);if (T==S) return EOF;}
    return *S++;
}
#endif
#ifdef Fading
#define gc getchar
#endif
void exgcd(ll a,ll b,ll &x,ll &y){
    if (!b) return (void)(x=1,y=0);
    exgcd(b,a%b,x,y);
    ll tmp=x;x=y;y=tmp-a/b*y;
}
ll gcd(ll a,ll b){
    if (b==0) return a;
    return gcd(b,a%b); 
}
inline ll INV(ll a,ll p){
    ll x,y;
    exgcd(a,p,x,y);
    return (x+p)%p;
}
inline ll lcm(ll a,ll b){
    return a/gcd(a,b)*b;
}
inline ll mabs(ll x){
    return (x>0?x:-x);
}
inline ll fast_mul(ll a,ll b,ll p){
    ll t=0;a%=p;b%=p;
    while (b){
        if (b&1LL) t=(t+a)%p;
        b>>=1LL;a=(a+a)%p;
    }
    return t;
}
inline ll fast_pow(ll a,ll b,ll p){
    ll t=1;a%=p;
    while (b){
        if (b&1LL) t=(t*a)%p;
        b>>=1LL;a=(a*a)%p;
    }
    return t;
}
inline ll read(){
    ll x=0,f=1;char ch=gc();
    while (!isdigit(ch)) {if (ch=='-') f=-1;ch=gc();}
    while (isdigit(ch)) x=x*10+ch-'0',ch=gc();
    return x*f;
}
inline ll F(ll n,ll P,ll PK){
    if (n==0) return 1;
    ll rou=1;//循环节
    ll rem=1;//余项 
    for (ll i=1;i<=PK;i++){
        if (i%P) rou=rou*i%PK;
    }
    rou=fast_pow(rou,n/PK,PK);
    for (ll i=PK*(n/PK);i<=n;i++){
        if (i%P) rem=rem*(i%PK)%PK;
    }
    return F(n/P,P,PK)*rou%PK*rem%PK;
}
inline ll G(ll n,ll P){
    if (n<P) return 0;
    return G(n/P,P)+(n/P);
}
inline ll C_PK(ll n,ll m,ll P,ll PK){
    ll fz=F(n,P,PK),fm1=INV(F(m,P,PK),PK),fm2=INV(F(n-m,P,PK),PK);
    ll mi=fast_pow(P,G(n,P)-G(m,P)-G(n-m,P),PK);
    return fz*fm1%PK*fm2%PK*mi%PK;
}
ll A[1001],B[1001];
//x=B(mod A)
inline ll exLucas(ll n,ll m,ll P){
    ll ljc=P,tot=0;
    for (ll tmp=2;tmp*tmp<=P;tmp++){
        if (!(ljc%tmp)){
            ll PK=1;
            while (!(ljc%tmp)){
                PK*=tmp;ljc/=tmp;
            }
            A[++tot]=PK;B[tot]=C_PK(n,m,tmp,PK);
        }
    }
    if (ljc!=1){
        A[++tot]=ljc;B[tot]=C_PK(n,m,ljc,ljc);
    }
    ll ans=0;
    for (ll i=1;i<=tot;i++){
        ll M=P/A[i],T=INV(M,A[i]);
        ans=(ans+B[i]*M%P*T%P)%P;
    }
    return ans;
}
signed main(){
    ll n=read(),m=read(),P=read();
    printf("%lld\n",exLucas(n,m,P));
    return 0;
}

数论函数

线性筛法

线性筛法求素数:

void Prime(){
	memset(vis,false,sizeof(vis));
	memset(prime,0,sizeof(prime));
	for(int i=2;i<=maxn;i++){
		if(!vis[i]){
			prime[++cnt]=i;
		}
		for(int j=1;j<=cnt&&i*prime[j]<=maxn;j++){
			vis[i*prime[j]]=true;
			if(i%prime[j]==0){
				break;
			}
		}
	}
}

线性筛法求欧拉函数:

void euler(){
    m = 0;
    memset(v, 0, sizeof(v));
    for(int i = 2; i < maxn; i++){
        if(v[i] == 0){
            v[i] = i;
            p[m++] = i;
            phi[i] = i - 1;
        }
        for(int j = 0; j < m; j++){
            if(p[j] > v[i] || p[j] > maxn / i)
                break;
            v[i * p[j]] = p[j];
            phi[i * p[j]] = phi[i] * (i % p[j] ? p[j] - 1 : p[j]);
        }
    }
}

线性筛法求莫比乌斯函数:

void pre() {
  mu[1] = 1;
  for (int i = 2; i <= 1e7; ++i) {
    if (!v[i]) mu[i] = -1, p[++tot] = i;
    for (int j = 1; j <= tot && i <= 1e7 / p[j]; ++j) {
      v[i * p[j]] = 1;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }

线性筛法求约数个数:

void pre() {
  d[1] = 1;
  for (int i = 1; i <= n; ++i) {
    if (!v[i]) v[i] = 1, p[++tot] = i, d[i] = 2, num[i] = 1;
    for (int j = 1; j <= tot && i <= n / p[j]; ++j) {
      v[p[j] * i] = 1;
      if (i % p[j] == 0) {
        num[i * p[j]] = num[i] + 1;
        d[i * p[j]] = d[i] / num[i * p[j]] * (num[i * p[j]] + 1);
        break;
      } else {
        num[i * p[j]] = 1;
        d[i * p[j]] = d[i] * 2;
      }
    }
  }
}

线性筛法求约数和:

void pre() {
  g[1] = f[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!v[i]) v[i] = 1, p[++tot] = i, g[i] = i + 1, f[i] = i + 1;
    for (int j = 1; j <= tot && i <= n / p[j]; ++j) {
      v[p[j] * i] = 1;
      if (i % p[j] == 0) {
        g[i * p[j]] = g[i] * p[j] + 1;
        f[i * p[j]] = f[i] / g[i] * g[i * p[j]];
        break;
      } else {
        f[i * p[j]] = f[i] * f[p[j]];
        g[i * p[j]] = 1 + p[j];
      }
    }
  }
  for (int i = 1; i <= n; ++i) f[i] = (f[i - 1] + f[i]) % Mod;
}
非常见积性函数的线性筛法

这里主要讨论求两个积性函数狄利克雷卷积的情况,这个函数显然也是积性函数。
l o w ( i ) low(i) low(i)表示 p 1 a 1 p_1^{a_1} p1a1,对于线性筛法最为关键的 i % p j = 0 i\%p_j=0 i%pj=0;
有了 l o w ( i ) low(i) low(i),我们可以分两种情况来讨论。
1、若 l o w ( i ) = i low(i)=i low(i)=i,此时 i i i一定是某个素数的幂次的形式,(否则会被 b r e a k break break),
此时只要我们可以快速地从 f ( p k ) f(p^k) f(pk)来更新 f ( p k + 1 ) f(p^{k+1}) f(pk+1)
2、若 l o w ( i ) ! = i {low(i)}!=i low(i)!=i,此时 i / l o w ( i ) i/low(i) i/low(i)一定与 l o w ( i ) ∗ p j low(i)*p_j low(i)pj是互质的,
我们可以直接利用积性函数的性质去更新。
下面是伪代码:

vis[1] = low[1] = 1; H[1] = 初始化 
for(int i = 2; i <= N; i++) {
    if(!vis[i]) prime[++tot] = i, mu[i] = -1, H[i] = 质数的情况, low[i] = i;
    for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
        vis[i * prime[j]] = 1;
        if(!(i % prime[j])) {
            low[i * prime[j]] = (low[i] * prime[j]); 
            if(low[i] == i) H[i * prime[j]] = 特殊判断;
            else H[i * prime[j]] = H[i / low[i]] * H[prime[j] * low[i]];
            break;
        } 
        H[i * prime[j]] = H[i] * H[prime[j]];
        low[i * prime[j]] = prime[j];
    }
}

莫比乌斯反演

数论分块
引理1

∀ a , b , c ∈ Z , ⌊ a b c ⌋ = ⌊ ⌊ a b ⌋ c ⌋ \forall a,b,c\in Z,\lfloor \frac{a}{bc} \rfloor=\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}\rfloor a,b,cZ,bca=cba

引理2

∀ n ∈ N , ∣ { ⌊ n d ∣ d ∈ N ⌋ } ∣ ≤ ⌊ 2 n ⌋ \forall n\in N,\vert\{\lfloor\frac{n}{d}|d\in N \rfloor\}\vert \leq \lfloor 2\sqrt{n} \rfloor nN,{dndN}2n
∣ V ∣ 表 示 集 合 V 的 元 素 个 数 \vert V \vert表示集合V的元素个数 VV

数论分块

对于含有 ⌊ n i ⌋ \lfloor \frac{n}{i} \rfloor in的求和式子(n为常数)
对于任意一个 i ( i ≤ n ) i(i \leq n) i(in),我们需要找到一个最大的 j ( i ≤ j ≤ n ) j(i \leq j \leq n) j(ijn),使得 ⌊ n i ⌋ = ⌊ n j ⌋ . \lfloor\frac{n}{i}\rfloor=\lfloor \frac{n}{j} \rfloor. in=jn.则有 j = ⌊ n ⌊ n i ⌋ ⌋ j=\lfloor \frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor j=inn.
所以我们可以每次以 [ i , j ] [ i,j] [i,j]为一块,分块求和即可。

积性函数
定义:

若 g c d ( x , y ) = 1 且 f ( x y ) = f ( x ) f ( y ) , 则 f ( n ) 为 积 性 函 数 若gcd(x,y)=1且f(xy)=f(x)f(y),则f(n)为积性函数 gcd(x,y)=1f(xy)=f(x)f(y),f(n)

性质:

若 f ( x ) 和 g ( x ) 均 为 积 性 函 数 , 则 以 下 函 数 也 为 积 性 函 数 : 若f(x)和g(x)均为积性函数,则以下函数也为积性函数: f(x)g(x)
h ( x ) = f ( x p ) \quad\quad\quad\quad h(x)=f(x^p) h(x)=f(xp)
h ( x ) = f p ( x ) \quad\quad\quad\quad h(x)=f^p(x) h(x)=fp(x)
h ( x ) = f ( x ) g ( x ) \quad\quad\quad\quad h(x)=f(x)g(x) h(x)=f(x)g(x)
h ( x ) = ∑ d ∣ x f ( d ) g ( x d ) \quad\quad\quad\quad h(x)=\sum\limits_{d|x}f(d)g(\frac{x}{d}) h(x)=dxf(d)g(dx)

常见积性函数:

1. 幺 元 函 数 : ϵ ( n ) = [ n = 1 ] 1.幺元函数: \epsilon(n)=[n=1] 1.ϵ(n)=[n=1]
2. 常 函 数 : o n e ( x ) = 1 2.常函数:one(x)=1 2.one(x)=1
3. 标 号 函 数 : i d ( x ) = x 3.标号函数:id(x)=x 3.id(x)=x
4. 欧 拉 函 数 : φ ( n ) = ∑ i = 1 n [ g c d ( i , n ) = 1 ] 4.欧拉函数:\varphi(n)=\sum_{i=1}^n[gcd(i,n)=1] 4.φ(n)=i=1n[gcd(i,n)=1]
5. 除 数 函 数 : 5.除数函数: 5.
σ ( k , x ) = ∑ d ∣ x d k \quad\quad\quad \sigma(k,x)=\sum_{d|x}d^k σ(k,x)=dxdk
当 k = 1 时 , 表 示 x 的 因 子 之 和 , 简 记 为 σ ( n ) \quad当k=1时,表示x的因子之和,简记为\sigma(n) k=1x,σ(n)
当 k = 0 时 , 表 示 x 的 因 子 个 数 , 简 记 为 d ( n ) \quad当k=0时,表示x的因子个数,简记为d(n) k=0x,d(n)
\quad k省略时默认为1
6. 莫 比 乌 斯 函 数 : μ ( n ) = { 1 n = 1 0 ∃ d : d 2 ∣ n ( − 1 ) ω ( n ) o t h e r w i s e 6.莫比乌斯函数:\mu(n)=\begin{cases}1&n=1\\0&\exists d:d^2|n\\(-1)^{\omega(n)}&otherwise\end{cases} 6.μ(n)=10(1)ω(n)n=1d:d2notherwise其中 ω ( n ) 表 示 n 的 本 质 不 同 的 质 因 子 个 数 , 是 一 个 加 性 函 数 \omega(n)表示n的本质不同的质因子个数,是一个加性函数 ω(n)n

Dirichlet 卷积
定义

两 个 数 论 函 数 f , g 的 D i r i c h l e t 卷 积 为 ( f ∗ g ) ( n ) = ∑ d ∣ n f ( d ) g ( n d ) 两个数论函数f,g的Dirichlet卷积为(f*g)(n)=\sum\limits_{d|n}f(d)g(\frac{n}{d}) f,gDirichlet(fg)(n)=dnf(d)g(dn)

性质

D i r i c h l e t 卷 积 满 足 交 换 律 、 结 合 律 和 分 配 率 . Dirichlet卷积满足交换律、结合律和分配率. Dirichlet.
交 换 律 : f ∗ g = g ∗ f ; 交换律:f*g=g*f; fg=gf;
结 合 律 : ( f ∗ g ) ∗ h = f ∗ ( g ∗ h ) ; 结合律:(f*g)*h=f*(g*h); (fg)h=f(gh);
分 配 率 : ( f + g ) ∗ h = f ∗ h + g ∗ h ; 分配率:(f+g)*h=f*h+g*h; (f+g)h=fh+gh;
ε 为 D i r i c h l e t 卷 积 中 的 单 位 元 ( 任 何 函 数 卷 ε 都 是 其 本 身 ) \varepsilon 为Dirichlet卷积中的单位元(任何函数卷\varepsilon都是其本身) εDirichletε

例子

ε = μ ∗ 1 ⇔ ε ( n ) = ∑ d ∣ n μ ( d ) \varepsilon =\mu*1\Leftrightarrow\varepsilon(n)=\sum\limits_{d|n}\mu(d) ε=μ1ε(n)=dnμ(d)
d = 1 ∗ 1 ⇔ d ( n ) = ∑ d ∣ n 1 d =1*1 \Leftrightarrow d(n)=\sum\limits_{d|n}1 d=11d(n)=dn1
σ = d ∗ 1 ⇔ σ ( n ) = ∑ d ∣ n d \sigma=d*1\Leftrightarrow\sigma(n)=\sum\limits_{d|n}d σ=d1σ(n)=dnd
φ = μ ∗ i d ⇔ φ ( n ) = ∑ d ∣ n d ⋅ μ ( n d ) \varphi =\mu*id\Leftrightarrow\varphi(n)=\sum\limits_{d|n}d\cdot \mu(\frac{n}{d}) φ=μidφ(n)=dndμ(dn)
i d = φ ∗ 1 ⇔ n = ∑ d ∣ n φ ( d ) id=\varphi*1\Leftrightarrow n=\sum\limits_{d|n}\varphi(d) id=φ1n=dnφ(d)
g = f ∗ μ ⇔ f = g ∗ 1 g=f*\mu\Leftrightarrow f=g*1 g=fμf=g1

莫比乌斯反演

μ ( n ) = { 1 n = 1 0 n 含 有 平 方 因 子 ( − 1 ) k k 为 n 的 本 质 不 同 质 因 子 个 数 \mu(n)=\begin{cases}1&n=1\\0&n含有平方因子\\(-1)^k&k为n的本质不同质因子个数\end{cases} μ(n)=10(1)kn=1nkn
反演重要结论:
ε ( n ) = ∑ d ∣ n μ ( d ) = [ n = 1 ] \varepsilon(n)=\sum\limits_{d|n}\mu(d)=[n=1] ε(n)=dnμ(d)=[n=1]
[ g c d ( i , j ) = 1 ] ⇔ ∑ d ∣ g c d ( i , j ) μ ( d ) [gcd(i,j)=1] \Leftrightarrow \sum\limits_{d|gcd(i,j)}\mu(d) [gcd(i,j)=1]dgcd(i,j)μ(d)

void getMu() {
  mu[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!flg[i]) p[++tot] = i, mu[i] = -1;
    for (int j = 1; j <= tot && i * p[j] <= n; ++j) {
      flg[i * p[j]] = 1;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }
}
莫比乌斯反演有两种形式:

一、
F(n)= ∑ d ∣ n f ( d ) ⇒ f ( n ) = ∑ d ∣ n μ ( d ) F ( n d ) \sum\limits_{d|n}f(d)\Rightarrow f(n)=\sum\limits_{d|n}\mu(d)F(\frac{n}{d}) dnf(d)f(n)=dnμ(d)F(dn)
二、
F(n)= ∑ n ∣ d f ( d ) ⇒ f ( n ) = ∑ n ∣ d μ ( d n ) F ( d ) \sum\limits_{n|d}f(d)\Rightarrow f(n)=\sum\limits_{n|d}\mu(\frac{d}{n})F(d) ndf(d)f(n)=ndμ(nd)F(d)

例题:

1、Luogu 2303
题意:求 ∑ i = 1 N g c d ( i , N ) , ( 1 ≤ N ≤ 2 32 ) \sum\limits_{i=1}^{N}gcd(i,N),(1\leq N\leq2^{32}) i=1Ngcd(i,N),(1N232)
解法: ∑ i = 1 N g c d ( i , N ) = ∑ d ∣ N d ∑ i = 1 N [ g c d ( i , N ) = d ] = ∑ d ∣ N d ∑ i = 1 ⌊ N d ⌋ [ g c d ( i , ⌊ N d ⌋ ) = 1 ] = ∑ d ∣ N d φ ( ⌊ N d ⌋ ) \sum\limits_{i=1}^Ngcd(i,N)=\sum\limits_{d|N}d\sum\limits_{i=1}^N[gcd(i,N)=d]=\sum\limits_{d|N}d\sum\limits_{i=1}^{\lfloor \frac{N}{d}\rfloor}[gcd(i,\lfloor \frac{N}{d}\rfloor)=1]=\sum\limits_{d|N}d\varphi(\lfloor \frac{N}{d}\rfloor) i=1Ngcd(i,N)=dNdi=1N[gcd(i,N)=d]=dNdi=1dN[gcd(i,dN)=1]=dNdφ(dN)
2、求 ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = 1 ] \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(i,j)=1] i=1nj=1m[gcd(i,j)=1]
解法: ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = 1 ] = ∑ i = 1 n ∑ j = 1 m ∑ d ∣ g c d ( i , j ) μ ( d ) = ∑ d = 1 m i n ( n , m ) μ ( d ) ⌊ n d ⌋ ⌊ m d ⌋ \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(i,j)=1]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d|gcd(i,j)}\mu(d)=\sum\limits_{d=1}^{min(n,m)}\mu(d)\lfloor \frac{n}{d}\rfloor\lfloor \frac{m}{d}\rfloor i=1nj=1m[gcd(i,j)=1]=i=1nj=1mdgcd(i,j)μ(d)=d=1min(n,m)μ(d)dndm

杜教筛

问题:求积性函数 f ( i ) f(i) f(i)的前缀和 ∑ i = 1 n f ( i ) \sum\limits_{i=1}^nf(i) i=1nf(i)
构造两个积性函数 h , g h,g h,g,使得 h = f ∗ g h=f*g h=fg
∑ i = 1 n h ( i ) = ∑ i = 1 n ∑ d ∣ i g ( d ) f ( i d ) \sum\limits_{i=1}^n h(i)=\sum\limits_{i=1}^n\sum\limits_{d|i}g(d)f(\frac{i}{d}) i=1nh(i)=i=1ndig(d)f(di)
= ∑ d = 1 n g ( d ) ∑ i = 1 ⌊ n d ⌋ f ( i ) \quad\quad\quad=\sum\limits_{d=1}^{n}g(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i) =d=1ng(d)i=1dnf(i)
∑ i = 1 n h ( i ) = ∑ i = 1 n g ( d ) S ( ⌊ n d ⌋ ) \sum\limits_{i=1}^{n}h(i)=\sum\limits_{i=1}^n g(d)S(\lfloor\frac{n}{d}\rfloor) i=1nh(i)=i=1ng(d)S(dn)
∑ i = 1 n h ( i ) = g ( 1 ) S ( n ) + ∑ d = 2 n g ( d ) S ( ⌊ n d ⌋ ) \sum\limits_{i=1}^{n}h(i)=g(1)S(n)+\sum\limits_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) i=1nh(i)=g(1)S(n)+d=2ng(d)S(dn)
g ( 1 ) S ( n ) = ∑ i = 1 n h ( i ) − ∑ d = 2 n g ( d ) S ( ⌊ n d ⌋ ) g(1)S(n)=\sum\limits_{i=1}^{n}h(i)-\sum\limits_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) g(1)S(n)=i=1nh(i)d=2ng(d)S(dn)
构 造 出 的 h ( i ) 要 使 得 ∑ i = 1 n h ( i ) 能 较 快 的 求 出 来 , 构造出的h(i)要使得\sum\limits_{i=1}^{n}h(i)能较快的求出来, h(i)使i=1nh(i)
先 用 线 性 筛 把 前 m 项 筛 出 来 , 当 m = n 2 3 时 , 复 杂 度 为 O ( n 2 3 ) ; 先用线性筛把前m项筛出来,当m=n^{\frac{2}{3}}时,复杂度为O(n^{\frac{2}{3}}); 线mm=n32O(n32);
当 m = n 1 2 时 , 复 杂 度 为 O ( n 3 4 ) ; 当m=n^{\frac{1}{2}}时,复杂度为O(n^{\frac{3}{4}}); m=n21O(n43);

ll GetSum(int n){//f的前缀和
	if(n<=1000000)return sum[n];//线性筛预处理
	if(mp[n]!=0)return mp[n];//记忆化
	ll ans = f_g_sum(n)%mod; //h的前缀和
	for(ll l = 2, r; l <= n; l = r + 1) { //整除分块
		r = (n / (n / l)); 
		ans = (ans - (g_sum(r) - g_sum(l - 1)) * GetSum(n / l) + mod)%mod;
		//g_sum是g的前缀和,GetSum递归求
	}
	mp[n]=ans;
	return ans; 
}

Min_25筛

前提:一个积性函数 f ( i ) f(i) f(i),且 f ( p ) f(p) f(p)是一个关于 p p p的项数较少的多项式或可以快速求值;即 f ( p c ) f(p^c) f(pc)可以快速求值。
用法:在 O ( n 3 4 log ⁡ n ) O(\frac{n^{\frac{3}{4}}}{\log n}) O(lognn43)的时间复杂度下求 ∑ i = 1 n f ( i ) \sum\limits_{i=1}^nf(i) i=1nf(i)
一些定义:
p k p_k pk:第 k k k个质数,特别地,令 p 0 = 1 p_0=1 p0=1
l p f ( n ) : lpf(n): lpf(n): n n n的最小质因子,特别地,令 l p f ( 1 ) = 1 lpf(1)=1 lpf(1)=1
F p r i m e ( n ) : = ∑ 2 ≤ p ≤ n f ( p ) F_{prime}(n):=\sum_{2\le p\le n}f(p) Fprime(n):=2pnf(p)
F k ( n ) : = ∑ i = 2 n [ p k ≤ l p f ( i ) ] f ( i ) F_k(n):=\sum\limits_{i=2}^{n}[p_k\le lpf(i)]f(i) Fk(n):=i=2n[pklpf(i)]f(i)
根据定义:可以发现 ∑ i = 1 n f ( i ) = F 1 ( n ) + f ( 1 ) \sum\limits_{i=1}^nf(i)=F_1(n)+f(1) i=1nf(i)=F1(n)+f(1)
F k ( n ) = ∑ i = 2 n [ p k ≤ l p f ( i ) ] f ( i ) ① = ∑ k ≤ i , p i ≤ n f ( p i ) + ∑ k ≤ i , p i 2 ≤ n ∑ c ≥ 1 , p i c ≤ n f ( p i c ) ( [ c > 1 ] + F i + 1 ( n p i c ) ) ② = F p r i m e ( n ) − F p r i m e ( p k − 1 ) + ∑ k ≤ i , p i 2 ≤ n ∑ c ≥ 1 , p i c ≤ n f ( p i c ) ( [ c > 1 ] + F i + 1 ( n p i c ) ) ③ = F p r i m e ( n ) − F p r i m e ( p k − 1 ) + ∑ k ≤ i , p i 2 ≤ n ∑ c ≥ 1 , p i c + 1 ≤ n ( f ( p i c ) F i + 1 ( n p i c ) + f ( p i c + 1 ) ) ④ \begin{aligned} F_k(n) &=\sum\limits_{i=2}^{n}[p_k\le lpf(i)]f(i)&①\\ &=\sum\limits_{k\le i,p_i\le n}f(p_i)+\sum\limits_{k\le i,p_i^2\le n}\sum\limits_{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1}(\frac{n}{p_i^c}))& ②\\ &=F_{prime}(n)-F_{prime}(p_{k-1})+\sum\limits_{k\le i,p_i^2\le n}\sum\limits_{c\ge 1,p_i^c\le n}f(p_i^c)([c>1]+F_{i+1}(\frac{n}{p_i^c}))& ③\\ &=F_{prime}(n)-F_{prime}(p_{k-1})+\sum\limits_{k\le i,p_i^2\le n}\sum\limits_{c\ge 1,p_i^{c+1}\le n}(f(p_i^c)F_{i+1}(\frac{n}{p_i^c})+f(p_i^{c+1}))& ④\end{aligned} Fk(n)=i=2n[pklpf(i)]f(i)=ki,pinf(pi)+ki,pi2nc1,picnf(pic)([c>1]+Fi+1(picn))=Fprime(n)Fprime(pk1)+ki,pi2nc1,picnf(pic)([c>1]+Fi+1(picn))=Fprime(n)Fprime(pk1)+ki,pi2nc1,pic+1n(f(pic)Fi+1(picn)+f(pic+1))
② ② 推导:枚举素因子 p i p_i pi,前面部分不解释了,后面部分是 p i p_i pi作为最小素因子时的合数,考虑合数中 p i p_i pi出现了 c c c次, F i + 1 ( n p i c ) F_{i+1}(\frac{n}{p_i^c}) Fi+1(picn)即是剩余部分,因为剩余部分最小的质因子也要大于 p i p_i pi,所以从 p i + 1 p_{i+1} pi+1开始。
③ ③ 推导是傻子问题,不解释了。
④ ④ 推导:对于满足 p i c ≤ n < p i c + 1 p_i^c\le n<p_i^{c+1} picn<pic+1 c c c,有 p i c + 1 > n    ⟺    n p i c < p i < p i + 1 p_i^{c+1}>n\iff \frac{n}{p_i^c}<p_i<p_{i+1} pic+1>npicn<pi<pi+1,所以 F i + 1 ( n p i c ) = 0 F_{i+1}(\frac{n}{p_i^c})=0 Fi+1(picn)=0
然后又有 ∑ c ≥ 1 , p i c + 1 ≤ n f ( p i c + 1 ) = ∑ c ≥ 1 , p i c ≤ n f ( p i c ) [ c > 1 ] \sum\limits_{c\ge 1,p_i^{c+1}\le n}f(p_i^{c+1})=\sum\limits_{c\ge 1,p_i^c\le n}f(p_i^c)[c>1] c1,pic+1nf(pic+1)=c1,picnf(pic)[c>1],其边界值为 F k ( n ) = 0 ( p k > n ) F_k(n)=0(p_k>n) Fk(n)=0(pk>n)
那假设现在已经求出了所有的 F p r i m e ( n ) F_{prime}(n) Fprime(n),那么有两种方式可以求出所有的 F k ( n ) F_k(n) Fk(n)
1.直接按照递推式计算。(一般数据不是特别大的话用这个)
2.从大到小枚举 p p p的转移,仅当 p 2 < n p^2<n p2<n时转移增加值不为0,按照递推式后缀和优化即可。
如何计算 F p r i m e ( n ) F_{prime}(n) Fprime(n)呢?
容易发现 F p r i m e F_{prime} Fprime有且仅有 1 , 2 , … , ⌊ n ⌋ , n n , … , n 2 , n 1,2,…,\lfloor\sqrt{n}\rfloor,\frac{n}{\sqrt{n}},…,\frac{n}{2},n 1,2,,n ,n n,,2n,n处的值是有用的。
Min_25可用的前提是 f ( p ) f(p) f(p)是一个关于 p p p的低次多项式,可以表示为 f ( p ) = ∑ a i p c i f(p)=\sum a_ip^{c_i} f(p)=aipci,那么每个 p c i p_{c_i} pci F p r i m e ( n ) F_{prime}(n) Fprime(n)的贡献是 a i ∑ 2 ≤ p ≤ n p c i a_i\sum_{2\le p\le n}p^{c_i} ai2pnpci
问题就转变成了:给定 n , s , g ( p ) = p s n,s,g(p)=p^s n,s,g(p)=ps
构造一个新的函数 G k ( n ) = ∑ i = 1 n [ p k < l p f ( i ) ∣ i s p r i m e ( i ) ] g ( i ) G_k(n)=\sum\limits_{i=1}^{n}[p_k<lpf(i)|isprime(i)]g(i) Gk(n)=i=1n[pk<lpf(i)isprime(i)]g(i)
对于一个合数 x x x,必定有 l p f ( x ) ≤ x lpf(x)\le \sqrt{x} lpf(x)x ,则 F p r i m e ( n ) = G ⌊ n ⌋ F_{prime}(n)=G_{\lfloor\sqrt{n}\rfloor} Fprime(n)=Gn (n),因为不存在 l p f ( i ) > p ⌊ n ⌋ lpf(i)>p_{\lfloor\sqrt{n}\rfloor} lpf(i)>pn
考虑G的边界值,有 G 0 ( n ) = ∑ i = 2 n g ( i ) , p 0 = 1 G_0(n)=\sum\limits_{i=2}^{n}g(i),p_0=1 G0(n)=i=2ng(i),p0=1
对于转移,有:
对于 n < p k 2 n<p_k^2 n<pk2的部分, G k ( n ) = G k − 1 ( n ) G_k(n)=G_{k-1}(n) Gk(n)=Gk1(n)
对于 p k 2 ≤ n p_k^2\le n pk2n的部分,被筛掉的数必有质因子 p k p_k pk,即 − g ( p k ) G k − 1 ( n p k ) -g(p_k)G_{k-1}(\frac{n}{p_k}) g(pk)Gk1(pkn);但是会有 l p f ( i ) < p k lpf(i)<p_k lpf(i)<pk i i i被减去,应该加回来,即 g ( p k ) G k − 1 ( p k − 1 ) g(p_k)G_{k-1}(p_{k-1}) g(pk)Gk1(pk1)
所以,有: G k ( n ) = G k − 1 ( n ) − [ p k 2 ≤ n ] g ( p k ) ( G k − 1 ( n p k ) − G k − 1 ( p k − 1 ) ) G_k(n)=G_{k-1}(n)-[p_k^2\le n]g(p_k)(G_{k-1}(\frac{n}{p_k})-G_{k-1}(p_{k-1})) Gk(n)=Gk1(n)[pk2n]g(pk)(Gk1(pkn)Gk1(pk1))
代码实现方案:
对于 F p r i m e ( n ) F_{prime}(n) Fprime(n),直接按照递推式实现;
对于 p k 2 ≤ n p_k^2\le n pk2n,可以用线性筛预处理 s k = F p r i m e ( p k ) s_k=F_{prime}(p_k) sk=Fprime(pk)来代替 F k F_k Fk递推式中的 F p r i m e ( p k − 1 ) F_{prime}(p_{k-1}) Fprime(pk1) G G G递推式中的 G k − 1 ( p k − 1 ) = ∑ i = 1 k − 1 g ( p i ) G_{k-1}(p_{k-1})=\sum_{i=1}^{k-1}g(p_i) Gk1(pk1)=i=1k1g(pi)也可以预处理。
使用Min_25时,首先要明确:
如何线性筛出前 n \sqrt{n} n f f f值;
f ( p ) f(p) f(p)的多项式表示;
如何快速求出 f ( p c ) f(p^c) f(pc)
然后按照顺序实现:
1、筛出 [ 1 , n ] [1,\sqrt{n}] [1,n ]内的质数与前 n \sqrt{n} n f f f值;
2、对 f ( p ) f(p) f(p)多项式表示中的每一项筛出对应的 G G G,合并得到 F p r i m e F_{prime} Fprime的所有 O ( n ) O(\sqrt{n}) O(n )个有用点值;
3、按照 F k F_k Fk的递推式实现递归,求出 F 1 ( n ) F_1(n) F1(n)
以下模板以dengtesla大佬的模板进行地修改(完全是因为码风不同,所以改的),大佬模板链接dengtesla大佬的Min_25模板

#include<bits/stdc++.h>
#define ll long long
using namespace std;

const int maxn = 2000;
const int N = 710000;
const int mod = 1e9+7;

int b[maxn+1],c[maxn+1][maxn+1],Inv[maxn+1];
ll sqr,n; /// sqr为sqrt(n)
ll w[N],id1[N],id2[N];
///w[i]记录了形如n/k的第i大的取值是多少

int tot; ///记录对于要筛的n,sqrt(n)以内质数的个数
int isp[N],p[N];
ll zh[N][3]; ///zh[i][k]记录(p[1])^k + (p[2])^k + ... + (p[i])^k
ll g[N][3];
///这个k要根据题目改范围,就是这个3要改成多项式的最高次幂
ll qpow(ll a,ll b){
    ll ans=1;
	a=a%mod;
    while(b){
        if(b&1)
            ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}


ll sigma_f(ll n,int k){ ///得到∑i^k, i:1~n
    if(k==0) return n;
    n++;
    n%=mod;
    ll tmp = n;
    ll ans=0;
    for (int i=1;i<=k+1;i++){
        ans += 1LL*c[k+1][i]*b[k+1-i]%mod*n%mod;
        ans %= mod;
        n *= tmp%mod;
        n %= mod;
    }
    ans = ((ans*Inv[k+1])%mod+mod)%mod;
    return ans;
}

void get_p(int n,int w){
    tot = 0;
    memset(isp,1,sizeof(isp));
    isp[0] = 0;
    isp[1] = 0;
    for(int i=2;i<=n;i++){
        if(isp[i]){
            p[++tot] = i;
            ll wait = 1;
            for(int j=0;j<=w;j++){
                zh[tot][j] = zh[tot-1][j] + wait;
                zh[tot][j] %= mod;
                wait *= i;
            }
        }
        for(int j=1;p[j]*i<=n&&j<=i;j++){
            isp[i*p[j]] = 0;
            if(i%p[j]==0) break;
        }
    }
}

void get_g(ll n,int t){
///对每个x=n/i,求出∑[i是质数](i^t) (i from 1 to x)。每个对应的值存储在g[x][t]中
    int m = 0;
    ll i=1,r;
    while(i<=n){
        ll len = n/i;
        r = n/len;
        if(len<=sqr) id1[len] = ++m;
        else id2[r] = ++m;
        for(int ww=0;ww<=t;ww++){
            g[m][ww] = sigma_f(len,ww)-1;
            g[m][ww] %= mod;
            g[m][ww] += mod;
            g[m][ww] %= mod;
        }
        w[m] = len; ///w[i]记录了形如n/k的第i大的取值是多少
        i = r+1;
    }
    for(int i=1;i<=tot;i++){
        for(int j=1;j<=m;j++){
            if(1LL*p[i]*p[i]>w[j]) break;
            else{
                int op;
                if(w[j]/p[i]<=sqr) op = id1[w[j]/p[i]];
                else op = id2[n/(w[j]/p[i])];
                for(int ww=0;ww<=t;ww++){
                    g[j][ww] = g[j][ww] - qpow(p[i],ww)*((g[op][ww]-zh[i-1][ww])%mod);
                    g[j][ww] %= mod;
                    g[j][ww] += mod;
                    g[j][ww] %= mod;
                }
            }
        }
    }
}

inline ll get_value(int wz,int k){
    ll w = (g[wz][2]+2*g[wz][1]-g[wz][0])-(zh[k-1][2]+2*zh[k-1][1]-zh[k-1][0]);
    w %= mod;
    w += mod;
    w %= mod;
    ///比如f(x) = x - 1,
    ///ll w = (g[wz][1]-g[wz][0])-(zh[k-1][1]-zh[k-1][0]);
    return w; ///自己填写f(x)的表达式(在质数时)
    ///比如f(x) = x^2 + 2*x - 1,
    ///就写(g[wz][2]+2*g[wz][1]-g[wz][0])-(zh[k-1][2]+2*zh[k-1][1]-zh[k-1][0])
}

ll f(ll p,ll k){ ///计算f(p^k)处的值
    if(p==1) return f(1)%mod;//f(1)的值
    else if(k==1) return ...%mod;
    return ...%mod; ///自己填写
}

ll get_s(ll x,int k){
    if(x<=1||p[k]>x) return 0;
    int wz;
    if(x<=sqr) wz = id1[x];
    else wz = id2[n/x];
    ll ans = get_value(wz,k);
    //if(k==1) ans += 2;
    for(int i=k;i<=tot&&1LL*p[i]*p[i]<=x;++i){
        for(ll l=p[i],e=1;l*p[i]<=x;l=l*p[i],++e){
            ans = ans + (get_s(x/l,i+1)*f(p[i],e)%mod)%mod+f(p[i],e+1);
            ans %= mod;
        }
    }
    ans += mod;
    ans %= mod;
    return ans;
}

void init(){
    c[0][0]=1;
    for (int i=1;i<maxn;i++){//求组合数
        for (int j=1;j<=i;j++) 
            c[i][j]=(c[i-1][j-1]+c[i-1][j]) % mod;
        c[i][0]=1;
    }
    Inv[1]=1;//线性求逆元
    for (int i=2;i<maxn;i++) 
        Inv[i]=1LL*Inv[mod % i] * (mod-mod/i) % mod;
    b[0]=1;
    for (int i=1;i<maxn;i++){//伯努利数
        b[i]=0;
        for (int k=0;k<i;k++) 
            b[i]=(b[i]+1LL*c[i+1][k]*b[k] % mod) % mod;
        b[i]=(1LL*b[i]*(-Inv[i+1]) % mod+mod)%mod;
    }
}

void solve(ll n){   
    sqr = sqrt(n);
    get_p(sqr,2);
    get_g(n,2);
    ll ans = get_s(n,1)+f(1,1);
    cout << ans << endl;
}

int main(){
    init();
    while(cin >> n){
        solve(n);
    }
}

代码中的一些说明:
∑ i = 1 n i k \sum\limits_{i=1}^{n}i^k i=1nik时用了伯努利数,即生成函数。
∑ i = 1 n i k = ∑ j = 0 k ( − 1 ) j C k + 1 j B j n k + 1 − j k + 1 \sum\limits_{i=1}^{n}i^k=\frac{\sum\limits_{j=0}^{k}(-1)^jC_{k+1}^{j}B_jn^{k+1-j}}{k+1} i=1nik=k+1j=0k(1)jCk+1jBjnk+1j,其中 B j B_j Bj是伯努利数。
B n = [ n = 0 ] − ∑ k = 0 n − 1 C n + 1 k B k n + 1 B_n=[n=0]-\sum\limits_{k=0}^{n-1}C_{n+1}^{k}\frac{B_k}{n+1} Bn=[n=0]k=0n1Cn+1kn+1Bk,且 B 0 = 1 B_0=1 B0=1

  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值