math--primes

数论 - 质数（素数）与约数知识

素数相关

概念

对于一个正整数$n$，除$1$与$n$自身外，不能被其他数整除的数就叫做质数（素数），质数在数论学习里面是一个十分基础重要的知识，务必掌握

补充:

typedef long long ll;

素数判断

朴素版

此版本可以用于 $int$ 范围内的素数判断

bool isprime(ll n) {
    for(ll i = 2; i < n / i; i++) {
        if(n % i == 0) return false;
    }
    return true;
}

优化版

此版本可以优化循环次数，$long$范围内可用

bool isprime(ll n) {
    if (n <= 1) return false; 
    if (n <= 3) return true;  // 对比较小的数进行特判
    //对于比较容易判别的数（被2， 3 整除的数）先判断
    if (n % 2 == 0 || n % 3 == 0) return false; 
    //之后可以优化循环
    for (ll i = 5; i <= n / i; i += 6) 
      	if (n % i == 0 || n % (i+2) == 0) 
         	return false; 
  	return true; 
}

进阶版—素数性测试

此版本可用于$10^{18}$以下的素数判断，采用费马小定理以及快速幂进行优化，将判别时间压到常数级别，极大程度优化素数判定

对于费马小定理，当被检测数$n$是素数时，必有$a^{n - 1}\equiv 1(mod \ n)$

并且对于素数(除二外，并且小的素数可以通过特判快速排除) $n$ ，$n$ 必是奇数，故可以将$n - 1$分解：

​ $n - 1 \equiv 2^s \times d$ $(d \equiv odd)$

故而可以以此分解费马小定理的方程：

​ $a^{n-1} \equiv 1\ (mod \ n)\ \leftrightarrow a^{2^s\times d}\ - 1\equiv 0\ (mod \ n)$

​ $ \leftrightarrow (a^{2^{s-1} \times d}\ +1)(a^{2^{s-1} \times d\ } - 1) \equiv 0\ (mod\ n)$

​ $ \leftrightarrow (a^{2^{s-1} \times d}\ +1)(a^{2^{s-2} \times d} + 1)(a^{2^{s-2} \times d} - 1) \equiv 0\ (mod\ n)$

​ …….

​ $ \leftrightarrow (a^{2^{s-1} \times d}\ +1)(a^{2^{s-2} \times d} + 1)\ …\ (a^d+1)(a^d-1) \equiv 0\ (mod\ n)$

由分解可知，如果 $n$ 是一个质数，那么必符合上诉公式，则必有：

​ $a^d \equiv 1\ (mod\ n)$

或：

​ $a^{2^r\times d} \equiv -1\ (mod\ n)$ $(0\leq r\leq s-1)$

综上，我们可以选取一个随机数 $a$ $(a \leq n-1)$ ，通过在模 $n$ 情况下进行快速幂，最后检验两种结果的方法判断 $n$ 是否是素数，当然，有一些非素数对于小于也是符合上诉公式的，故需要进行多次素数检验（一般20次）才可判断 $n$ 是素数，并且因为是模以大素数下的计算，快速幂会破 $long\ long$ 故需要开 $uint64_t$ 以及 $__unit128_t $ 作为中间值

using u64 = uint64_t;
using u128 = __uint128_t;
// typedef unit64_t ud4;
// typedef __uint128_t u128;

//使用快速幂求出次方
u64 qmi(u64 a, u64 k, u64 mod) {
    // 计算a的k次方在模mod意义下的值
    u64 ans = 1;
    a %= mod;	//防止溢出，控制a小于mod
    while(k) {
        if(k & 1) ans = (u128)ans * a % mod;
        a = (u128)a * a % mod;
        k >>= 1;
    }
    return ans;
}

//对于某一个a,检查n是不是显示出质数的性质(是否符合公式推导),不符合返回true
bool check_prime(u64 n, u64 a, u64 d, int s) {
    u64 x = qmi(a, d, n);
    if(x == 1 || x == n - 1) return false;
    for(int r = 1; r < s; r++) {
        x = (u128)x * x % n;
        if(x == n - 1)
            return false;
    }
    return true;
}

