数论退役帖

acm

数论退役帖

基本知识

最大公因数

找寻a,b的最大公因数, 利用辗转相除法：

1
2
3
4
int gcd(int a, int b){
    return b?gcd(b,a%b):a;
}
// 或者直接使用自带函数 __gcd(a, b)

常用定理

  1. $ gcd(a^n - 1, a^m - 1) = a^{gcd(n, m)} -1$
  2. 扩展 $a>b, gcd(a, b)=1 时 gcd(a^m-b^m, a^n- b^n) = a^{gcd(m, n)} - b^{gcd(m, n)}$
  3. $G = gcd(C^1_n, C^2_n, …, C^{n-1}_n)$
    • n为质数时G=n
    • n有多个质因子时，G=1
    • n有一个质因子p，G=p
  4. $F_n 为斐波那契额数列， gcd(F_n, F_m) = F_{gcd(n, m)}$
  5. 给定两个互素的正整数A,B， 那么他们组合的数C=p*A+q*B>0， C最大为AB-A-B,不能组合的个数为$\frac{(A-1)*(B-1)}{2}$
  6. $(n+1) lcm(C_n^0, C_n^1, .., C^n_n) = lcm(1, 2, 3, .. , n+1)$
  7. 对于质数p,$(x+y+…+w)^p \%p == (x^p + y^p + … + w^p) % p$

素数测定

素数测定常用Miller-Rabin素数测试，是基于费马小定理 $a^p \equiv a(mod \ p ) p为质数$

反过来，如果满足该式子，p大多数为质数

Miller-Rabin 素数测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#define ll long long
using namespace std;
const int Times = 10;

ll multi(ll a, ll b ,ll m){
    ll ans = 0; a%=m;
    while(b){
        if(b&1)ans = (ans+a)%m;
        a = (a+a)%m; b >>= 1;
    } return ans;
}

ll quick_mod(ll a, ll b, ll m){
    ll ans = 1; a %= m;
    while(b){
        if(b&1)ans = multi(ans, a, m);
        a = multi(a,a,m); b>>=1;
    } return ans;
}

bool Miller_Rabin(ll n){
    if(n==2)return true;
    if((n<2) || !(n&1))return false;
    ll m = n-1;
    int k = 0;
    while((m&1)==0){
        k++; m>>=1;
    }
    for(int i = 0; i < Times; i++){
        ll a = rand()%(n-1)+1;
        ll x = quick_mod(a,m,n);
        ll y = 0;
        for(int j = 0; j < k; j++){
            y = multi(x,x,n);
            if(y==1 && x!=1 && x!= n-1)return false;
            x = y;
        } if(y!=1)return false;
    } return true;
}

51Nod 质数检测：

1
2
3
4
5
6
7
8
9
10
11
12
import java.util.Scanner;
import java.math.*;

public class Main {
	public static void main(String[] args) {
		Scanner cin = new Scanner(System.in);
		BigInteger x;
		x=cin.nextBigInteger();
		if(x.isProbablePrime(1)) System.out.println("Yes"); // java 自带Miller_Rabin
		else System.out.println("No");
	}
}

扩展欧几里得

对于二元一次方程 ax+by=gcd(a, b) 利用扩展欧几里得可以求得一组可行解。

1
2
3
4
5
6
7
8
9
void e_gcd(int a, int b, int &d, int &x, int &y){
    if(!b){
        d=a,x=1,y=0;
        return;
    }else{
        e_gcd(b,a%b,d,y,x);
        y -= (a/b)*x;
    }
}

a>b时x为正数。可以令y=(y%a + a) % a 把y变成正数，再利用原始计算x

毕达哥拉斯三元组本原解

对于$x^2 + y ^2 = z^2$ 我们只考虑本原解，即x,y,z互质。

这样的话，我们不妨假设x为偶数，y，z为奇数，那么就存在互质的数n,m

而且n,m一奇一偶且m>n, 有$x=2 m n, y = m^2 - n ^ 2 , z= m^2 + n^2$

POJ - 1305

题意： 给一个正整数n求解，求解n的范围内gcd值为1的勾股数和不是勾股数的个数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include <bits/stdc++.h>
const int N = 1e6+9;
using namespace std;

int biao[N], n;
bool flag[N];

int gcd(int a, int b){return b?gcd(b,a%b):a;}

void init(){
    for(int i = 2; i*i < N; i++){
        for(int j = 1; j < i; j++){
            if(gcd(i,j)!=1)continue;
            if((i&1)&&(j&1))continue;
            if(i*i+j*j<N)biao[i*i+j*j]++; // 标记i^2+j^2即斜边
        }
    } for(int i = 1; i < N; i++) biao[i]+=biao[i-1];
}

int main(){
    init();
    while(~scanf("%d",&n)){
        memset(flag,0,sizeof flag);
        int ans = 0;
        for(int i = 2; i*i < N; i++){
            for(int j = 1; j < i; j++){
                if(gcd(i,j)!=1)continue;
                if((i&1)&&(j&1))continue;
                for(int k = 1; k*(i*i+j*j)<=n; k++){ // 枚举斜边长并标记
                    flag[2*i*j*k]=flag[k*(i*i-j*j)]=flag[k*(i*i+j*j)]=1; // 标记使用过的
                }
            }
        } for(int i = 1; i <= n; i++) if(flag[i])ans++;
        printf("%d %d\n",biao[n],n-ans); 
    } return 0;
}

欧拉函数与欧拉定理

我们定义欧拉函数$\phi (n) = 小于n且与n互质的数， \phi(1) = 1, \phi(2) = 1$

  1. m与n互素， $\phi(mn)=\phi(n)*\phi(m)$
  2. $n = p_1^{a_1} p_2^{a_2} p_3^{a_3}… p_n^{a_n}, \phi(n) = n(1 - \frac{1}{p_1})(1 - \frac{1}{p_2})(1 - \frac{1}{p_3})…*(1 - \frac{1}{p_n}) $
  3. $\sum_{d|n} \phi(d) = n$ 即n的因子的欧拉函数之和为n
  4. a,m互质时, $x \equiv va^{\phi(m) - 1} 等价于 ax \equiv b(mod \ m)$ 事实上就是指数循环节

容斥计算欧拉函数，复杂度$O(\sqrt{n})$：

1
2
3
4
5
6
7
8
9
10
int phi(int n){
    int ans = n;
    for(int i = 2; i*i <= n; i++){
        if(n%i==0){
            ans = ans - ans/i; // 减去i的倍数
            while(n%i==0) n/=i;
        }
    } if(n>1)ans = ans-ans/n; 
    return ans;
}

欧拉筛

欧拉筛是可以做到O(n) 预处理1-n之间所有的质因数， 通过对其修改，我们也能预处理出欧拉函数和莫比乌斯函数以及因数个数。一般欧拉函数不需要修改，可以作为模板.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
const int N = 5e4+5;
// mu为莫比乌斯函数，d为因数个数， phi为欧拉函数， cnt为1-n质数个数
int mu[N],prime[N],d[N],ti[N],phi[N], cnt;
bool vis[N];

void init(){
    mu[1] = d[1] = phi[1] = 1; 
    for(int i = 2 ; i < N; i ++){
        if(!vis[i]){
            prime[cnt++] = i; mu[i] = -1;
            d[i] = 2; ti[i] = 1; phi[i] = i-1;
        } for(int j = 0; j < cnt; j++){
            if(i*prime[j] >= N)break;
            vis[i*prime[j]] = 1;
            if(i%prime[j]==0){
                ti[i*prime[j]] = ti[i] + 1;
                d[i*prime[j]] = d[i]/(ti[i]+1)*(ti[i]+2);
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            } mu[i*prime[j]] = -mu[i];
            d[i*prime[j]] = d[i]*2;
            ti[i*prime[j]] = 1;
            phi[i*prime[j]] = phi[i]*(prime[j]-1);
        }
    }
}

欧拉例题

二次筛法 POJ - 2689

题意: 给出一个区间[L,U]，求给定区间内的质数距离最小的一对和质数距离最大的一对。U<=1e9, U-L <= 1e6

方法: 预处理$\sqrt{1e9} $ 的质数, 然后对[L, U] 区间内去晒(即标记)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <cstdio>
#include <iostream>
#include <cstring>
#define uint unsigned int
const int N = 5e4;
int prime[N], cnt;
bool judge[21*N], vis[N];
using namespace std;

void init(){
    for(int i = 2; i<N; i++){
        if(!vis[i]) prime[cnt++] = i;
        for (int j = 0; j < cnt; j++){
            if(i*prime[j] >= N) break;
            vis[i*prime[j]] = 1;
            if(i%prime[j]==0) break;
        }
    }
}

