目录
- 目录
- 素数的判定
- 素数的筛选构造素数表
- PollardRho大整数分解
- 幂乘快速计算二分思想
- 矩阵快速幂
- GCDLCM的相关
- 扩展欧几里德算法
- 威尔逊定理
- 欧拉定理
- 中国剩余定理
- 费马小定理
素数的判定
一个数 n 如果是合数,那么它的所有的因子不超过
sqrt(n)
,复杂度是
o(n*sqrt(n))
:
bool prime(int x){//判断x是不是质数,是返回true,不是返回false
if(x <= ) return false;
for(int i = ; i <= sqrt(x + ); i ++){//0.5是防止根号的精度误差
if(x % i == ) return false;
}
return true;
}
或者是这样写也行:
bool prime(int x){//判断x是不是质数,是返回true,不是返回false
if(x <= ) return false;
for(int i = ; i * i <= x; i ++){//用乘法避免根号的精度误差
if(x % i == ) return false;
}
return true;
}
//根据题目不同,如果i*i会爆int,记得开longlong
引入素数表来进行维护:
bool IsPrime(unsigned n)
{
if ( n < )
{ // 小于2的数即不是合数也不是素数
throw ;
}
static unsigned aPrimeList[] = { // 素数表
, , , , , , , , , , , , , ,
, , , , , , , , , , , , ,
, , , , , , , , , , ,
, , , , , , , , , ,
,, , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, , , , , , , , ,
, ,
};
const int nListNum = sizeof(aPrimeList)/sizeof(unsigned);//计算素数表里元素的个数
for (unsigned i=;i<nListNum;++i )
{
if(n/+<aPrimeList[i])
{
return true;
}
if(==n%aPrimeList[i])
{
return false;
}
}
for (unsigned i=aPrimeList[nListNum-];i<n/+;i++ )
{
if (==n%i)
{ // 除尽了,合数
return false;
}
}
return true;
}
二次探测定理:如果p是素数,且
0<x<p
,则方程
x^2=1(mod p)
的解为1或p-1。
最后一种情况,面对大数时采用Miller-Rabbin素数测试法:
Miller-Rabbin素数测试是一个不确定的算法,只能从概率意义上判定一个数可能是素数,但并不能确保。算法流程如下:
1.选择T个随机数A,并且有A
int Miller_Rabbin(long long n)
{
long long i,s,a;
s=; //s的值可以根据需要变大
// randomize();
for(i=;i <s;i++)
{
a=long long(rand()%(n-)+); //自动生成受限
if(modular_exp(a,n-,n)>)
return ;
}
return ;
}
long long modular_exp(long long a,long long b,long long c)//求a^b%c该函数受限
{
if(a==)
return ;
if(b==)
return ;
if(b==)
return a%c;
return (a*modular_exp(a,b-,c))%c;
}
完整版代码:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
using namespace std;
const int Times = ;
typedef long long LL;
LL multi(LL a, LL b, LL m)
{
LL ans = ;
a %= m;
while(b)
{
if(b & )
{
ans = (ans + a) % m;
b--;
}
b >>= ;
a = (a + a) % m;
}
return ans;
}
LL quick_mod(LL a, LL b, LL m)
{
LL ans = ;
a %= m;
while(b)
{
if(b & )
{
ans = multi(ans, a, m);
b--;
}
b >>= ;
a = multi(a, a, m);
}
return ans;
}
bool Miller_Rabin(LL n)
{
if(n == ) return true;
if(n < || !(n & )) return false;
LL m = n - ;
int k = ;
while((m & ) == )
{
k++;
m >>= ;
}
for(int i=; i<Times; i++)
{
LL a = rand() % (n - ) + ;
LL x = quick_mod(a, m, n);
LL y = ;
for(int j=; j<k; j++)
{
y = multi(x, x, n);
if(y == && x != && x != n - ) return false;
x = y;
}
if(y != ) return false;
}
return true;
}
int main()
{
int T;
scanf("%d",&T);
while(T--)
{
LL n;
scanf("%I64d",&n);
if(Miller_Rabin(n)) puts("Yes");
else puts("No");
}
return ;
}
质数检测Java大数版:
import java.io.*;
import java.util.*;
import java.math.BigInteger;
public class Main{
public static final int Times = ;
public static BigInteger quick_mod(BigInteger a,BigInteger b,BigInteger m){
BigInteger ans = BigInteger.ONE;
a = a.mod(m);
while(!(b.equals(BigInteger.ZERO))){
if((b.mod(BigInteger.valueOf())).equals(BigInteger.ONE)){
ans = (ans.multiply(a)).mod(m);
b = b.subtract(BigInteger.ONE);
}
b = b.divide(BigInteger.valueOf());
a = (a.multiply(a)).mod(m);
}
return ans;
}
public static boolean Miller_Rabin(BigInteger n){
if(n.equals(BigInteger.valueOf())) return true;
if(n.equals(BigInteger.ONE)) return false;
if((n.mod(BigInteger.valueOf())).equals(BigInteger.ZERO)) return false;
BigInteger m = n.subtract(BigInteger.ONE);
BigInteger y = BigInteger.ZERO;
int k = ;
while((m.mod(BigInteger.valueOf())).equals(BigInteger.ZERO)){
k++;
m = m.divide(BigInteger.valueOf());
}
Random d = new Random();
for(int i=;i<Times;i++){
int t = ;
if(n.compareTo(BigInteger.valueOf()) == ){
t = ;
}else{
t = n.intValue() - ;
}
int a = d.nextInt(t) + ;
BigInteger x = quick_mod(BigInteger.valueOf(a),m,n);
for(int j=;j<k;j++){
y = (x.multiply(x)).mod(n);
if(y.equals(BigInteger.ONE) && !(x.equals(BigInteger.ONE)) && !(x.equals(n.subtract(BigInteger.ONE)))) return false;
x = y;
}
if(!(y.equals(BigInteger.ONE))) return false;
}
return true;
}
public static void main(String[] args){
Scanner cin = new Scanner(System.in);
while(cin.hasNextBigInteger()){
BigInteger n = cin.nextBigInteger();
if(Miller_Rabin(n)) System.out.println("Yes");
else System.out.println("No");
}
}
}
素数的筛选(构造素数表)
不采用筛素:假设我们已经我素数序列: p1, p2, .. pn
现在要判断pn+1是否是素数, 则需要(1, sqrt(pn+1)]范围内的所有素数序列,
而这个素数序列显然已经作为p1, p2, .. pn的一个子集被包含了
void makePrimes(int primes[], int num)
{
int i, j, cnt;
primes[] = ;
primes[] = ;
for(i = , cnt = ; cnt < num; i += )
{
int flag = true;
for(j = ; primes[j]*primes[j] <= i; ++j)
{
if(i%primes[j] == )
{
flag = false; break;
}
}
if(flag) primes[cnt++] = i;
}
}
可以使用筛素法,复杂度
O(nloglogn)
:
bool prime[N];
void init(){
for(int i = ; i < N; i ++) prime[i] = true;
for(int i = ; i*i < N; i ++){//判断改成i*i<N
if(prime[i]){
for(int j = i*i; j < N; j += i){//从i*i开始就可以了
prime[j] = false;
}
}
}
}
或者采用线性筛素法
O(n)
:
#include<cstdio>
const int N = + ;
bool prime[N];//prime[i]表示i是不是质数
int p[N], tot;//p[N]用来存质数
void init(){
for(int i = ; i < N; i ++) prime[i] = true;//初始化为质数
for(int i = ; i < N; i++){
if(prime[i]) p[tot ++] = i;//把质数存起来
for(int j = ; j < tot && i * p[j] < N; j++){
prime[i * p[j]] = false;
if(i % p[j] == ) break;//保证每个合数被它最小的质因数筛去
}
}
}
int main(){
init();
}
Meisell-Lehmer 模板,求1-n范围内的素数个数:
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = + ;
bool np[N];
int prime[N], pi[N];
int getprime() {
int cnt = ;
np[] = np[] = true;
pi[] = pi[] = ;
for(int i = ; i < N; ++i) {
if(!np[i]) prime[++cnt] = i;
pi[i] = cnt;
for(int j = ; j <= cnt && i * prime[j] < N; ++j) {
np[i * prime[j]] = true;
if(i % prime[j] == ) break;
}
}
return cnt;
}
const int M = ;
const int PM = * * * * * * ;
int phi[PM + ][M + ], sz[M + ];
void init() {
getprime();
sz[] = ;
for(int i = ; i <= PM; ++i) phi[i][] = i;
for(int i = ; i <= M; ++i) {
sz[i] = prime[i] * sz[i - ];
for(int j = ; j <= PM; ++j) {
phi[j][i] = phi[j][i - ] - phi[j / prime[i]][i - ];
}
}
}
int sqrt2(LL x) {
LL r = (LL)sqrt(x - );
while(r * r <= x) ++r;
return int(r - );
}
int sqrt3(LL x) {
LL r = (LL)cbrt(x - );
while(r * r * r <= x) ++r;
return int(r - );
}
LL getphi(LL x, int s) {
if(s == ) return x;
if(s <= M) return phi[x % sz[s]][s] + (x / sz[s]) * phi[sz[s]][s];
if(x <= prime[s]*prime[s]) return pi[x] - s + ;
if(x <= prime[s]*prime[s]*prime[s] && x < N) {
int s2x = pi[sqrt2(x)];
LL ans = pi[x] - (s2x + s - ) * (s2x - s + ) / ;
for(int i = s + ; i <= s2x; ++i) {
ans += pi[x / prime[i]];
}
return ans;
}
return getphi(x, s - ) - getphi(x / prime[s], s - );
}
LL getpi(LL x) {
if(x < N) return pi[x];
LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - ;
for(int i = pi[sqrt3(x)] + , ed = pi[sqrt2(x)]; i <= ed; ++i) {
ans -= getpi(x / prime[i]) - i + ;
}
return ans;
}
LL lehmer_pi(LL x) {
if(x < N) return pi[x];
int a = (int)lehmer_pi(sqrt2(sqrt2(x)));
int b = (int)lehmer_pi(sqrt2(x));
int c = (int)lehmer_pi(sqrt3(x));
LL sum = getphi(x, a) + LL(b + a - ) * (b - a + ) / ;
for (int i = a + ; i <= b; i++) {
LL w = x / prime[i];
sum -= lehmer_pi(w);
if (i > c) continue;
LL lim = lehmer_pi(sqrt2(w));
for (int j = i; j <= lim; j++) {
sum -= lehmer_pi(w / prime[j]) - (j - );
}
}
return sum;
}
int main() {
init();
LL n;
while(cin >> n) {
cout << lehmer_pi(n) << endl;
}
return ;
}
PollardRho大整数分解
Pollard rho算法的原理就是通过某种方法得到两个整数a和b,而待分解的大整数为n,计算p=gcd(a-b,n),直到p不为1,或者a,b出现循环为止。然后再判断p是否为n,如果p=n成立,那么返回n是一个质数,否则返回p是n的一个因子,那么我们又可以递归的计算Pollard(p)和Pollard(n/p),这样,我们就可以求出n的所有质因子。
模板如下:
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <stdio.h>
const int Times = ;
const int N = ;
using namespace std;
typedef long long LL;
LL ct, cnt;
LL fac[N], num[N];
LL gcd(LL a, LL b)
{
return b? gcd(b, a % b) : a;
}
LL multi(LL a, LL b, LL m)
{
LL ans = ;
a %= m;
while(b)
{
if(b & )
{
ans = (ans + a) % m;
b--;
}
b >>= ;
a = (a + a) % m;
}
return ans;
}
LL quick_mod(LL a, LL b, LL m)
{
LL ans = ;
a %= m;
while(b)
{
if(b & )
{
ans = multi(ans, a, m);
b--;
}
b >>= ;
a = multi(a, a, m);
}
return ans;
}
bool Miller_Rabin(LL n)
{
if(n == ) return true;
if(n < || !(n & )) return false;
LL m = n - ;
int k = ;
while((m & ) == )
{
k++;
m >>= ;
}
for(int i=; i<Times; i++)
{
LL a = rand() % (n - ) + ;
LL x = quick_mod(a, m, n);
LL y = ;
for(int j=; j<k; j++)
{
y = multi(x, x, n);
if(y == && x != && x != n - ) return false;
x = y;
}
if(y != ) return false;
}
return true;
}
LL pollard_rho(LL n, LL c)
{
LL i = , k = ;
LL x = rand() % (n - ) + ;
LL y = x;
while(true)
{
i++;
x = (multi(x, x, n) + c) % n;
LL d = gcd((y - x + n) % n, n);
if( < d && d < n) return d;
if(y == x) return n;
if(i == k)
{
y = x;
k <<= ;
}
}
}
void find(LL n, int c)
{
if(n == ) return;
if(Miller_Rabin(n))
{
fac[ct++] = n;
return ;
}
LL p = n;
LL k = c;
while(p >= n) p = pollard_rho(p, c--);
find(p, k);
find(n / p, k);
}
int main()
{
LL n;
while(cin>>n)
{
ct = ;
find(n, );
sort(fac, fac + ct);
num[] = ;
int k = ;
for(int i=; i<ct; i++)
{
if(fac[i] == fac[i-])
++num[k-];
else
{
num[k] = ;
fac[k++] = fac[i];
}
}
cnt = k;
for(int i=; i<cnt; i++)
cout<<fac[i]<<"^"<<num[i]<<" ";
cout<<endl;
}
return ;
}
幂乘快速计算(二分思想)
逆向思维,不断地二分,累积底数的平方即可:
typedef long long LL;
LL pow_mod(LL a, LL b, LL p){//a的b次方求余p
LL ret = ;
while(b){
if(b & ) ret = (ret * a) % p;
a = (a * a) % p;
b >>= ;
}
return ret;
}
衍生出ll的快速乘:
LL mul(LL a, LL b, LL p){//快速乘,计算a*b%p
LL ret = ;
while(b){
if(b & ) ret = (ret + a) % p;
a = (a + a) % p;
b >>= ;
}
return ret;
}
关于计算结果过除以MOD取余的问题:
1、计算加法每相加一次执行一次取余;
2、计算减法给被减数加上MOD之后先减法后取余;
3、计算乘法每相乘一次取余一次。
矩阵快速幂
算一个
n*n
矩阵的的x次方对1e9+7求余:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = ;
const int MOD = +;
#define mod(x) ((x)%MOD)
int n;
struct mat{
int m[maxn][maxn];
}unit;
mat operator * (mat a,mat b){
mat ret;
ll x;
for(int i = ;i < n;i++){
for(int j = ;j < n;j++){
x = ;
for(int k = ;k < n;k++)
x += mod((ll)a.m[i][k]*b.m[k][j]);
ret.m[i][j] = mod(x);
}
}
return ret;
}
void init_unit(){
for(int i = ;i < maxn;i++)
unit.m[i][i] = ;
return ;
}
mat pow_mat(mat a,ll n){
mat ret = unit;
while(n){
if(n&) ret = ret*a;
a = a*a;
n >>= ;
}
return ret;
}
int main(){
ll x;
init_unit();
while(cin>>n>>x){
mat a;
for(int i = ;i < n;i++)
for(int j = ;j < n;j++)
cin>>a.m[i][j];
a = pow_mat(a,x);
for(int i = ;i < n;i++){
for(int j = ;j < n;j++){
if(j + == n) cout << a.m[i][j] << endl;
else cout<<a.m[i][j]<<" ";
}
}
}
return ;
}
GCD&LCM的相关
计算gcd循环写法和递归写法:
LL gcd(LL a, LL b){
if(b == ) return a;
else return gcd(b, a%b);
}
LL gcd(LL a, LL b){
return b ? gcd(b, a%b) : a;
}
LL gcd(LL a, LL b){
LL t;
while(b){
t = b;
b = a % b;
a = t;
}
return a;
}
相关的公式:
a*b = gcd(a,b) * lcm(a,b)
gcd(ka, kb) = k * gcd(a, b)
lcm(ka, kb) = k * lcm(a, b)
lcm(S/a, S/b) = S/gcd(a, b)
扩展欧几里德算法
那么已知 a,b 求 一组解 x,y 满足 ax+by = gcd(a, b) 这个公式:
代码模板:
#include<cstdio>
typedef long long LL;
void ex_Eulid(LL a, LL b, LL &x, LL &y, LL &d){
if (!b) {d = a, x = , y = ;}
else{
extend_Eulid(b, a % b, y, x, d);
y -= x * (a / b);
}
}
int main(){
LL a, b, d, x, y;
while(~scanf("%lld%lld", &a, &b)){
ex_Eulid(a, b, x, y, d);
printf("%lld*a + %lld*b = %lld\n", x, y, d);
}
}
威尔逊定理
当且仅当p为素数时:( p -1 )! ≡ -1 ( mod p )
或者这么写( p -1 )! ≡ p-1 ( mod p )
若p为质数,则p能被(p-1)!+1整除
在初等数论中,这是威尔逊给出了判定一个自然数是否为素数的充分必要条件。
欧拉定理
欧拉定理,也称费马-欧拉定理
若n,a为正整数,且n,a互质,即gcd(a,n) = 1,则a^φ(n) ≡ 1 (mod n)
φ(n) 是欧拉函数,欧拉函数是求小于等于n的数中与n互质的数的数目
欧拉函数是求 (小于n的数 )中 (与n互质的数 )的数目或者说欧拉函数是求 1到n-1 中 与n互质的数的数目。
如果n是质数那么1到n-1所有数都是与n互质的,所以φ(n) = n-1
欧拉函数模板:
#include<cstdio>
const int N = + ;
int phi[N];
void Euler(){
phi[] = ;
for(int i = ; i < N; i ++){
if(!phi[i]){
for(int j = i; j < N; j += i){
if(!phi[j]) phi[j] = j;
phi[j] = phi[j] / i * (i-);
}
}
}
}
int main(){
Euler();
}
另一种,比上面更快的方法需要用到如下性质
p为质数时有:
phi(p)=p-1 因为质数p除了1以外的因数只有p,故1至p的整数只有p与p不互质 如果i mod p = 0, 那么 phi(i * p)=phi(i) * p
若i mod p ≠0, 那么 phi( i * p )=phi(i) * ( p-1 )
#include<cstdio>
using namespace std;
const int N = + ;
int phi[N], prime[N];
int tot;//tot计数,表示prime[N]中有多少质数
void Euler(){
phi[] = ;
for(int i = ; i < N; i ++){
if(!phi[i]){
phi[i] = i-;
prime[tot ++] = i;
}
for(int j = ; j < tot && l*i*prime[j] < N; j ++){
if(i % prime[j]) phi[i * prime[j]] = phi[i] * (prime[j]-);
else{
phi[i * prime[j] ] = phi[i] * prime[j];
break;
}
}
}
}
int main(){
Euler();
}
中国剩余定理
CRT问题代码描述如下:
void extend_Euclid(LL a, LL b, LL &x, LL &y)
{
if(b == )
{
x = ;
y = ;
return;
}
extend_Euclid(b, a % b, x, y);
LL tmp = x;
x = y;
y = tmp - (a / b) * y;
}
int CRT(int a[],int m[],int n)
{
int M = ;
int ans = ;
for(int i=; i<=n; i++)
M *= m[i];
for(int i=; i<=n; i++)
{
int x, y;
int Mi = M / m[i];
extend_Euclid(Mi, m[i], x, y);
ans = (ans + Mi * x * a[i]) % M;
}
if(ans < ) ans += M;
return ans;
}
费马小定理
假如p是质数:
若p不能整除a,则 a^(p-1) ≡1(mod p);
若p能整除a,则a^(p-1) ≡0(mod p)。
或者说,若p是质数,且a,p互质,那么 a的(p-1)次方除以p的余数恒等于1。