//判断n是不是素数，是返回true,不是返回false,测试次数iter一般在10以内即可
bool MillerRabin(u64 n, int iter = 10) {
    if(n < 4)
        return n == 2 || n == 3;
    if(n % 2 == 0 || n % 3 == 0) return false;
    int s = 0;
    u64 d = n - 1;
    //获得s还有d
    while((d & 1) == 0) {
        d >>= 1;
        s++;
    }
    
    for(int i = 0; i < iter; i++) {
        u64 a = 2 + rand() % (n - 3);
        if(check_prime(n, a, d, s))
            return false;
    }
    return true;
}

当然，因为概率性原因，米勒~罗宾素数判定法会存在误差，但通过概率公式推导可以发现，误差概率极小，故可以忽略不计，并且根据素数的性质判断，我们可以先选取排名前几的素数进行素数推理，可以将误差进一步压低为近似为0，可以放心使用

确定性版本为：

bool MillerRabin(ud4 n) {
     if(n < 4)
        return n == 2 || n == 3;
    if(n % 2 == 0 || n % 3 == 0) return false;
    
    int s = 0;
    u64 d = n - 1;
    //获得s还有d
    while((d & 1) == 0) {
        d >>= 1;
        s++;
    }
    
    //以下基于c++11以上写法，是可以修改成c的
    for(u64 a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        if(n == a)
            return true;
        if(check_prime(n, a, d, r))
            return false;
    }
    return true;
}

素数筛法

补充：在下列筛法中，存在数组元素的定义

const int N = 1e6 + 10;	 //自定义设置，一般不超过1e6
bool st[N]; 			//真值表，否为素数，真为合数、
vector<ll> primes;		//素数表

朴素筛法（时间复杂度过高，不建议用）

朴素法为对于$1-n$的所有数，进行素数判断，如果是素数就加入素数表

时间复杂度$n\sqrt{n}$

vector<ll> get_primes(ll n) { //此时n为需要筛的1--n的区间大小
    vector<ll> primes;
   	for(ll i = 2; i <= n; i++) {
        if(isprime(i)) primes.push_back(i);
    }
    return primes;
}

埃式筛法

埃式筛法时间复杂度为$nlog_2{log_2n}$ 速度较朴素法大幅提高

该算法通过用质因子对其合数进行筛选，故筛选时需要打真值表(注：真值表长度需要大于n)判别

vector<ll> get_primes(ll n) {
    vector<ll> primes;
    for(ll i = 2; i <= n; i++) {
        if(!st[i]) {
            primes.push_back(i);
            for(ll j = i * 2; j <= n; j += i) 
                st[j] = true;
        }
    }
    return primes;
}

线性筛法

此法通过数学推导可证，每一个合数只会被其最小质因数筛除一次，故时间复杂度是$n$,属于极佳筛法

vector<ll> get_primes(ll n) {
    vector<ll> primes;
    for(ll i = 2; i <= n; i++) {
        if(!st[i]) primes.push_back(i);
        for(int j = 0; primes[j] * i <= n; j++) {
            st[primes[j] * i] = true;
            if(i % primes[j] == 0) break;
        }
    }
    return primes;
}

区间筛法

这个筛法比较特殊，一般使用于筛大数的一段区间的质数，区间在$1e6-1e7$之间, 最大值不超过int范围，采用埃筛性质，时间复杂度与埃筛接近

vector<ll> pres = get_primes(1000010);
bool st[1000010]; //真值表大小对应区间大小
vector<ll> get_section_primes(ll l, ll r) {
    memset(st, 0, sizeof st);
    for(auto p : pres) {
        for(ll i = max((l - 1 + p) / p * p, p * 2); i <= r; i += p) {
            st[i - l] = 1;
        }
    }
    vector<ll>primes;
    for(int i = 0; i < r - l + 1; i++) {
        if(!st[i] && i + l > 1) primes.push_back(i + l);
    }
    return primes;
}

约数相关

概念

对于任意一个数字$n$，都拥有约数，即小于等于$n$并且可$n$其整除的数，特别的，质数的约数个数为2个

约数个数

补充：

const int N = 1600 + 10; //因为在1e9里面的数，约数个数不会超过1600,(暴力求解过)
vector<int>divs;		//存每个数的约数个数的数组

朴素筛法

对于任意一个元素$n$,朴素筛法只会循环到最多$\sqrt{n}$ 过，故时间复杂度为$\sqrt{n}$ ,一般适用

此方法可以求出$n$的所有约数还有约数个数

vector<ll> get_divs(ll n) {
    vector<ll> divs;
    for(ll i = 1; i <= n / i; i++) {
        if(n % i == 0) {
            divs.push_back(i);
            if(n / i != i) divs.push_back(n / i);
        }
    }
    return divs;
}