int main(){
    init(); uint u,v;
    while(~scanf("%d%d",&u,&v)){
        if(u==1)u=2; // 1不是质数
        memset(judge,0,sizeof judge);
        for(int i = 0; i < cnt && prime[i] <= v; i++){
            uint j=((u-1)/prime[i]+1)*prime[i]; // 第一个>=u的prime[i]的倍数
            if(j == prime[i]) j += prime[i]; // j为质数跳过
            for(;j <= v; j+=prime[i]) judge[j-u]=1; // 晒去prime[i]的倍数
        } int st,stlen=N,ed,edlen=0,cot=0,len=0;
        for(int i = 0; i <= v-u; i++){
            if(!judge[i]){  // 即u+i为质数
                if(len && len<stlen){
                    st = i-len; stlen = len; //更新最近的质数对
                } if(len && len>=edlen){
                    ed = i-len; edlen = len; //更新最远的质数对
                } cot++; len = 0; // 统计区间质数个数
            } len ++ ;
        } if(cot>1)printf("%d,%d are closest, %d,%d are most distant.\n",u+st,u+st+stlen,u+ed,u+ed+edlen);
        else puts("There are no adjacent primes.");
    } return 0;
}

HDU - 3939 直角三角形个数

题意: 求a,b,c<=L的三角形个数且满足a<b<c.abc互质. L < 1e12

利用三元组本原解枚举n,m实现sqrt(1e12), 容斥去重

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include <bits/stdc++.h>
#define ll long long
const int N = 1e6+9;
using namespace std;

bool vis[N];
int prime[N/10], phi[N], check[40], cnt, num;
ll ans, L;

void init(){ // 欧拉筛
    for (int i = 2; i < N; i++) {
        if(!vis[i]){prime[cnt++]=i;phi[i]=i-1;}
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j] >= N) break; vis[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);
        }
    }
}

void prime_check(int n){ // 质因数分解
    num = 0;
    if(!vis[n]){ check[num++]=n; return ; } //n 为质数直接返回
    for(int i = 0; i < cnt && prime[i]*prime[i] <= n; i++){
        if(n%prime[i]==0){
            check[num++]=prime[i];
            while(n%prime[i]==0) n/=prime[i];
        }
    } if(n>1) check[num++]=n;
}

//容斥原理计算无关的数fun(j)=j-(j/c1+j/c2+j/c3)+(j/(c1*c2)+j/(c1*c3)+j/(c2*c3))-(j/(c1*c2*c3))。
void dfs(int k, int r, int s, int n){
    if(k==num){
        if(r&1) ans -= n/s;
        else ans+= n/s;
        return;
    } dfs(k+1, r, s, n);
    dfs(k+1, r+1, s*check[k], n);
}

int main(){
    init(); int t;scanf("%d",&t);
    while(t--){
        ans = 0; scanf("%lld",&L);
        int m = (int)sqrt(1.0*L + 0.5);
        for(int i = m; i > 0; i--) { // 枚举m ， n<m 只要计算有多少符合n就行
            int p = (int)sqrt(L - (ll)i*i + 0.5); // n的上界
            if ( i<=p ) {
                if(i&1){
                    prime_check(i); dfs(0,0,1,i>>1);
                } else ans += phi[i]; // 对于偶数，与其互质的数都为奇数可以直接
            } else {
                // i>p时，j所取的范围在[1,p]
                prime_check(i);
                if(i&1)dfs(0,0,1,p>>1);
                else dfs(0,0,1,p);
            }
        } printf("%lld\n",ans);
    }
}

逆元

对于$am \equiv 1 (mod \ m)$, 这个同余方程中最小的x正整数解叫做a模m的逆元;

也就是说$\frac{1}{a} \equiv x (mod \ m)$

由费马小定理得a对于模数mod为质数的逆元$a^{mod - 2}$可以用快速幂log求得

1
2
3
4
5
6
7
8
9
10
11
12
ll qp(ll a, int k){
    ll ans = 1; a %= mod;
    while(k) {
        if(k & 1) ans = ans * a % mod;
        a = a*a % mod;
        k >> = 1;
    } return ans;
}

ll inv(ll a) {
    return qp(a%mod, mod-2);
}

预处理逆元:

我们可预处理1-n的逆元是O(n)的复杂度

1
2
3
inv[1]=1;
for (int i=2;i<=n;++i) 
    inv[i]=(mod-mod/i)*inv[mod%i]%mod;

预处理n! 的逆元inv

1
2
3
4
5
6
7
ll fac[N], inv[N];
void init() {
    fac[0] = inv[0] = 1;
    for (int i = 1; i < N; i++) fac[i]=fac[i-1]*i%mod;
    inv[N-1] = qp(fac[N-1], mod-2);
    for (int i = N-2; i; i--) inv[i]=(i+1)*inv[i+1]%mod;
}

等比数列求和

$S_n = (a + a^2 + .. + a^n) mod M$

如果使用高中的公式那么就是a = 1 时$s_n = n \mod \ M$, 如果a!=1, $S_n = \frac{a - a ^{n+1}}{1 - a}$

这里还有一种方法: 二分

对于n%2 == 0 $S_n = (1 + a ^{\frac{n}{2}}) S_{\frac{n}{2}}$

否则 $S_n = (1 + a^{\frac{n+1}{2}})S_{\frac{n+1}{2}} + a^{\frac{n+1}{2}}$

POJ - 3233 等比矩阵求和

如果利用上面的递推式子就可以解决矩阵的i次求和问题.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
int n,k;
int mod = 1e9+7;

struct Matrix{
    int mp[35][35];
    Matrix(){ memset(mp,0,sizeof mp); }
    Matrix(int v){
        memset(mp,0,sizeof mp);
        for(int i = 1; i <= n; i++){ mp[i][i] = v; }
    }
}a;

inline Matrix multi(Matrix a, Matrix b){
    Matrix c;
    for(int i = 1; i <= n; i++){
        for(int k = 1; k <= n; k++){ if(a.mp[i][k])
            for(int j = 1; j <= n; j++){if(b.mp[k][j])
                c.mp[i][j] = (c.mp[i][j] + a.mp[i][k]*b.mp[k][j])%mod;
            }
        }
    } return c;
}

inline Matrix pls(Matrix a, Matrix b){
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            a.mp[i][j] = (a.mp[i][j]+b.mp[i][j])%mod;
        }
    } return a;
}

inline Matrix powmul(Matrix a, int k){
    Matrix c = Matrix(1);
    while(k){
        if(k&1)c = multi(c,a);
        a = multi(a,a); k>>=1;
    } return c;
}

inline Matrix S(Matrix a, int k){
    if(k==1) return a;
    Matrix temp = Matrix(1);
    if(k&1){
        Matrix temp2 = powmul(a,(k+1)/2);
        return pls(multi(pls(temp, temp2), S(a,(k-1)/2)), temp2);
    } else {
        Matrix temp2 = powmul(a,k/2);
        return multi(pls(temp, temp2), S(a,k/2));
    }
}

int main(){
    while(~scanf("%d%d%d",&n,&k,&mod)){
        for(int i = 1; i <= n; i++){
            for(int j = 1; j <= n; j++){
                scanf("%d",&a.mp[i][j]);
            }
        } a = S(a,k);
        for(int i = 1 ; i <= n; i++){
            for(int j = 1; j < n; j++){
                printf("%d ",a.mp[i][j]);
            } printf("%d\n", a.mp[i][n]);
        }
    } return 0;
}

这个复杂度时有点高的,还可以构造矩阵利用矩阵快速幂一步到位。

我们要求的矩阵和为$S_n$

那么 $$ \begin{bmatrix} S_n \ A^n \end{bmatrix} = \begin{bmatrix}1 & A \ 0 & A \end{bmatrix} \begin{bmatrix} S_{n-1} \ A^{n-1} \end{bmatrix} $$

那么复杂度为O(60^3 * log k)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
int n,k;
int mod = 1e9+7;

struct Matrix{
    int mp[70][70];
    Matrix(){ memset(mp,0,sizeof mp); }
    Matrix(int v){
        memset(mp,0,sizeof mp);
        for(int i = 1; i <= n; i++){ mp[i][i] = v; }
    }
};

inline Matrix multi(Matrix a, Matrix b){
    Matrix c;
    for(int i = 1; i <= n; i++){
        for(int k = 1; k <= n; k++){ if(a.mp[i][k])
            for(int j = 1; j <= n; j++){if(b.mp[k][j])
                c.mp[i][j] = (c.mp[i][j] + a.mp[i][k]*b.mp[k][j])%mod;
            }
        }
    } return c;
}

inline Matrix pls(Matrix a, Matrix b){
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            a.mp[i][j] = (a.mp[i][j]+b.mp[i][j])%mod;
        }
    } return a;
}

inline Matrix powmul(Matrix a, int k){
    Matrix c = Matrix(1);
    while(k){
        if(k&1)c = multi(c,a);
        a = multi(a,a); k>>=1;
    } return c;
}

int main(){
    while(~scanf("%d%d%d",&n,&k,&mod)){
        Matrix a; int g;
        for(int i = 1; i <= n; i++){
            a.mp[i][i] = 1;
            for(int j = 1; j <= n; j++){
                scanf("%d",&g);
                a.mp[i][j+n]=a.mp[i+n][j+n]=g%mod;
            }
        } n *= 2; a = powmul(a, k); n /= 2;
        for(int i = 1; i <= n; i++){
            for(int j = 1+n; j < n*2; j++){
                printf("%d ",(a.mp[i][j])%mod);
            } printf("%d\n", (a.mp[i][2*n])%mod);
        }
    } return 0;
}

组合数计算

从n个物品中选取k个物品有$C_n^{k}$ 种方法。

组合数表预处理 mod任意:

1
2
3
4
5
6
7
void init(){
    c[0][0] = 1;
    for (int i = 1; i < N; i++) {
        c[i][0] = c[i][i] = 1;
        for (int j = 1; j < i; j++) c[i][j] = (c[i-1][j] + c[i-1][j-1]) % mod;
    }
}

预处理阶乘， mod只能为质数

1
2
3
4
5
6
7
8
void init() {
    fac[0] = inv[0] = 1;
    for (int i = 1; i < N; i++) fac[i]=fac[i-1]*i%mod;
    inv[N-1] = qp(fac[N-1], mod-2);
    for (int i = N-2; i; i--) inv[i]=(i+1)*inv[i+1]%mod;
}

ll C(int n, int m) { return fac[n]*inv[m]%mod*inv[n-m]%mod; }

Lucas 模数必须为小质数

1
2
3
4
int Lucas(int n, int m){
    if(m==0)return 1;
    return C(n%p,m%p)*Lucas(n/p,m/p)%p;
}

中国剩余定理(CRT)

有n组同余方程
$$
x \equiv a_1 (mod \ m_1) \ x \equiv a_2 (mod \ m_2) \ x \equiv a_3 (mod \ m_3) \ … \ x \equiv a_n (mod \ m_n) \
$$
利用CRT既可以解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void e_gcd(int a, int b, int &d, int &x, int &y){
    if(!b){
        d = a; x = 1; y = 0;
        return ;
    }
    e_gcd(b, a%b, d, y, x);
    y -= x*(a/b);
}
//中国剩余定理
int CRT(int a[], int m[], int n){
    int M = 1;
    int ans = 0;
    for(int i = 0; i < n; i++)M *= m[i];
    for(int i = 0; i < n; i++){
        int x,y,d;
        int Mi = M/m[i];
        e_gcd(Mi, m[i], d, x, y);
        ans = (ans + Mi*x*a[i])%M;
    }
    if(ans<0)ans+= M;
    return ans;
}

对于普通中国剩余定理要求的$m_1,m_2,…,m_k$互素.但如果发生不互素时,需要采用两两合并.

反素数

f(n)为n的因数个数，对于0<i<n，f(i) < f(n), 则称n为反素数。

反素数一般满足质因数的幂次随着质因数的大小增大而减小

即 $若 n = 2^{t_1}3^{t_2}, 则必有t_1 \ge t_2$

求最小的数使其因数为n CF27E

n<=1000，我们贪心的认为其质因数不会超过前50个质数，而且幂次不会超过63。

dfs暴力枚举加剪枝

1
2
3
4
5
6
7
8
9
void dfs(int depth, ll tmp, int num){  //depth 为第几个质数, tmp为当前的值，um为因数个数
    if(num > n || depth>=cnt)return ;
    if(num == n && ans > tmp){
        ans = tmp; return; // 满足的情况下更新
    } for(int i= 1; i <= 63; i++){ // 枚举幂次， 事实上可以记录上一次的幂次来优化剪枝
        if(ans/prime[depth]<tmp || num*(i+1)>n) return; // 用除法防止溢出
        dfs(depth+1, tmp*=prime[depth], num*(i+1));
    }
}

求1-n因数个数最大的数 URAL 1748

同样使用dfs暴力枚举

1
2
3
4
5
6
7
8
9
void dfs(int depth, ll ansa, int ansb, int limit){ //ansa为值, ansb为因数个数
    if(depth >= tot || ansa > n)return;
    if(ansb > b|| ((ansb == b) && a>ansa)){
        a = ansa; b = ansb;
    } for(ll i = 1; depth < tot && i <= limit; i++){
        if(ansa > n/prime[depth])break;
        dfs(depth+1, ansa*=prime[depth], ansb*(i+1), i);
    } return ;
}

HDU - 4542

题意： X mod ai = bi 对于ai在区间 [1, X] 范围内每个值取一次时，有K个ai使bi等于0 。

告诉你有K个ai使得(不使得) bi等于0

找最小的X

解析：K个相等的时候我们去dfs暴力枚举，K个不相等的话这个数一定不大，先预处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
void init() {
    for(int i = 1; i < N; i++) d[i] = i; 
    for(int i = 1; i < N; i++){
        for(int j = i; j <N; j+=i) d[j]--; // 容斥减去因子个数
        if(!d[d[i]]) d[d[i]] = i; // 更新因子个数为d[i]位置的数为i， 节约了空间,这样也保是最小的
        d[i] = 0
    }
}
void dfs(int depth, ll ansa, int ansb, int limit){
    if(ansb > k || depth > 16)return ;
    if(ansb == k && ansa<ans){
        ans = ansa; return ;
    } for(int i = 1; i <= limit; i++){
        if(ansa > pos/prime[depth] || ansb*(i+1)>k) break;
        ansa*=prime[depth];
        if(k%(ansb*(i+1))==0) dfs(depth+1, ansa, ansb*(i+1), i);
    }
}
int main(){
    init(); int T,t=0; scanf("%d",&T);
    while(t++<T) {
        scanf("%d%d",&type,&k); ans = flag = 0;
        if(type) {
            if(d[k]) printf("Case %d: %d\n",t,d[k]);
            else printf("Case %d: Illegal\n",t);
        } else {
            ans = ~0ull; dfs(0,1,1,63);
            if(ans > inf) printf("Case %d: INF\n",t);
            else printf("Case %d: %llu\n",t,ans);
        }
    } return 0;
}

指数循环节

$$a^b \equiv a^{b \% \phi(c)+\phi(c)}(mod \ c),b>\phi(c) $$

例题 C - What is N?

这里求$$n^{n!} \equiv b(mod \ p)$$有多少个成立的n，对于指数%phi(p)为0后的数，是存在循环节的。只要记录一下一次循环的个数p，最后多余的部分去跑一下记忆化的就知道有多少组了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include <bits/stdc++.h>
#define ll unsigned long long
const int maxn = 1e5+5;
using namespace std;

ll phi(ll n){
    ll res = n;
    for(ll i = 2; i <= sqrt(n); i++){
        if(n%i==0){
            res = res - res/i;
            while(n%i==0) n/=i;
        }
    } if(n > 1) res = res - res/n;
    return res;
}

ll quickmod(ll a, ll b, ll c){
    ll ans = 1;
    while(b) {
        if(b&1) ans = ans*a%c;
        a = a*a%c;
        b >>= 1;
    } return ans;
}

ll b, p, m, fa[maxn];
int main(){
    int t=0,T;scanf("%d",&T);
    while(t++<T){
        scanf("%I64u%I64u%I64u", &b, &p, &m); printf("Case #%d: ",t);
        if(p==1){
            if(m==18446744073709551615ull) printf("18446744073709551616\n"); // 防止溢出
            else printf("%I64u\n",m+1); continue;
        }else{
            ll i = 0, ans = 0, fac = 1;
            ll phic = phi(p);
            for(i = 0; i <= m && fac<=phic; i++) {
                if(quickmod(i, fac, p) == b) ans ++; // 计算n! <= phi(p)的个数
                fac *= i+1;
            } fac = fac%phic; //指数循环节公式
            for(;i <= m && fac; i++){
                if(quickmod(i, fac + phic, p) == b) ans++;
                fac = (fac*(i+1))%phic;
            } // 至多跑一个循环或者不超过M
            if(i <= m){
                ll cnt = 0;
                for(int j = 0; j < p; j++){
                    fa[j] = quickmod(i+j, phic, p); // 记忆化p长度的循环
                    if(fa[j] == b) cnt ++;
                }
                ll idx = (m-i+1)/p; // 循环个数
                ans += cnt * idx;   // 循环的种数
                ll remain = (m-i+1)%p; // 不足一个循环的数
                for(int j = 0; j < remain; j++) if(fa[j] == b) ans ++ ;
            } printf("%I64u\n",ans);
        }
    } return 0;
}

牛客练习赛 22E

题意：给一个长为n的序列，m次操作，每次操作：

1.区间img加img

2.对于区间img，查询imgimg img，一直到$a_r$

请注意每次的模数不同。

区间加值我们可以考虑使用线段树lazy标记，查询就直接使用指数循环节公式倒过来计算就行了。

这道题因为数据太多，用文件读入才能过

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
const int N = 5e5+9;
const int M = 2e7+9;
ll val[N];
bool judge[M];
int prime[M],p[M],mod[N],tot;