约数进阶求法1

根据题目不同，可以对代码进行不同程度优化，以下为例

求一个数的所有质因数

vector<ll> get_divs(ll n) {
    if(isprime(n)) return {n};
    vector<ll> divs;
    for(ll i = 2; i <= n / i; i++) {
        if(n % i == 0) {
            divs.push_back(i);
            while(n % i == 0) n /= i;
        }
    }
    if(n > 1) divs.push_back(n);
    return divs;
}

求一个数的所有质因数以及其个数

注：$p$为质数，$s$为质数的个数

typedef struct Factor{
    ll p;
    ll s;
}factor;

法1：单此查询时可用，速度较快，效率较高

vector<factor> get_factor(ll n) {
    vector<factor> facs;
    for(ll i = 2; i <= n / i; i++) {
        if(n % i == 0) {
            int cnt = 0;
            while(n % i == 0) cnt++, n /= i;
        }
        facs.push_back({i, cnt});
    }
    if(n > 1) facs.push_back({n, 1});
    return facs;
}

法2：多次查询时用，需要额外打素数表，因为打了素数表的原因，故多次查询下比法一快多

vector<ll> primes = get_primes(1000010);	// 此处的数字要比多次查询下的最大根号n略大
vector<factor> get_factor(ll n) {
    vector<factor> facs;
    for(int i = 0; primes[i] <= n / primes[i]; i++) {
         if(n % primes[i] == 0) {
            int cnt = 0;
            while(n % primes[i] == 0) cnt++, n /= primes[i];
        }
        facs.push_back({primes[i], cnt});
    }
    if(n > 1) facs.push_back({n, 1});
    return facs;
}

利用以上信息快速求解约数

基于已经拥有了一个数的所有质因子以及个数，故可用$dfs$排列组合出所有的因子

void dfs(vector<factor> &facs, vector<ll> &divs, ll u, ll p) {
    if(u == facs.size()) {
        divs.push_back(p);
        return;
    }
    for(int i = 0; i <= facs[u].s; i++) {
        dfs(facs, divs, u + 1, p);
        p *= facs[u].p;
    }
}


vector<ll> get_divs(ll n) {
    vector<factor> facs = get_factor(n);
    vector<ll> divs;
    dfs(facs, divs, 0, 1);
    return divs;
}

例题：多次询问不同的$n$ 的约数

地址：https://www.acwing.com/problem/content/description/202/

#include <vector>
#include <iostream>
#include <algorithm>

using namespace std;
typedef long long ll;

typedef struct Factor{
    ll p;
    ll s;
}factor;

bool st[100020];

vector<ll> get_primes(ll n) {
    vector<ll> primes;
    for(ll i = 2; i <= n; i++) {
        if(!st[i]) primes.push_back(i);
        for(int j = 0; primes[j] * i <= n; j++) {
            st[primes[j] * i] = true;
            if(i % primes[j] == 0) break;
        }
    }
    return primes;
}

vector<ll> primes;

vector<factor> get_factor(ll n) {
    vector<factor> facs;
    for(int i = 0; primes[i] <= n / primes[i]; i++) {
         if(n % primes[i] == 0) {
            ll cnt = 0;
            while(n % primes[i] == 0) cnt++, n /= primes[i];
            facs.push_back({primes[i], cnt});
        }
    }
    if(n > 1) facs.push_back({n, 1});
    return facs;
}

void dfs(vector<factor> &facs, vector<ll> &divs, ll u, ll p) {
    if(u == facs.size()) {
        divs.push_back(p);
        return;
    }
    for(int i = 0; i <= facs[u].s; i++) {
        dfs(facs, divs, u + 1, p);
        p *= facs[u].p;
    }
}


vector<ll> get_divs(ll n) {
    vector<factor> facs = get_factor(n);
    vector<ll> divs;
    dfs(facs, divs, 0, 1);
    return divs;
}

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

int main() {
    primes = get_primes(100010);
    int T;
    cin >> T;
    while(T --) {
        ll a, b, c, d;
        cin >> a >> b >> c >> d;
        vector<ll> divs = get_divs(d);
        ll res = 0;
        for(auto x : divs) {
            if(gcd(a, x) == b && c * x / gcd(c, x) == d) res++;
        }
        cout << res << endl;
    }
    return 0;
}