namespace fastIO{
#define BUF_SIZE 100000
#define OUT_SIZE 100000
  //fread->read
  bool IOerror=0;
  inline char nc(){
    static char buf[BUF_SIZE],*p1=buf+BUF_SIZE,*pend=buf+BUF_SIZE;
    if (p1==pend){
      p1=buf; pend=buf+fread(buf,1,BUF_SIZE,stdin);
      if (pend==p1){IOerror=1;return -1;}
      //{printf("IO error!\n");system("pause");for (;;);exit(0);}
    }
    return *p1++;
  }
  inline bool blank(char ch){return ch==' '||ch=='\n'||ch=='\r'||ch=='\t';}
  inline void read(int &x){
    bool sign=0; char ch=nc(); x=0;
    for (;blank(ch);ch=nc());
    if (IOerror)return;
    if (ch=='-')sign=1,ch=nc();
    for (;ch>='0'&&ch<='9';ch=nc())x=x*10+ch-'0';
    if (sign)x=-x;
  }
  inline void read(ll &x){
    bool sign=0; char ch=nc(); x=0;
    for (;blank(ch);ch=nc());
    if (IOerror)return;
    if (ch=='-')sign=1,ch=nc();
    for (;ch>='0'&&ch<='9';ch=nc())x=x*10+ch-'0';
    if (sign)x=-x;
  }
  inline void read(double &x){
    bool sign=0; char ch=nc(); x=0;
    for (;blank(ch);ch=nc());
    if (IOerror)return;
    if (ch=='-')sign=1,ch=nc();
    for (;ch>='0'&&ch<='9';ch=nc())x=x*10+ch-'0';
    if (ch=='.'){
      double tmp=1; ch=nc();
      for (;ch>='0'&&ch<='9';ch=nc())tmp/=10.0,x+=tmp*(ch-'0');
    }
    if (sign)x=-x;
  }
  inline void read(char *s){
    char ch=nc();
    for (;blank(ch);ch=nc());
    if (IOerror)return;
    for (;!blank(ch)&&!IOerror;ch=nc())*s++=ch;
    *s=0;
  }
  inline void read(char &c){
    for (c=nc();blank(c);c=nc());
    if (IOerror){c=-1;return;}
  }
  //fwrite->write
  struct Ostream_fwrite{
    char *buf,*p1,*pend;
    Ostream_fwrite(){buf=new char[BUF_SIZE];p1=buf;pend=buf+BUF_SIZE;}
    void out(char ch){
      if (p1==pend){
        fwrite(buf,1,BUF_SIZE,stdout);p1=buf;
      }
      *p1++=ch;
    }
    void print(int x){
      static char s[15],*s1;s1=s;
      if (!x)*s1++='0';if (x<0)out('-'),x=-x;
      while(x)*s1++=x%10+'0',x/=10;
      while(s1--!=s)out(*s1);
    }
    void println(int x){
      static char s[15],*s1;s1=s;
      if (!x)*s1++='0';if (x<0)out('-'),x=-x;
      while(x)*s1++=x%10+'0',x/=10;
      while(s1--!=s)out(*s1); out('\n');
    }
    void print(ll x){
      static char s[25],*s1;s1=s;
      if (!x)*s1++='0';if (x<0)out('-'),x=-x;
      while(x)*s1++=x%10+'0',x/=10;
      while(s1--!=s)out(*s1);
    }
    void println(ll x){
      static char s[25],*s1;s1=s;
      if (!x)*s1++='0';if (x<0)out('-'),x=-x;
      while(x)*s1++=x%10+'0',x/=10;
      while(s1--!=s)out(*s1); out('\n');
    }
    void print(double x,int y){
      static ll mul[]={1,10,100,1000,10000,100000,1000000,10000000,100000000,
        1000000000,10000000000LL,100000000000LL,1000000000000LL,10000000000000LL,
        100000000000000LL,1000000000000000LL,10000000000000000LL,100000000000000000LL};
      if (x<-1e-12)out('-'),x=-x;x*=mul[y];
      ll x1=(ll)floor(x); if (x-floor(x)>=0.5)++x1;
      ll x2=x1/mul[y],x3=x1-x2*mul[y]; print(x2);
      if (y>0){out('.'); for (size_t i=1;i<y&&x3*mul[i]<mul[y];out('0'),++i); print(x3);}
    }
    void println(double x,int y){print(x,y);out('\n');}
    void print(char *s){while (*s)out(*s++);}
    void println(char *s){while (*s)out(*s++);out('\n');}
    void flush(){if (p1!=buf){fwrite(buf,1,p1-buf,stdout);p1=buf;}}
    ~Ostream_fwrite(){flush();}
  }Ostream;
  inline void print(int x){Ostream.print(x);}
  inline void println(int x){Ostream.println(x);}
  inline void print(char x){Ostream.out(x);}
  inline void println(char x){Ostream.out(x);Ostream.out('\n');}
  inline void print(ll x){Ostream.print(x);}
  inline void println(ll x){Ostream.println(x);}
  inline void print(double x,int y){Ostream.print(x,y);}
  inline void println(double x,int y){Ostream.println(x,y);}
  inline void print(char *s){Ostream.print(s);}
  inline void println(char *s){Ostream.println(s);}
  inline void println(){Ostream.out('\n');}
  inline void flush(){Ostream.flush();}
#undef ll
#undef OUT_SIZE
#undef BUF_SIZE
};
using namespace fastIO;

inline void init(){
    p[1] = 1;
    for(int i = 2; i < M; i++){
        if(!judge[i]){ prime[tot++]=i; p[i] = i-1; }
        for(int j = 0; j < tot; j++){
            if(i*prime[j]>=M)break; judge[i*prime[j]] = 1;
            if(i%prime[j] == 0){ p[i*prime[j]] = p[i]*prime[j]; break;} 
            p[i*prime[j]]=p[i]*(prime[j]-1);
        }
    }
}

inline int phi(int n){ return p[n]; }
inline int quickmod(ll a, int b, int c){
    ll ans = 1; a %= c;
    while(b) {
        if(b&1) ans = ans*a%c;
        a = a*a%c; b >>= 1;
    } return ans;
}

inline int slove(ll a, int b, int c){ // 指数循环节公式
    if(b*log(a) >= log(c)) return quickmod(a, b, c)+c;
    else return quickmod(a, b, c);
}
int n,m,op;

// 线段树模板
namespace sgt {
    long long chg[1000010];
    int son[1000010][2],cnt,root;
    void build(int &x,int l,int r) {
        x=++cnt;
        if(l==r) return ;
        int mid=(l+r)>>1;
        build(son[x][0],l,mid);
        build(son[x][1],mid+1,r);
    } void update(int a,int b,int k,int l,int r,long long v) {
        if(a>b || l>r)return;
        if(a<=l && b>=r)chg[k]+=v;
        else {
            int mid=(l+r)>>1;
            if(b<=mid)update(a,b,son[k][0],l,mid,v);
            else if(a>mid)update(a,b,son[k][1],mid+1,r,v);
            else update(a,mid,son[k][0],l,mid,v),update(mid+1,b,son[k][1],mid+1,r,v);
        }
    } long long query(int a,int k,int l,int r,long long v) {
        if(l==r)return v+chg[k]+val[a];
        int mid=(l+r)>>1;
        if(a<=mid)return query(a,son[k][0],l,mid,v+chg[k]);
        else return query(a,son[k][1],mid+1,r,v+chg[k]);
    }
}

int main(){ 
    init(); read(n); read(m);
    for (int i = 1; i <= n; i++) read(val[i]);
    sgt::build(sgt::root,1,n);
    for (int t = 0; t < m; t++) {
        read(op); ll x; int l, r;
        if(op==1){
            read(l); read(r); read(x);
            sgt::update(l, r, sgt::root, 1, n, x); // 区间加值
        } else {
            read(l); read(r); read(x); int ans = 1; mod[l] = x;
            for(int i = l+1; i <= r; i++){
                mod[i] = phi(mod[i-1]); // 预处理幂次的幂次的模数
                if(mod[i] <= 2) r = i; // 某幂次<2那么基本就是0了可以停止了
            } for(int i = r; i >= l; i--){
                ans = slove(sgt::query(i, sgt::root, 1, n, 0), ans, mod[i]); //单点查询，倒着计算
            } printf("%d\n",ans%x);
        }
    } return 0;
}

积性函数

积性函数的定义

  1. 若$f(n)$的定义域为正整数域，值域为复数，即$f:Z+→Cf:Z+→C$，则称f(n)为数论函数。
  2. 若f(n)为数论函数，且f(1)=1，对于互质的正整数p,q有f(p⋅q)=f(p)⋅f(q)，则称其为积性函数。
  3. 若f(n)为积性函数，且对于任意正整数p,q都有f(p⋅q)=f(p)⋅f(q)，则称其为完全积性函数。

常见的积性函数

  1. 除数函数$\sigma_k(n) = \sum_{d|n} d^k$, 表示n的约数的k次幂和
  2. 约数个数函数$d(n) = \sigma_0(n) = \sum_{d|n} 1$ , 表示约数个数
  3. 约数和函数$\sum_{d|n} d$
  4. 欧拉函数$\phi(n)=\sum_{i=1}^{n} [(n, i) == 1]$, 对于n>=2 有 $\phi(n)=\sum_{i=1}^{n} [(n, i) == 1] · i = n\phi(n)/2 $
  5. 莫比乌斯函数$u(n)$, 与恒等函数互为逆元
  6. 元函数 $e(n) = [n == 1]$
  7. 恒等函数$I(n) = 1$
  8. 单位函数$id(n) =n$
  9. 幂函数$id^k(n) = n^k$

两个常用的公式

  1. $e(n) = \sum_{d|n}u(d)$, 可以将u(d)看作是容斥的系数
  2. $n = \sum_{d|n} \phi(d)$

积性函数性质

若f(n)为积性函数，对于正整数$n = p_{1}^{k_1}p_{2}^{k_2}…p_{n}^{k_n}$

有$f(n)=\prod_{i=1}^{t}{f(p_i^{k_i})} $, 对于完全积性函数有$f(n)=\prod_{i=1}^{t}{f(p_i)^{k_i}}$

分块技巧

UVALive - 7426

题意：有一份包含个 bug 的n( 1≤𝑛≤$ 10^6 $)行代码，运一次到崩溃需要的时间为 r( 1≤𝑟≤$ 10^9 $)。你可以任意行添加 printf 语句来输出调试,即你知道是否在执行 printf 语句前就崩溃了。每设置一个 printf语句需要花费 p( 1≤𝑝≤$10^9$)

问最坏的情况下，最少需要多少时间可以定位bug的行数

解释：我们设置f(n)为n行代码debug需要的最少时间。那么对于f(n)我们时可以由前面的状态转移过来的，我们对于n行代码分成2块，那么我们就需要f ( (n+1)/2 ) + r + p的时间， 我们分成k块，那么最大的块的大小就为 (n+k-1)/k 向下取整个， 所以总体时间就是 f((n+k-1)/k) + r + (k-1)*p;

这样我们从1-n递推过来就是O(n^2) 的算法， 我们要考虑优化， 随着k增大，我们会发现有很多(n+k-1)/k 是相同，但是我们应该取相同种k最小的，所以我们反过来枚举(n+k-1)/k的大小，也就是每个块最大为l， 那么最多有(n+l-1)/l 块， 所以这样总体时间就是f(l) + r +( (n+l-1)/l )* p;

对于n√n的算法还是很吃力的，我们考虑使用记忆化搜索来降低复杂度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include <bits/stdc++.h>
using ll = long long;
using namespace std;
const int N = 1e6+9;
int n, r, p, vis[N], t;
ll f[N];
ll dfs(int n){
    if(vis[n] == t || n == 1) return f[n];  // 如果已经计算过直接返回
    ll ans = 1e18+7; vis[n] = t; // 表示t组样例的时候f[n]被跟新了
    for(int j = 1; j <= sqrt(n+0.5); j++) {
        ans = min(ans, dfs((n+j)/(j+1)) + r + 1LL*j*p); // 枚举j+1块
        ans = min(ans, dfs(j)+r+((n+j-1)/j-1LL)*p);     // 枚举每块最大为j个
    } return f[n] = ans; // 记忆化
}

int main(){
    while(~scanf("%d%d%d", &n, &r, &p)) {
        t++; printf("%lld\n", dfs(n));   //每组，记忆话搜索一下
    } return 0;
}

迪利克雷卷积

设f,g为两个数论函数，则满足$h(n) = \sum_{d|n} f(d)g(\frac{n}{d})$ 称为f与g的迪利克雷卷积。

积性函数例题 HDU - 5528

题意：f(n)=为[0, n)种选2个数a,b且ab不为n的倍数的方案数。

求$g(n) = \sum_{d|n} f(d) \% 2^{64}$

$1 \le T \le 2e4, 1 \le n \le 1e9$

解析：我们先考虑f(n) ， f(n) = n*n - ab为n的倍数的方案数

即$f(n) = n^2 - \sum_{d|n}\phi(\frac{n}{d})d$

令$h(n) = \sum_{d|n}\phi(\frac{n}{d})d$

那么原式=$\sum_{d|n}d^2 - \sum_{d|n}h(d)$, 很显然这两个函数都是积性函数，我们如果直接对其质因数分解是可以做的。不过处理后面的还是很麻烦。

$\sum_{d|n}\sum_{w|d}\phi(w) * \frac{d}{w} = \sum_{w|n}w· \sum_{\frac{d}{w} |\frac{n}{w} } \phi(\frac{d}{w}) = \sum_{w|n} n = n·d(n)$

这样处理就方便了很多

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include <bits/stdc++.h>
const int N = 1e5+8;
using namespace std;
using ll = unsigned long long;
using lll = unsigned __int128;
ll a[N], num[N], cnt, prime[N], tot, n;
bool vis[N];

void init(){ // 欧拉筛
    for(int i = 2; i < N; i++) {
        if(!vis[i]) prime[tot++] = i;
        for (int j = 0; j < tot; j++) {
            if(i*prime[j]>=N) break;
            vis[i*prime[j]] = 1;
            if(i%prime[j]==0) break;
        }
    }
}

void d(int n){ // 质因数分解
    #define j prime[i]
    for(int i = 0; 1LL*j*j <= n && i < tot; i++){
        if(n%j==0) {
            num[cnt] = 0; a[cnt] = j;
            while(n%j==0) {
                num[cnt]++; n /= j;
            } cnt++;
        }
    } if(n>1) a[cnt]=n,num[cnt++]=1;
    #undef j
}

lll qp(lll a, int k){ // 快速幂
    lll ans = 1;
    while(k) {
        if(k&1) ans *= a;
        a *= a; k >>= 1;
    } return ans;
}

ll p(ll n){
    if(n == 1) return 1; ll ans = 1;
    for (int i = 0; i < cnt; i++) ans *= (qp(a[i], num[i]*2+2)-1)/(a[i]*a[i]-1);
    return ans;
}

ll q(ll n) {
    ll ans = n;
    for (int i = 0; i < cnt; i++) ans *= num[i] + 1;
    return ans;
}

int main(){
    int t;scanf("%d", &t); init();
    while(t--) {
        scanf("%lld", &n); cnt = 0; d(n);
        printf("%llu\n", p(n) - q(n));
    } return 0;
}

积性函数习题 757E

题意: $f_0 (n)$ 为 (p, q) == 1, 且p*q = n的个数, $f_{r+1}(n) = \sum_{u * v=n} \frac{f_{r}(u) + f_{r}(v)}{2} $

求$f_r(n) \%1e9+7$

解释： 显然$f_0(n) = 2^{n的质因子个数}$，这是一个积性函数，我们可以对其化简$f_0(n) = f_0(p_1^{k_1})f_0(p_2^{k_2})…f_0(p_t^{k_t}) = \prod_{i=1}^t(1+k_i)$

f(p^k)只与k有关，而k又是logn的，由$f_r(n)$定义可知这是一个积性函数，所以做r次前缀和就可以了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
const int N = 1e6+9;
ll f[N][21];
const int mod = 1e9+7;
int prime[N], cnt;
bool vis[N];

void init(){
    for (int i = 2; i < N; i++) {
        if(!vis[i]) prime[cnt++] = i;
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j] >= N) break; vis[i*prime[j]] = 1;
            if(i%prime[j]==0) break;
        }
    } f[0][0] = 1;
    for(int i = 1; i < 21; i++) f[0][i] = 2;
    for(int i = 1; i < N; i++) {
        f[i][0] = 1;
        for (int j = 1; j < 21; j++) f[i][j] = (f[i][j-1]+f[i-1][j])%mod; //前缀和
    }
}
int q, r, n;

ll d(int n){ // 质因数分解
    #define j prime[i]
    ll ans = 1;
    for (int i = 0; i < cnt && j*j <= n; i++) {
        if(n%j==0) {
            int t = 0;
            while(n%j==0)n/=j,t++; // 算p^k
            ans = (ans * f[r][t])%mod;  //计算每个质因数的贡献
        }
    } if(n>1) ans = (ans * f[r][1]) % mod; 
    return ans;
    #undef j
}

int main() {
    init(); scanf("%d", &q);
    while(q--) {
        scanf("%d%d", &r, &n);
        printf("%lld\n", d(n));
    } return 0;
}

莫比乌斯反演

$$
f(n) = \sum_{d|n} g(d) <=> g(n) = \sum_{d|n}u(d)f(\frac{n}{d})
$$

这就是莫比乌斯函数，对于已知的f(n)，我们都可以利用该等式来反过来求g(n);

经典例题-求互质的数对数 HDU - 1695

题意：求1<=i<=b, 1<=j<=d 由多少对(i,j)满足gcd(i, j) = k

解析：我们先解决1<=i<=n, 1<=j<=m. gcd(i, j)== 1的问题

令F(d) 为有多少对i，j满足gcd(i, j) == d的倍数

f(d) 为有多少对gcd(i, j) == d， 那么我们就有$F(d) = \sum _{}f(id)$, $f(1) = \mu(1)F(1) + \mu(2)F(2)+…$

显然$f(1) = \sum_{i=1}^{min(n, m)} \mu(i) F(i)= \sum_{i=1}^{min(n, m)} \mu(i) \frac{n}{i} \frac{m}{i}$,

利用分块就可以做到$\sqrt{n}+\sqrt{m}$ 查询

所以对于原问题，b,d/=k就转成了gcd(i,j)==1, 再去减去不符合的就可以了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <bits/stdc++.h>
#define ll long long
const int N = 1e5+10;
using namespace std;
bool vis[N];
int mu[N], prime[N], cnt;
inline void init(){ // 欧拉筛
    mu[1] = 1;
    for(int i = 2; i < N; i++){
        if(!vis[i]){ prime[cnt++] = i; mu[i] = -1; }
        for(int j = 0; j < cnt; j++){
            if(i*prime[j]>N)break;
            vis[i*prime[j]]= 1;
            if(i%prime[j] == 0){ mu[i*prime[j]] = 0; break; }
            mu[i*prime[j]] = -mu[i];
        }
    } for (int i = 2; i < N; i++) mu[i] += mu[i-1];
}

inline ll slove(int l, int r){ // 分块
    ll ans = 0;
    for(int i = 1, tmp; i <= l; i=tmp+1) {
        tmp = min(l/(l/i), r/(r/i));
        ans += (ll)(mu[tmp]-mu[i-1])*(l/i)*(r/i);
    } return ans;
}

int main(){
    int b,d,k,T=0, t; scanf("%d",&t); init();
    while(t--) {
        scanf("%*d%d%*d%d%d",&b,&d,&k);
        if(!k) {
            printf("Case %d: 0\n",++T); continue;
        } b/=k;d/=k; if(b > d) swap(b, d);
        printf("Case %d: %lld\n",++T,slove(b,d) - slove(b, b)/2);
    } return 0;
}

GCD表中的质数

题意：有一个M * N的表格，行与列分别是1 - M和1 - N，格子中间写着行与列的最大公约数Gcd(i, j)(1 <= i <= M, 1 <= j <= N)。 求质数个数。

解析：我们枚举p,利用上面的结论就有$ans = \sum_{p}\sum_{i=1}^{\frac{n}{p}} \lfloor \frac{n}{ip} \rfloor \lfloor \frac{m}{ip} \rfloor$ , 根据这个公式的复杂度$O(T (\frac{n}{log\ n}·(\sqrt{n}+\sqrt{m}) ) )$ , 这样复杂度有点高，我们对原始做一些优化，交换两个求和符号就可以得到$\sum_{i=1}^{n}\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor \sum_{j|i}\mu(\frac{i}{j}) \&\& j 为质数$

如果我们预处后面的一个求和，那么就可以做到分块查询

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include <bits/stdc++.h>
using namespace std;
const int N = 5e6+9;
using ll = long long;
int prime[N/5], mu[N], cnt, n, m;
ll f[N]; // 为后一个求和
bool vis[N];

inline void init(){
    mu[1] = 1;
    for (int i = 2; i < N; i++) {
        if(!vis[i]) {prime[cnt++] = i; mu[i]=-1; f[i]=1;} // 对于质数f[i] = mu[1] = 1;
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j]>=N) break;
            vis[i*prime[j]] = 1;
            if(i%prime[j]==0) { f[i*prime[j]] = mu[i]; break; } //prime[j]次数为2时，之前的u都变成了0，只剩下p=prime[j],所以f[i*prime[j]]=mu[i];
            f[i*prime[j]] = -f[i] + mu[i]; mu[i*prime[j]] = -mu[i]; //次数为1时，u变成原来的相反数，再多上新增的。
        }
    } for (int i = 2; i < N; i++) f[i] += f[i-1]; // 前缀和
}

inline ll solve(int n, int m) { // 分块查询
    ll ans = 0;
    for (int i = 1, j; i <= n; i=j+1) {
        j = min(n/(n/i), m/(m/i));
        ans += (f[j]-f[i-1])*(n/i)*(m/i);
    } return ans;
}

int main(){
    int T; scanf("%d", &T); init();
    while(T--){
        scanf("%d%d", &n, &m); if(n > m) swap(n, m);
        printf("%lld\n", solve(n, m));
    }
}

诸如此类的还有

  1. $\sum_{a=1}^{N} \sum_{b=1}^M gcd(a, b) = \sum_{d=1}^{min(N, M)} \phi(d) \lfloor \frac{N}{d} \rfloor \lfloor \frac{M}{d} \rfloor$
  2. $\sum_{a=1}^{N} \sum_{b=1}^M lcm(a, b) = \frac{1}{4}\sum_{d=1}^{min(N, M)} \lfloor \frac{N}{d} \rfloor \lfloor \frac{M}{d} \rfloor(\lfloor \frac{N}{d} \rfloor + 1)(\lfloor \frac{M}{d} \rfloor+1) d \sum_{d’|d}d’u(d’) $
  3. $\sum_{a=1}^{N} \sum_{b=1}^N lcm(a, b) = \sum_{i=1}^N (-i+2\sum_{j=1}^ilcm(i,j))$, 预处理后可以O(1)输出
  4. $\sum_{i=1}^n gcd(i, n) = \sum_{d|n}\phi(d) \frac{n}{d}$
  5. $\sum_{i=1}^n e(gcd(i, n)) = \sum_{d|n}\mu(d) \frac{n}{d}$
  6. 约束和之和$\sum_{i=1}^n = \sum_{i=1}^n i \lfloor \frac{n}{i} \rfloor$
  7. $g(n)=\frac{n}{2}(\phi(n)+e(n))$
  8. $\sum_{i=1}^n lcm(i, n) = n\sum_{d|n}g(\frac{n}{d})$
  9. $\sum_{i=1}^n\sum_{j=1}^mij = \frac{n(n+1)m(m+1)}{4}$
  10. $\sum_{i=1}^n\sum_{j=1}^mij · e(gcd(i,j)) = \sum_{d=1}^n \mu(d) \sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}ij$

杜教筛

设$S(n) = \sum_{i=1}^n f(i)$, 若存在数论函数g,设$h=f*g$, 则有$\sum_{i=1}^{n}h(i)=\sum_{i=1}^ng(i)S(\lfloor \frac{n}{i} \rfloor)$ , 化简得

$g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)S(\lfloor \frac{n}{i} \rfloor)$,

如果我们可以$O(sqrt{n})$计算 $\sum_{i=1}^n h(i)$, O(1)计算g(n)前缀和，那么我们就可以分块查询S(n), 复杂度为$O(n^{\frac{3}{4}})$, 如果预处理前$n^{\frac{3}{4}}$ 项就可以优化到$O(n^{\frac{3}{4}})$,

欧拉函数前缀和 51nod 1239

题意：求$\sum_{i=1}^n \phi(n)$ , n < 1e11;

解析：$S（n）= \sum_{i=1}^n i - \sum_{i=2}^n 1 *S(\lfloor \frac{n}{i} \rfloor)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
const int N = 1.5e7+9;
const int mod = 1e9+7;
int prime[N], phi[N], cnt;
bool vis[N];
unordered_map<ll, ll> mp; // 记忆化

void init(){ // 欧拉筛
    phi[1] = 1;
    for (int i = 2; i < N; i++) {
        if(!vis[i]) {
            prime[cnt++] = i; phi[i] = i - 1;
        } for (int j = 0; j < cnt; j++) {
            if(i*prime[j] >= N) break;
            vis[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);
        }
    } for (int i = 2; i < N; i++) phi[i] = (phi[i-1] + phi[i])%mod;
}

ll qp(ll a, int k){
    ll ans = 1;
    while(k) {
        if(k & 1) ans = ans * a % mod;
        a = a* a %mod; k >>= 1;
    } return ans;
}ll inv2 = qp(2, mod-2);

ll dfs(ll n) {
    if(n < N) return phi[n]; // 预处理
    if(mp.find(n) != mp.end()) return mp[n];
    ll ans = n%mod*(n%mod+1)%mod*inv2%mod; // n*(n+1)/2
    for (ll i = 2,j; i <= n; i=j+1) { // 分块
        j = n/(n/i);
        ans = (ans - (j-i+1)%mod*dfs(n/i)%mod) %mod; // 1 从i加到j = (j-i+1)
    } return mp[n]=(ans%mod + mod)%mod;
}

ll n;
int main() {
    init(); scanf("%lld", &n);
    return printf("%lld\n", dfs(n)), 0;
}

莫比乌斯函数前缀和

题意：求$\sum_{i=1}^n \mu(i)$, i <= 1e10

解析：u*I=e, 所以$S(n) = 1 - \sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor )$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
ll a, b;
const int N = 1e7+9;
int prime[N], cnt;ll mu[N];
bool vis[N];
unordered_map<ll, ll> mp;

void init(){ // 欧拉筛
    mu[1] = 1;
    for (int i = 2; i < N; i++) {
        if(!vis[i]) { prime[cnt++] = i; mu[i] = -1; }
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j] >= N) break; vis[i*prime[j]] = 1;
            if(i%prime[j]==0) break; mu[i*prime[j]] = - mu[i];
        }
    } for (int i = 2; i < N; i++) mu[i] += mu[i-1];
}

ll dfs(ll n){
    if(n < N) return mu[n];
    if(mp.count(n)) return mp[n];
    ll ans = 1; // e
    for (ll i = 2,j; i <= n; i=j+1){
        j = n/(n/i); ans -= dfs(n/i)*(j-i+1); // 和上一题一样
    } return mp[n]=ans;
}

int main() {
    scanf("%lld%lld", &a, &b); init();
    return printf("%lld\n", dfs(b)-dfs(a-1)), 0;
}

一道习题 ：

$N^2-3N+2=\sum_{d|N} f(d)$, 求$\sum_{i=1}^Nf(i)\ mod \ 1e9+7, N \le 1e9$

这道题很容易用杜教筛做，不过他的前1e7前缀预处理可以用容斥来写

1
2
3
4
5
6
7
8
inline void init(){
    for(int i = 1; i < N; i++){
        f[i] = (i-1LL)*(i-2LL)%mod;
    } for (int i = 1; i < N; i++) {
        for (int j = i+i; j < N; j+=i)
            f[j] = ((f[j]-f[i])+mod)%mod; // 容斥
    } for (int i = 2; i < N; i++) f[i]=(f[i]+f[i-1])%mod; // 前缀和
}

Extended Eratosthenes Sieve洲阁筛

洲阁筛时一种在$O(\frac{n^{\frac{3}{4}}}{log \ n})$ 时间内计算大多数积性函数的前缀和的方法。

F(x)是一个积性函数，要求在低于线性的时间内求出。当p为质数时，$F(p^c)$是关于p的低阶多项式。

对于1-n中的所有数，我们按照是否含有>$\sqrt{n}$ 的质因子分为两类，则显然有
$$
\sum_{i=1}^n F(i) = \sum_{i=1 \&\& i没有大于\sqrt{n} 的质因数}^n F(i)(1 + \sum_{\sqrt{n}\le j \le \lfloor \frac{n}{i} \rfloor \& \& j为质数} F(j)\ \ \ )
$$

事实上，我们可以预处理 $ 1 \le i < \sqrt{n} $ 的F, 那么现在的问题就是要解决

  1. $\sum_{\sqrt{n} < j < \le \lfloor \frac{n}{i} \&\&j为质数\rfloor} F(j)$

我们令g(i, j) 表示[1, j]中与前i各质数互质的数的k次幂和。显然有
$$
g(i, j) = g(i-1, j)-p^k_ig(i-1, \lfloor \frac{j}{p_i}\rfloor) \
p_i^2 > j时 g(i, j) = g(i-1, j) - p^k_i
$$

  1. $$\sum_{\sqrt{n} \le i \le n \&\& i 没有质因数>\sqrt{n}} F(i)$$

设f(i, j) 表示[1, j]中仅由i个质数个组成的数的F(x)之和
$$
f(i, j) = f(i-1, j) + \sum_{c \ge 1}F(p^c_i)f(i-1, \lfloor \frac{j}{p_i^c} \rfloor) \
当p^2_i > j 时， f(i, j) = f(i-1, j) + F(p_i)
$$

素数个数

题意：求1-n中素数个数，n<1e11;

解析：对于这个虽然不是积性函数，但是我们依然可以分为两块1-$\sqrt{n}$ 和$(\sqrt{n}+1 ) - n$

前者可以O(n)预处理，后者就可以利用上面第一块的$g(\sqrt{n},n)$递推

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include <bits/stdc++.h>

const int N = 4e6+9;
using ll = long long;
using namespace std;
int prime[N], cnt, res[N];
bool vis[N];

void init() { //欧拉筛
    for (int i = 2; i < N; i++) {
        res[i]=res[i-1]; if(!vis[i]) {prime[cnt++]=i; res[i]++;}
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j]>=N) break; vis[i*prime[j]] = 1;
            if(i%prime[j]==0) break;
        }
    }
}

ll n, v[N], k, last[N], g[N];
int solve(){
    if(n < N) return printf("%d\n", res[n]), 0;
    int tot = 0, sn = sqrt(n+0.5); // sqrt(n)
    int pos = upper_bound(prime, prime+cnt, sn) - prime; // 第一个大于sqrt(n)的质数
    for (ll i = n; i >= 1; i = n/(n/i+1) ) v[++tot] = n/i; // 分块
    for (int i = 1; i <= tot; i++) g[i]=v[i], last[i] = 0; // 初始化
    for (int i = 0; i < pos; i++) {
        for (int j = tot; j; j--) {
            k = v[j]/prime[i]; if(k < prime[i]) break; // 忽略j<p^2
            k = k < sn? k: tot - n/k + 1; // 找到在v中的下标
            g[j] -= g[k] - (i - last[k]); // 减去g[k]中已经减过的
            last[j] = i + 1; // 上一次处理的是第i+1个质数(从1开始)
        }
    } printf("%lld\n", res[sn]*1LL + g[tot] - 1);
}

int main(){
    init(); scanf("%lld", &n); return solve(), 0;
}

洲阁筛实现起来较为复杂，现在以及被Min_25筛替代。

另外一些题

题意：定义w(n) = n的质因子个数，$g(n)=2^{w(n)}$, 求$S(n) = \sum_{i=1}^{n}g(i)\ mod\ 1e9+7$

CCPC 杭州2016J

$$
\begin{eqnarray} ans&=&\sum_{i=1}^ng(i)\ &=&\sum_{i=1}^n\sum_{d|i}\mu^2(d)\ &=&\sum_{i=1}^n\sum_{d|i}\sum_{k^2|d}\mu(k)\ &=&\sum_{k=1}^n\mu(k)\sum_{k^2|d}\lfloor\frac{n}{d}\rfloor\ &=&\sum_{k=1}^n\mu(k)\sum_{i=1}^{\lfloor\frac{n}{k^2}\rfloor}\lfloor\frac{n}{k^2i}\rfloor\ &=&\sum_{k=1}^{\sqrt{n}}\mu(k)S(\lfloor\frac{n}{k^2}\rfloor) \end{eqnarray}
$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <bits/stdc++.h>
const int N = 1e6+9;
const int mod = 1e9+7;
using namespace std;
using ll = long long;
int prime[N], cnt, mu[N];
bool vis[N];

ll n;
inline void init(){
    mu[1] = 1;
    for (int i = 2; i < N; i++) {
        if(!vis[i]) {prime[cnt++] = i; mu[i]=-1;}
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j] >= N) break; vis[i*prime[j]]=1;
            if(i%prime[j] == 0) break; mu[i*prime[j]] = -mu[i];
        }
    }
}
int f[N];
inline int S(ll n) {
    if(n < N && f[n]) return f[n];
    ll ans = 0;
    for (ll i = 1, j; i <= n; i = j+1) {
        j = n/(n/i); ans = (ans + (n/i)*(j-i+1))%mod;
    } if(n < N) f[n] = ans;
    return ans;
}

int main() {
    int T, ks=1; scanf("%d", &T); init();
    while(T--) {
        scanf("%lld", &n);
        ll ans = 0;
        for(ll k = 1; k * k <= n; k++) {
            if(mu[k]) ans = (ans + mu[k]*S(n/k/k)) % mod;
        } printf("Case #%d: %lld\n",ks++, (ans+mod)%mod);
    }
}

2018 四川省赛

题意：求$\sum_{i=1}^n\sum_{j=1}^i n \%(ij), n \le 1e11$

原式=$\sum_{i=1}^n\sum_{j=1}^i n - \lfloor \frac{n}{ij} \rfloor = nn(n+1)/2 -\sum_{i=1}^n\sum_{j=1}^i \lfloor \frac{n}{ij} \rfloor$

= $nn(n+1)/2 -1/2\sum_{i=1}^n\sum_{j=1}^n \lfloor \frac{n}{ij} \rfloor - 1/2\sum_{i=1}^n \lfloor \frac{n}{i^2} \rfloor $

欧拉筛预处理前n^{2/3}项后其他的暴力算，总复杂度是O(n^{2/3})

由于最后答案较大需要使用大数或者__int128

预处理的话，可以考虑$f[n] = \sum_{i=1}^i \lfloor \frac{n}{i}\rfloor i$, 那么$f[n]-f[n-1] = \sum_{i=1}^n i*( \lfloor \frac{n}{i}\rfloor - \lfloor \frac{n-1}{i}\rfloor) = \sum_{i|n}i$, 这个在欧拉筛中可以预处理,然后求一遍前缀和

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <bits/stdc++.h>
#define ll long long
typedef __int128 dll;
const int N = 3e7;
using namespace std;
 
ll n;
int prime[N], cnt, t[N];
dll f[N];
bool vis[N];
 
void init(){ // 欧拉筛
    f[1] = t[1] = 1;
    for (int i = 2; i < N; i++) {
        if(!vis[i]) {prime[cnt++]=i; f[i]=i+1; t[i]=i;}
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j]>=N) break; vis[i*prime[j]] = 1;
            if(i%prime[j] == 0) {
                f[i*prime[j]] = f[i/t[i]] + f[i]*prime[j];
                t[i*prime[j]] = t[i] * prime[j]; break;
            } f[i*prime[j]] = f[i]*(prime[j]+1); t[i*prime[j]] = prime[j];
        }
    } for (int i = 2; i < N; i++) f[i] += f[i-1];
}
 
inline dll F(ll n) { // 分块求\sum_ n/i * i
    if(n < N) return f[n];
    dll ans = 0;
    for (ll i = 1,j; i <= n; i=j+1) {
        j = n/(n/i);
        ans += ((dll)(j-i+1))*(i+j)/2*(n/i);
    } return ans;
}
 
void print(dll x){ // int128位数输出
    if(x > 9) print(x/10);
    putchar(x%10 + '0');
}
 
int main() {
    int t; scanf("%d", &t); init();
    while(t--) {
        scanf("%lld", &n);
        dll ans = n; ans *= n; ans += ans * n; // ans = n*n*(n+1);
        for (int i = 1; 1LL * i * i <= n; i++) ans -= n/i/i*i*i;
        for (ll i = 1,j; i <= n; i=j+1) {
            j=n/(n/i);
            ans -= ((dll)(j-i+1))*(i+j)/2*F(n/i);
        } print(ans/2); puts("");
    } return 0;
}

rng_58-clj等式

$$
\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^cd(ijk) = \sum_{gcd(i,j)=gcd(j,k)=gcd(i,k)=1} \lfloor \frac{a}{i} \rfloor \lfloor \frac{b}{j} \rfloor \lfloor \frac{c}{k} \rfloor
$$

事实上这个等式可以扩展至任意维。

题意：求$\sum_{i=1}^a \sum_{j=1}^b$d(ij), T, a, b < 5e4;

原式化简为$\sum_{gcd(i,j)=1}\lfloor \frac{a}{i} \rfloor \lfloor \frac{b}{j} \rfloor = \sum_{g=1}^{min(a,b)} \mu(g) \sum_{i=1}^{a/g} \sum_{j=1}^{b/g}\lfloor \frac{a}{i} \rfloor \lfloor \frac{b}{j} \rfloor$

分块查询即可,代码略

题意：$\sum\limits_{i=1}^a\sum\limits_{j=1}^b\sum\limits_{k=1}^c d(ijk)$, a,b,c < 2000

化简：$\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{k=1}^t {\lfloor {n\over i}\rfloor}{\lfloor {m\over j}\rfloor}{\lfloor {t\over k}\rfloor} [(i,j)=1][(j,k)=1][(i,k)=1]$

=>$\sum\limits_{k=1}^t {\lfloor {t\over k}\rfloor} \sum\limits_{d=1}^n \mu(d) \sum\limits_{i=1}^ {\lfloor {n\over i}\rfloor} {\lfloor {n\over id}\rfloor}[(di,k)=1] \sum\limits_{j=1}^ {\lfloor {m\over j}\rfloor} {\lfloor {m\over id}\rfloor} [(dj,k)=1]$

暴力算一下就行了(gcd可以记忆化一下)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include <bits/stdc++.h>
#define gcd __gcd
const int N = 2e3+5;
const int mod = (1<<30) - 1;
using namespace std;

int mu[N], prime[N], a, b, c, gd[N][N], cnt;
bool vis[N];

inline int cal(int d, int x){
    int ans = d;
    for(int i = 2; i <= d; i++) if(gd[i][x] == 1) ans += d/i;
    return ans;
}

inline void init(){
    mu[1] = 1;
    for(int i = 2; i < N; i++) {
        if(!vis[i]) { prime[cnt++] = i; mu[i] = -1; }
        for(int j = 0; j < cnt; j++) {
            if(i*prime[j] >= N) break; vis[i*prime[j]] = 1;
            if(i%prime[j] == 0) break; mu[i*prime[j]] = - mu[i];
        }
    } for (int i = 1; i < N; i++) {
        for (int j = 1; j <= i; j++) gd[i][j]=gd[j][i]=gcd(i,j);
    }
}

int main(){
    scanf("%d%d%d",&a,&b,&c); init();
    int ans = 0; if(b > c) swap(b, c);
    for(int i = 1; i <= a; i++) { for(int j = 1; j <= b; j++)
        if(gcd(i, j) == 1 && mu[j]) ans += (a/i)*mu[j]*cal(b/j, i)*cal(c/j, i);
    } return printf("%u\n",ans&mod), 0;
}

BZOJ4176

题意：$\sum_{i=1}^n \sum_{j=1}^nd(ij), n \le 1e9$

化简得$= \sum_{g=1}^{n} \mu(g) \sum_{i=1}^{n/g} \sum_{j=1}^{n/g}\lfloor \frac{n}{i} \rfloor \lfloor \frac{n}{j} \rfloor$

分块查询的难调在于u(g)的前缀和，这个使用杜教筛就可以解决。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int N = 1e7+9;
const int mod = 1e9+7;
int prime[N], cnt;ll mu[N];
bool vis[N];
unordered_map<ll, ll> mp;
ll n;

void init(){ // 欧拉筛
    mu[1] = 1;
    for (int i = 2; i < N; i++) {
        if(!vis[i]) { prime[cnt++] = i; mu[i] = -1; }
        for (int j = 0; j < cnt; j++) {
            if(i*prime[j] >= N) break; vis[i*prime[j]] = 1;
            if(i%prime[j]==0) break; mu[i*prime[j]] = - mu[i];
        }
    } for (int i = 2; i < N; i++) mu[i] += mu[i-1];
}

ll dfs(ll n){ // 莫比乌斯前缀和
    if(n < N) return mu[n];
    if(mp.count(n)) return mp[n];
    ll ans = 1; // e
    for (ll i = 2,j; i <= n; i=j+1){
        j = n/(n/i); ans -= dfs(n/i)*(j-i+1); // 和上一题一样
    } return mp[n]=ans;
}

ll F(ll n) { // 分块求和
    ll ans = 0;
    for (ll i = 1,j; i <= n; i=j+1) {
        j = n/(n/i);
        ans = (ans + (j-i+1)%mod*(n/i)) % mod;
    } return ans * ans % mod;
}

int main() {
    init(); scanf("%lld", &n); ll ans = 0;
    for (ll i=1,j; i <= n; i=j+1) {
        j = n/(n/i);
        ans = (ans + (dfs(j) - dfs(i-1) + mod)%mod*F(n/i)%mod) % mod;
    } printf("%lld\n", ans);
}

Min_25筛

Min_25筛是一个替代洲阁晒的新产物，也是用于求解积性函数求和的问题，时间空间复杂度都比洲阁晒要优秀。

同样考虑筛质数，对1-n内所有质数求和。

设函数$S(x, j)=\sum_{i=2}^xi*[i为质数或i的最小质因子大于p_j]$ , 质数从p1开始

S(x, 0)=x*(x+1)/2 - 1; 我们要求的就是S(n, p0)

考虑转移：

$ p^2_j > x, S(x, j) = S(x, j-1) \否则 S(x,j)=S(x,j-1)-f(p_j)*(S(\lfloor \frac{x}{p_j} \rfloor,j-1)-\sum_{i=1}^{j-1}f(p_i)) $

求具体积性函数和时

$G(x,j)=G(x,j+1)+\sum_{k=1}(G(\lfloor\frac{x}{p_j^k}\rfloor,j+1)-\sum_{i=1}^xf(p_i) )*f(p_j^k)+\sum_{k=2}f(p_j^k)$

能力有限，以后再学。

本文标题:数论退役帖

文章作者:呼哧

发布时间:2018-08-20, 02:40:42

最后更新:2018-08-20, 02:32:47

原始链接:http://hu-chi.github.io/2018/08/20/acm/数论/

许可协议: "署名-非商用-相同方式共享 4.0" 转载请保留原文链接及作者。

文章目录
  1. 1. 数论退役帖
    1. 1.1. 基本知识
      1. 1.1.1. 最大公因数
        1. 1.1.1.1. 常用定理
      2. 1.1.2. 素数测定
      3. 1.1.3. 扩展欧几里得
      4. 1.1.4. 毕达哥拉斯三元组本原解
        1. 1.1.4.1. 欧拉函数与欧拉定理
      5. 1.1.5. 欧拉筛
      6. 1.1.6. 欧拉例题
        1. 1.1.6.1. 二次筛法 POJ - 2689
        2. 1.1.6.2. HDU - 3939 直角三角形个数
      7. 1.1.7. 逆元
      8. 1.1.8. 等比数列求和
        1. 1.1.8.1. POJ - 3233 等比矩阵求和
      9. 1.1.9. 组合数计算
      10. 1.1.10. 中国剩余定理(CRT)
      11. 1.1.11. 反素数
        1. 1.1.11.1. 求最小的数使其因数为n CF27E
        2. 1.1.11.2. 求1-n因数个数最大的数 URAL 1748
        3. 1.1.11.3. HDU - 4542
    2. 1.2. 指数循环节
    3. 1.3. 积性函数
      1. 1.3.0.1. 积性函数的定义
    4. 1.3.1. 常见的积性函数
      1. 1.3.1.1. 两个常用的公式
      2. 1.3.1.2. 积性函数性质
    5. 1.3.2. 分块技巧
    6. 1.3.3. 迪利克雷卷积
      1. 1.3.3.1. 积性函数例题 HDU - 5528
      2. 1.3.3.2. 积性函数习题 757E
    7. 1.3.4. 莫比乌斯反演
      1. 1.3.4.1. 经典例题-求互质的数对数 HDU - 1695
      2. 1.3.4.2. GCD表中的质数
    8. 1.3.5. 杜教筛
      1. 1.3.5.1. 欧拉函数前缀和 51nod 1239
      2. 1.3.5.2. 莫比乌斯函数前缀和
    9. 1.3.6. Extended Eratosthenes Sieve洲阁筛
      1. 1.3.6.1. 素数个数
    10. 1.3.7. 另外一些题
    11. 1.3.8. rng_58-clj等式
    12. 1.3.9. Min_25筛
|