关于数论学习的一点点笔记以及手写的板子代码（一直更新）
//中间有些图片不想导入了，反正本地有：P
快速幂
int quickpow(int base, int power) {
int res = 1;
while (power > 0) {
if (power & 1) res = res * base % mod;
power >>= 1;
base = (base *base) % mod;
}
return res;
}
慢速乘
int mul(int a, int b, int mod)
{
int res = 0;
while (b != 0)
{
if (b & 1) res = (res + a) % mod;
b /= 2;
a = (a + a) % mod;
}
return res % mod;
}
取对数
int quickpow(int base, int power) {
int res = 1;
while (power > 0) {
if (power & 1) res = res * base % mod;
power >>= 1;
base = (base *base) % mod;
}
return res;
}
int log2(int x) {
for (int i = 0; i <= 64; i++) {
if (quickpow(2, i) >= x)
return i;
}
}
逆元
//费马小定理 & 欧拉定理
ll quickpow(ll a,ll b){
ll res = 1;
while(b > 0){
if(b & 1)res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}
ll getinv(int a,int mod)//求a在mod下的逆元，不存在逆元返回-1
{
return quickpow(a,mod - 2);
}
//a的逆元 inv[a] = quickpow(a,mod - 2);
//逆元线性筛 O(n)：
const int mod = 1000000009;
const int maxn = 10005;
int inv[maxn];
inv[1] = 1;
for(int i = 2; i < 10000; i++)
inv[i] = inv[mod % i] * (mod - mod / i) % mod;
//阶乘的逆元 O(logn+n)：
inv[maxn]=mod_pow(fac[maxn],mod-2);
for(ll i=maxn-1;i>=0;i--)
inv[i]=(inv[i+1]*(i+1))%mod;
位运算大数相乘（1e18)
//0 < a,b,p < 1e18 ;
//求a * b % p
//原理把乘法变成加法
ll quick_add(ll a,ll b,ll p)
{
ll res=0;
while(b)
{
if(b&1) res=(res+a)%p;
a=(a+a)%p;
b>>=1;
}
return res;
}
杨辉三角
#include <iostream>
#include <iomanip>
using namespace std;
/*每个数等于它上方两数之和*/
int main(){
const int n=11;
int i,j,a[n][n];
for (i=1;i<n;i++){//前两行全为1，拿出来单独处理
a[i][i]=1;//使最右侧边全为1
a[i][1]=1;//使最左侧边全为1
}//从第三行开始处理
for (i=3;i<n;i++)
for (j=2;j<=i-1;j++) //j始终慢i一步
a[i][j]=a[i-1][j-1]+a[i-1][j];
for (i=1;i<n;i++){
for (j=1;j<=i;j++)
cout<<setw(5)<<a[i][j]<<" ";
cout<<endl;
}
cout<<endl;
return 0;
}
卡特兰数
卡特兰数的应用都可以归结到一种情况：有两种操作，分别为操作一和操作二，它们的操作次数相同，都为 N，且在进行第 K 次操作二前必须先进行至少 K 次操作一，问有多少中情况？结果就Catalan(N)。
通项公式：
Catalan数的典型应用：
1.由n个+1和n个-1组成的排列中，满足前缀和>=0的排列有Catalan(N)种。
2.括号化问题。左括号和右括号各有n个时，合法的括号表达式的个数有Catalan(N)种。
3.有n+1个数连乘，乘法顺序有Catalan(N)种,相当于在式子上加括号。
4.n个数按照特定顺序入栈，出栈顺序随意，可以形成的排列的种类有Catalan(N)种。
5.给定N个节点，能构成Catalan(N)种种形状不同的二叉树。
6.n个非叶节点的满二叉树的形态数为Catalan(N)。
7.对于一个n*n的正方形网格，每次只能向右或者向上移动一格，不能穿越对角线，那么从左下角到右上角的不同种类有Catalan(N)种。
8.对于在n位的2进制中，有m个0，其余为1的catalan数为：C（n,m）-C(n,m-1)。
9.对凸n+2边形进行不同的三角形分割（只连接顶点对形成n个三角形）数为Catalan(N)。
10.将有2n个元素的集合中的元素两两分为n个子集，若任意两个子集都不交叉，那么我们称此划分为一个不交叉划分。此时不交叉的划分数是Catalan(N)。
11.n层的阶梯切割为n个矩形的切法数也是Catalan(N)。
12.在一个2*n的格子中填入1到2n这些数值使得每个格子内的数值都比其右边和上边的所有数值都小的情况数也是Catalan(N)。
// 1、前三十项卡特兰数表
[1,1,2,5,14,42,132,429,1430,4862,16796,58786,
208012,742900,2674440,9694845,35357670,129644790,
477638700,1767263190,6564120420,24466267020,
91482563640,343059613650,1289904147324,
4861946401452,18367353072152,69533550916004,
263747951750360,1002242216651368,3814986502092304]
//2、卡特兰数求模模板
const int C_maxn = 1e4 + 10;
LL CatalanNum[C_maxn];
LL inv[C_maxn];
inline void Catalan_Mod(int N, LL mod)
{
inv[1] = 1;
for(int i=2; i<=N+1; i++)///线性预处理 1 ~ N 关于 mod 的逆元
inv[i] = (mod - mod / i) * inv[mod % i] % mod;
CatalanNum[0] = CatalanNum[1] = 1;
for(int i=2; i<=N; i++)
CatalanNum[i] = CatalanNum[i-1] * (4 * i - 2) %mod * inv[i+1] %mod;
}
//3.求n<=35以内的卡特兰数
ll h[36];
void init()
{
int i,j;
h[0]=h[1]=1;
for(i=2;i<36;i++)
{
h[i]=0;
for(j=0;j<i;j++)
h[i]=h[i]+h[j]*h[i-j-1];
}
}
printf("%lld\n",h[n]);
// 4.快速求第n位卡特兰数模板（mod1e9+7版）
const long long M=1000000007;
long long inv[1000010];
long long last,now=1;
void init()
{
inv[1]=1;
for(int i=2;i<=n+1;i++)inv[i]=(M-M/i)*inv[M%i]%M;
}
int main()
{
scanf("%lld",&n);
init();
for(int i=2;i<=n;i++)
{
last=now;
now=last*(4*i-2)%M*inv[i+1]%M;
}
printf("%lld\n",last);
return 0;
}
//5.Java大数打表卡特兰数模板
import java.io.*;
import java.math.BigInteger;
import java.util.*;
public class Main {
public static void main(String[] args) {
Scanner cin=new Scanner(System.in);
BigInteger s[]=new BigInteger[105];
s[1]=BigInteger.ONE;
for(int i=2;i<105;i++){
s[i]=s[i-1].multiply(BigInteger.valueOf((4*i-2))).divide(BigInteger.valueOf(i+1));
}
while(cin.hasNext()){
int n=cin.nextInt();
System.out.println(s[n]);
}
}
}
//6.通用卡特兰数打表模板：
int a[105][250];
void ktl()
{
int i,j,yu,len;
a[1][0]=1;
a[1][1]=1;
len=1;
for(i=2;i<101;i++)
{
yu=0;
for(j=1;j<=len;j++)
{
int t=(a[i-1][j])*(4*i-2)+yu;
//如果是求考虑顺序的排列，如不同点敏感的n节点的二叉树种类则改成这句即可
//int t=(a[i-1][j])*(4*i-2)*i+yu;
yu=t/10;
a[i][j]=t%10;
}
while(yu)
{
a[i][++len]=yu%10;
yu/=10;
}
for(j=len;j>=1;j--)
{
int t=a[i][j]+yu*10;
a[i][j]=t/(i+1);
yu=t%(i+1);
}
while(!a[i][len])
{
len--;
}
a[i][0]=len;
}
}
int main()
{
ktl();
int n;
while(scanf("%d",&n)!=EOF)
{
for(int i=a[n][0];i>0;i--)
{
printf("%d",a[n][i]);
}
puts("");
}
return 0;
}
组合数学
1排列
1.1不可重排列
1.2可重排列
1.3圆排列
1.4不尽相异元素全排列
1.5多重集的排列
2组合
2.1不可重组合数
2.2可重组合
2.3不相邻组合
2.4多重集的组合
2.5常用组合数公式
(1)C ( n , m ) = C ( n − 1 , m ) + C ( n − 1 , m − 1 )
(2)C ( n , m ) = C ( n , n − m )
(3)C ( n , m + 1 ) = ( n − m ) / ( m + 1 ) ∗ C ( n , m )
组合数
int n=30;
c[1][0]=c[1][1]=1;
for(register int i=2;i<=n;i++){
c[i][0]=1;
for(register int j=1;j<=i;j++)
c[i][j]=c[i-1][j]+c[i-1][j-1];
}
int C(int x, int y)
{
int ans = 1;
for (int i = 1; i <= y; i++)
{
ans *= (x - i + 1);
ans /= i;
}
return ans;
}
int quickpow(int a, int b) {
int res = 1;
while (b > 0) {
if (b & 1) res = res * a % mod;
b >>= 1;
a = a * a % mod;
}
return res;
}
int C(int n, int m) {
if (n < m) return 0;
int a = 1, b = 1;
m = min(m, n - m);
for (int i = 1; i <= m; i++) {
(a *= n - i + 1) %= mod;
(b *= i) %= mod;
}
return a % mod * quickpow(b, mod - 2) % mod;
}
int lucas(int n, int m) {
if (!m) return 1;
return C(n % mod, m % mod) * lucas(n / mod, m / mod) % mod;
}
第二类斯特林数
用递推求
n个互不相同的求，放到k个互补区分的盒子里，每个盒子里至少要有一个球，求方案数
(看最后一个球和前n-1个球的比较状态)
f(n,k) = k * f(n - 1,k) + f(n - 1,k - 1);
通常写为S2(n,k)
2.6组合数取模(模板) lucas
3常用公式及定理
3.1二项式定理
3.2鸽巢原理
将n+1个物品放到n个抽屉中，有一个抽屉至少有两个物品
3.3常见恒等式
3.4帕斯卡恒等式
3.5卢卡斯定理推论 exlucas
3.6容斥原理
3.7错排问题
4常见数列及其性质
4.1斐波那契数列
*4.2卡特兰数列 *
5递推方程
5.1线性递推方程
5.2非线性递推方程
5.3求解递推方程(模板)
6母函数
6.1普通母函数
6.2指数型母函数
6.3整数拆分
7Polya计数
8快速傅里叶FFT
exgcd拓欧
int exgcd(int a, int b, int &x, int &y)
{
if (b == 0) return x = 1, y = 0, a;
int res = exgcd(b, a%b, x, y);
int z = x;
x = y; y = z - y * (a / b);
return res;
}
int exgcd(int a,int b,int &x,int &y){
if(b == 0) return x = 1,y = 0,a;
int r = exgcd(b, a % b, x y);
tie(x,y) = make_tuple(y, x - (a / b) * y); //组成元组
/*
make_tuple(y, x - (a / b) * y):
int tmp = x;
x = y;
y = tmp - (a / b) * y;
*/
return;
}
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
LL a,b,t,x,y;//ax+by=t
void gcd(LL a,LL b,LL &t,LL &x,LL &y){
if(!b){t=a;x=1;y=0;}//边界：gcd(a,0)=1*a+0*0=a
else{gcd(b,a%b,t,y,x);y-=x*(a/b);}
}
int main(){
cin>>a>>b;
gcd(a,b,t,x,y);
if(t!=1)cout<<"No answer"<<endl;
else cout<<(x%b+b)%b<<endl;
return 0;
}
lucas卢卡斯定理-组合数取模
在m，n很大的时候,快速求解\(C_n^m\%p\)
（等价)
时间复杂度为\(O(p ⋅ l o g p m ⋅ l o g 2 p)\)
long long Lucas(long long n,long long m){
if(m==0) return 1;
return Lucas(n/p,m/p)*C(n%p,m%p)%p;
}
///
int fact[N];
int quickpow(int a, int b) {
int res = 1;
while (b > 0) {
if (b & 1) res = (res * a) % mod;
b >>= 1;
a = (a * a) % mod;
}
return res;
}
int getinv(int a) {
return quickpow(a, mod - 2);
}
int C(int x, int y) {
return fact[x] * quickpow(fact[y], mod - 2) % mod * quickpow(fact[x - y], mod - 2) % mod;
}
void init() {
fact[0] = 1;
for (int i = 1; i <= m; i++) {
fact[i] = fact[i - 1] * i % mod;
}
return;
}
long long C(long long n,long long m){
if(n<m) return 0;
if(m>n-m) m=n-m;
long long a=1,b=1;
for(int i=0;i<m;i++){
a=(a*(n-i))%p;
b=(b*(i+1))%p;
}
return a*quickpow(b,p-2)%p;
}
#include <iostream>
#include <algorithm>
#include <cstdio>
using namespace std;
long long n,m,p;
long long quickpow(long long base,long long power){
long long ret=1;
while(power){
if(power%2)
ret=ret*base%p;
base=base*base%p;
power/=2;
}
return ret;
}
long long C(long long n,long long m){
if(n<m) return 0;
if(m>n-m) m=n-m;
long long a=1,b=1;
for(int i=0;i<m;i++){
a=(a*(n-i))%p;
b=(b*(i+1))%p;
}
return a*quickpow(b,p-2)%p;
}
long long Lucas(long long n,long long m){
if(m==0) return 1;
return Lucas(n/p,m/p)*C(n%p,m%p)%p;
}
int main()
{
int T;
scanf("%d",&T);
while(T--){
scanf("%lld%lld%lld",&n,&m,&p);
printf("%lld\n",Lucas(n,m));
}
return 0;
}
exlucas拓展卢卡斯
扩展卢卡斯定理用于求如下式子（其中p pp不一定是质数）：
$ C_n^m\ mod\ p$
结合代码方便理解：
ll fac(const ll n, const ll p, const ll pk)
{
if (!n)
return 1;
ll ans = 1;
for (int i = 1; i < pk; i++)
if (i % p)
ans = ans * i % pk;
ans = power(ans, n / pk, pk);
for (int i = 1; i <= n % pk; i++)
if (i % p)
ans = ans * i % pk;
return ans * fac(n / p, p, pk) % pk;
}
层次二：组合数模质数幂
回到这个式子
\(\frac{\frac{n!} {p^{k1}}}{\frac{m!}{p^{k2}}\times \frac{(n-m)!}{p^{k3}}}\times p^{k1-k2-k3}\)
可以很容易地把它转换成代码（注意i要开long long）：
ll C(const ll n, const ll m, const ll p, const ll pk)
{
if (n < m)
return 0;
ll f1 = fac(n, p, pk), f2 = fac(m, p, pk), f3 = fac(n - m, p, pk), cnt = 0;
for (ll i = n; i; i /= p)
cnt += i / p;
for (ll i = m; i; i /= p)
cnt -= i / p;
for (ll i = n - m; i; i /= p)
cnt -= i / p;
return f1 * inv(f2, pk) % pk * inv(f3, pk) % pk * power(p, cnt, pk) % pk;
}
层次一：原问题
完整代码（题目：洛谷4720【模板】扩展卢卡斯）：
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
#include <climits>
#include <cmath>
using namespace std;
namespace zyt
{
const int N = 1e6;
typedef long long ll;
ll n, m, p;
inline ll power(ll a, ll b, const ll p = LLONG_MAX)
{
ll ans = 1;
while (b)
{
if (b & 1)
ans = ans * a % p;
a = a * a % p;
b >>= 1;
}
return ans;
}
ll fac(const ll n, const ll p, const ll pk)
{
if (!n)
return 1;
ll ans = 1;
for (int i = 1; i < pk; i++)
if (i % p)
ans = ans * i % pk;
ans = power(ans, n / pk, pk);
for (int i = 1; i <= n % pk; i++)
if (i % p)
ans = ans * i % pk;
return ans * fac(n / p, p, pk) % pk;
}
ll exgcd(const ll a, const ll b, ll &x, ll &y)
{
if (!b)
{
x = 1, y = 0;
return a;
}
ll xx, yy, g = exgcd(b, a % b, xx, yy);
x = yy;
y = xx - a / b * yy;
return g;
}
ll inv(const ll a, const ll p)
{
ll x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
}
ll C(const ll n, const ll m, const ll p, const ll pk)
{
if (n < m)
return 0;
ll f1 = fac(n, p, pk), f2 = fac(m, p, pk), f3 = fac(n - m, p, pk), cnt = 0;
for (ll i = n; i; i /= p)
cnt += i / p;
for (ll i = m; i; i /= p)
cnt -= i / p;
for (ll i = n - m; i; i /= p)
cnt -= i / p;
return f1 * inv(f2, pk) % pk * inv(f3, pk) % pk * power(p, cnt, pk) % pk;
}
ll a[N], c[N];
int cnt;
inline ll CRT()
{
ll M = 1, ans = 0;
for (int i = 0; i < cnt; i++)
M *= c[i];
for (int i = 0; i < cnt; i++)
ans = (ans + a[i] * (M / c[i]) % M * inv(M / c[i], c[i]) % M) % M;
return ans;
}
ll exlucas(const ll n, const ll m, ll p)
{
ll tmp = sqrt(p);
for (int i = 2; p > 1 && i <= tmp; i++)
{
ll tmp = 1;
while (p % i == 0)
p /= i, tmp *= i;
if (tmp > 1)
a[cnt] = C(n, m, i, tmp), c[cnt++] = tmp;
}
if (p > 1)
a[cnt] = C(n, m, p, p), c[cnt++] = p;
return CRT();
}
int work()
{
ios::sync_with_stdio(false);
cin >> n >> m >> p;
cout << exlucas(n, m, p);
return 0;
}
}
int main()
{
return zyt::work();
}
大整数gcd
//（a,b）代表ab公因数
莫比乌斯反演
莫比乌斯函数的定义
$u(i) = \begin{cases} 1, & \text{if \(n\) = 1} \ (-1)^k, & \text{if \(n = p_1*p_2*...*p_k\) } \0, & \text{其它} \end{cases}$
根据上面的定义，n大于1时且n是平方因子数时莫比乌斯函数值为0，否则从n的唯一分解定理中根据素数的个数取奇偶即可
莫比乌斯函数的性质
1)$\sum_{d|n}u(d) = \begin{cases} 1, & \text{if \(n\) = 1} \ \0, & \text{其它} \end{cases}$
有了上述这个式子，可以直接简单地以\(O(nlogn)\)筛出1到n内所有数的莫比乌斯函数值
//O(n)
bool vis[N];
int primes[N], miu[N];
int init(int n) {
int tot = 0;
miu[1] = 1;
for (int i = 2; i <= n; i++) {
if (!vis[i]) {
primes[tot++] = i;
miu[i] = -1;
}
for (int j = 0; j < tot; j++) {
int k = i * primes[j];
if (k > n)break;
vis[k] = true;
if (i % primes[j]) miu[k] = -miu[i];
else break;
}
}
}
//O(nlogn)
void init(int n) {
miu[1] = 1;
int t = n >> 1;
for (int i = 1; i <= t; i++) if (miu[i]) {
for (int j = i << 1; j <= n; j += i) miu[j] -= miu[i];
}
}
2)\(\sum_{d|n} \frac{u(d)}{d} = \frac{φ(n)}{n}\)
莫比乌斯函数除了可以用在莫比乌斯反演中之外，还可以用来进行容斥。
举一个常见的例子，求取1到n这个区间内有多少个数与x互质，一般的做法是直接进行容斥，但是我们可以发现容斥的系数刚好是莫比乌斯函数，即\(ans = \sum_{d|x} u(d) * \frac{n}{d}\)，其实这两者从本质考虑是完全等价的。
莫比乌斯反演是一个这样的式子：
定义\(F(n) = \sum_{d|n} f(n)\)，那么可以得到\(f(n) = \sum_{d|n} u(\frac{n}{d})F(d)\),莫比乌斯反演还有一种更常见的形式：\(F(n) = \sum_{n|d} f(d)\), 那么有\(f(n) = \sum_{n|d} u(\frac{d}{n})F(d)\)。
一般应用的都是上述的第二种形式，证明可以通过归纳得出，也可以直接通过式子变换得出，还可以由狄利克雷卷积证明。\(f(n) = \sum_{d|n} u(\frac{n}{d})F(d) = \sum_{d|n} u(\frac{n}{d})\sum_{x|d}f(x) = \sum_{d|n} f(d) \sum_{x|\frac{n}{d}}u(x) = f(n)\)
上述的证明中给出了一种常见的和式变换技巧：交换求和顺序。通过交换求和顺序，我们往往可以将某些和式化简或者更容易求出该和式，后面公式的化简将多次用到
莫比乌斯反演还有一个推广式
//应用举例：POJ 3904
//题目给出n和n个正整数，问能找出多少个不同的四元组(a,b,c,d)使得该四元组的最大公因数为1.
using namespace std;
#define _for(i,a,b) for(int i=(a);i<(b);++i)
#define _rep(i,a,b) for(int i=(a);i<=(b);++i)
typedef long long LL;
const int INF = 1 << 30;
const int MOD = 1e9 + 7;
const int maxn = 10005;
int n, a[maxn], tot[maxn];
int mu[maxn], vis[maxn];
int primes[maxn], cnt;
void get_mu() {
memset(vis, 0, sizeof(vis));
memset(mu, 0, sizeof(mu));
cnt = 0; mu[1] = 1;
for (int i = 2; i <= maxn; ++i) {
if (!vis[i]) { primes[cnt++] = i; mu[i] = -1; }
for (int j = 0; j<cnt&&primes[j] * i <= maxn; ++j) {
vis[primes[j] * i] = 1;
if (i%primes[j] == 0)break;
mu[i*primes[j]] = -mu[i];
}
}
}
void get_tot() {
memset(tot, 0, sizeof(tot));
for (int i = 0; i<n; ++i) {
int x = a[i];
int m = sqrt(x);
for (int j = 1; j <= m; ++j) {
if (x%j == 0)tot[j]++, tot[x / j]++;
}
if (m*m == x)tot[m]--;
}
}
LL Cn4(int m) {
if (m == 0)return 0;
return 1ll * m*(m - 1)*(m - 2)*(m - 3) / 24;
}
int main()
{
get_mu();
while (~scanf("%d", &n)) {
for (int i = 0; i<n; ++i) scanf("%d", &a[i]);
get_tot();
LL ans = 0;
for (int i = 1; i<maxn; ++i) {
ans += 1ll * mu[i] * Cn4(tot[i]);
}
printf("%I64d\n", ans);
}
return 0;
}
积性函数
设\(f(x)\)为一个数论函数，如果对于任意正整数\(a\)、\(b\)满足\((a,b)=1\),有\(f(ab)=f(a)∗f(b)\)的话，我们称\(f(n)\)为积性函数；如果对于任意正整数啊\(a、b\)有\(f(ab)=f(a)∗f(b)\)的话，我们称\(f(n)\)为完全积性函数。
常见的积性函数:
\(φ(n)\) －欧拉函数
\(μ(n)\) －莫比乌斯函数，关于非平方数的质因子数目
\(gcd(n,k)\) －最大公因子，当k固定的情况
\(d(n)\) －n的正因子数目
\(σ(n)\)－n的所有正因子之和
\(e(n)=[n==1]\)—若\(n = 1\)，\(ε(n)=1\)；若 \(n > 1\)，\(ε(n)=0\)。别称为“对于狄利克雷卷积的乘法单位”（完全积性）
\(I(n)=1\)恒等函数（完全积性）
\(id(n)=n\)单位函数（完全积性）
\(λ(n)\) －刘维尔函数，关于能整除\(n\)的质因子的数目
性质 若将\(n\)表示成质因子分解式\(n = p_1^{a1}p_2^{a2}...p_n^{an}\),则有\(f(n) = f(p_1^{a1})f(p_2^{a2})...f(p_n^{an})\),若\(f\)为积性函数且有\(f(p^n) = f^n(p)\)，则\(f\)为完全积性函数
积性函数前缀和
积性函数的和也是积性函数，积性函数前缀和是一个常见的问题，常常需要低于线性时间内解决。
狄利克雷卷积
狄利克雷卷积是解决积性函数前缀和问题的一个重要工具。
定义
对两个算术函数f, g，定义其Dirichlet卷积为新函数f * g，满足\((f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})\)。
狄利克雷卷积满足以下定律：
\(交换律 f * g = g * f\)
\(结合律 (f * g) * h = f * (g * h)\)
\(分配率 f * (g + h) = f * g + f * h\)
\(单位元^1:f * ∈ = f\)
若f，g均为积性函数，则f*g也是积性函数。
可以\(O(nlogn)\)预处理两个函数的Dirichlet卷积
LL f[N], g[N], h[N];
void calc(int n) {
for (int i = 1; i * i <= n; i++) {
h[i * i] += f[i] * g[i];
for (int j = i + 1; i * j <= n; j++) h[i * j] += f[i] * g[j] + f[j] * g[i];
}
}
求积性函数前缀和常用结论：
①\(\sum_{d|n}u(d) = [n==1]\)
②\(\sum_{d|n}φ(d) = n\)
③\(\sum_{i=1}^{n}i*[(i,n)==1] = \frac{n * φ(n) + [n==1]}{2}\)
④\(d(n^2) = \sum_{d|n} 2^{w(d)} = \sum_{d|n}\sum_{i|d}u^2(i)\),w(d)表示d的不同质因子的个数。
⑤\(d(n*m)= \sum_{i|n}\sum_{j|m}[(i,j)==1]\)
⑥\(\sum_{i=1}^{n}[(i,n)==1]*(i-1,n)=\sum_{d|n}\sum_{i=1}^{n}[(i,n)==1]*[d|(i-1)] = \sum_{d|n}φ(d)*\frac{φ(n)}{φ(d)}=d(n)*φ(n)\)
⑦\(d|i*j == \frac{d}{(i,d)}|j\)
矩阵与高斯消元
矩阵快速幂
#include <bits/stdc++.h>
using namespace std;
#define inf 0x3f3f3f3f
#define llinf 0x3f3f3f3f3f3f3f3f
#define int long long
#define ull unsigned long long
#define PII pair<int,int>
#define endl '\n'
const int N = 40;
const int mod = 1e9 + 7;
const double pi = acos(-1.0);
typedef long long ll;
int t, n, m, M;
struct Matrix {//矩阵初始化，乘法重载
int a[N][N];//2^5只有32，开40*40足够
Matrix()
{
for (int i = 0; i < m; i++) {
for (int j = 0; j < m; j++) {
a[i][j] = 0;
}
}
}
Matrix operator * (const Matrix& Ma_) const
{
Matrix res;
for (int i = 0; i < m; ++i) {
for (int j = 0; j < m; ++j) {
for (int k = 0; k < m; ++k) {
res.a[i][j] = (res.a[i][j] + a[i][k] * Ma_.a[k][j] % mod) % mod;
}
}
}
return res;
}
};
Matrix quickpow(Matrix res,Matrix sta, ll b)//快速幂板子
{
while (b > 0)
{
if (b & 1) res = res * sta;
sta = sta * sta;
b >>= 1;
}
return res;
}
void solve()
{
m = (1 << M);
Matrix state;
Matrix trans;
for (int i = 0; i < m; i++) {
state.a[i][i] = 1;
}
for (int i = 0; i < m; i ++) {
for (int j = 0; j < m; j ++) {
if ((i | j)) {
int flag = 0;
for (int k = 0; k <= M; k ++) {
if (((1 << k) & i) & ((1 << k) & j)) {
flag = 1;
break;
}
}
if (flag == 1)continue;
trans.a[i][j] = 1;
}
}
}
state = quickpow(state, trans, (n - 1));
int ans = 0;
for (int i = 0; i < m; i++) {
for (int j = 0; j < m; j++) {
ans = (ans + state.a[i][j] + mod) % mod;
}
}
cout << (ans + mod) % mod << endl;
return;
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
while (cin >> M >> n) {
solve();
}
return 0;
}
高斯消元 取模&不取模
#include <bits/stdc++.h>
using namespace std;
#define inf 0x3f3f3f3f
#define llinf 0x3f3f3f3f3f3f3f3f
#define int long long
#define ull unsigned long long
#define PII pair<int,int>
#define endl '\n'
const int N = 100 + 10;
const int mod = 1e9 + 7;
const double pi = acos(-1.0);
typedef long long ll;
using namespace std;
int t, n;
int a[N][N];
string s;
int quickpow(int a, int b) {//快速幂
int res = 1;
while (b > 0) {
if (b & 1) res = res * a % mod;
b >>= 1;
a = a * a % mod;
}
return res;
}
int getinv(int a) { return quickpow(a, mod - 2);}//费马小定理求逆元
int gauss(int n) {//高斯消元带模数
int res = 1;
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
if (a[j][i]) {
for (int k = i; k <= n; k++) {
swap(a[i][k], a[j][k]);
}
if (i != j) {
res = -res;
}
break;
}
}
if (!a[i][i]) return 0;
for (int j = i + 1, inv = getinv(a[i][i]); j <= n; j++) {
int t = a[j][i] * inv % mod;
for (int k = i; k <= n; k++) {
a[j][k] = (a[j][k] - t * a[i][k] % mod + mod) % mod;
}
}
res = (res * a[i][i] % mod + mod) % mod;
}
return res;
}
void init() {
for (int i = 0; i <= n; i++) {
for (int j = 0; j <= n; j++) {
a[i][j] = 0;
}
}
return;
}
void solve() {
cin >> n;
cin >> s;//由于题面说了det_a是一个不超过1e4长度的大数，数值肯定存不下，所以这里要先读入字符串再用模数处理一下
init();//初始化
int det_a = s[0] - '0';
for (int i = 1; i <= s.length() - 1; i++) {
det_a = ((det_a * 10) % mod + (s[i] - '0')) % mod;
}//处理大数
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cin >> a[i][j];
}
}//读入矩阵
int sum = gauss(n);//高斯消元计算矩阵的值
if (sum == det_a) cout << "+" << endl;//正数绝对值与自身相同，输出+
else cout << "-" << endl;//负数绝对值与自身不同，输出-
return;
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
while (cin >> t) {
while (t--) {
solve();
}
}
return 0;
}//2021ICPC济南J，高斯消元取模
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <iostream>
#define int long long
using namespace std;
typedef vector<int> v;
typedef vector<v> matrix;
matrix ans(30, v(31));//构造矩阵
int t;
void gauss() {
for (int i = 0; i < 30; i ++) {
int judge = i;
for (; judge < 30; judge++) {//查找非0项
if (ans[judge][i]) break;
}
swap(ans[i], ans[judge]);//交换
for (int j = 0; j < 30; j ++) {//消元
if (i != j && ans[j][i]) {//亦或0不改，亦或1要改
for (int k = 0; k <= 30; k ++) {
ans[j][k] = ans[i][k] ^ ans[j][k];
}
}
}
}
}
void solve() {
for (int i = 0; i < 5; i ++) {
for (int j = 0; j < 6; j ++) {
cin >> ans[i * 6 + j][30];
ans[i * 6 + j][i * 6 + j] = 1;
if (j >= 1) ans[i * 6 + j][i * 6 + j - 1] = 1;
if (i >= 1) ans[i * 6 + j][i * 6 + j - 6] = 1;
if (j < 5) ans[i * 6 + j][i * 6 + j + 1] = 1;
if (i < 4) ans[i * 6 + j][i * 6 + j + 6] = 1;
}
}
gauss();
for (int i = 0; i < 5; i ++) {
for (int j = 0; j < 6; j ++) {
if(j != 5) cout << ans[i * 6 + j][30] << " ";
else cout << ans[i * 6 + j][30] << endl;
}
}
return;
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
while (cin >> t) {
for(int i = 1;i <= t;i ++){
cout << "PUZZLE #" << i << endl;
solve();
}
}
return 0;
}//POJ1222 高斯消元板子
线性基
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define inf 0x3f3f3f3f
#define PII pair<int,int>
#define endl '\n'
const int N = 5e5 + 10;
const int mod = 1e9 + 7;
const double pi = acos(-1.0);
typedef long long ll;
int t, n, m;
int base[65];//最高位为第i位的基
void insert(int x) {
for (int i = 32; i >= 0; i--) {
if (x & (1ll << i)) {//如果i大于31,1后面要加ll
if (!base[i]) {
base[i] = x;
return;
}
x ^= base[i];
}
}
}
int find(int x) {
for (int i = 32; i >= 0; i--) {
if (x & (1ll << i)) {
if (!base[i]) return 0;
x ^= base[i];
}
}
return 1;
}
void init() {
for (int i = 0; i <= 64; i++) {
base[i] = 0;
}
return;
}
void solve() {
init();
int x, y;
for (int i = 1; i <= n; i++) {
cin >> x;
insert(x);
}
int q;
cin >> q;
while (q--) {
cin >> x >> y;
if (find(x ^ y)) cout << "YES" << endl;
else cout << "NO" << endl;
}
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
while (cin >> n) {
solve();
}
return 0;
}
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define inf 0x3f3f3f3f
#define PII pair<int,int>
#define endl '\n'
const int N = 2e5 + 10;
const int mod = 1e9 + 7;
const double pi = acos(-1.0);
typedef long long ll;
int t, n, m;
int arr[N], base[65];//最高位为第i位的基
int ans;
void insert(int x) {
for (int i = 60; i >= 0; i--) {
if (!(x >> i)) {
continue;
}
if (base[i] == 0) {
base[i] = x;
break;
}
else x ^= base[i];
}
if (x) ans++;
return;
}
void init() {
for (int i = 0; i <= 64; i++) {
base[i] = 0;
}
for (int i = 0; i <= n; i++) {
arr[i] = 0;
}
ans = 0;
return;
}
void solve() {
init();
cin >> arr[1];
insert(arr[1]);
int spj = arr[1];
for (int i = 2; i <= n; i++) {
cin >> arr[i];
spj ^= arr[i];
insert(arr[i]);
}
if (spj == 0) cout << -1 << endl;
else cout << ans << endl;
return;
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
while (cin >> n) {
solve();
}
return 0;
}
素数筛
//【普通筛——埃拉托斯特尼(Eratosthenes)筛法】(又称埃氏筛法)
bool is_prime[MAX+5];//标记是否为素数
int prime[MAX/2];//存放素数
int cnt=0;//已存放的素数个数
void isprime2() {
memset(is_prime,true,sizeof is_prime);
for(int i=2; i<=MAX; i++) {
if(is_prime[i]) {
prime[cnt++]=i;
for(int j=i+i; j<=MAX; j+=i) //二次筛选法:i是素数，则i的倍数必不为素数，筛掉
is_prime[j]=false;
}
}
}
//【线性筛——欧拉Euler筛】
prime[]数组中的素数是递增的,当i%prime[j]==0，那么i*prime[j+1]这个合数肯定被prime[j]乘以某个数筛掉。因为i中含有因子prime[j],prime[j]比prime[j+1]小，即i=k*prime[j]，那么i*prime[j+1]=（k*prime[j]）*prime[j+1]=k*（prime[j]*prime[j+1]），接下去的素数同理。所以不用筛下去了。因此，在满足i%prime[j]==0这个条件之前以及第一次满足改条件时,prime[j]必定是prime[j]*i的最小因子。
bool is_prime[MAX+5];
void isprime3() {
int is_prime[MAX+5];
int i,j,cnt=0;
memset(is_prime,true,sizeof(is_prime));
for(i=2; i<=MAX; i++) {
if(is_prime[i])
prime[cnt++]=i;
for(j=0; j<cnt&&prime[j]*i<=MAX; j++) {
is_prime[prime[j]*i]=false;
if(i%prime[j]==0) //保证每个合数只会被它的最小质因数筛去，因此每个数只会被标记一次
break;
}
}
}
离散化
#include<algorithm> // 头文件
//n 原数组大小 num 原数组中的元素 lsh 离散化的数组 cnt 离散化后的数组大小
int lsh[MAXN] , cnt , num[MAXN] , n;
for(int i=1; i<=n; i++) {
scanf("%d",&num[i]);
lsh[i] = num[i];
}
sort(lsh+1 , lsh+n+1);
cnt = unique(lsh+1 , lsh+n+1) - lsh - 1;
for(int i=1; i<=n; i++)
num[i] = lower_bound(lsh+1 , lsh+cnt+1 , num[i]) - lsh;
注意事项：
1.去重并不是把数组中的元素删去，而是重复的部分元素在数组末尾，去重之后数组的大小要减一
2.二分的时候，注意二分的区间范围，一定是离散化后的区间
3.如果需要多个数组同时离散化，那就把这些数组中的数都用数组存下来
第二种
第二种方式其实就是排序之后，枚举着放回原数组
用一个结构体存下原数和位置，按照原数排序
我结构体里面写了个重载，也可以写一个比较函数
最后离散化后数在$r a n k [ ] $里面
#include<algorithm>
struct Node {
int data , id;
bool operator < (const Node &a) const {
return data < a.data;
}
};
Node num[MAXN];
int rank[MAXN] , n;
for(int i=1; i<=n; i++) {
scanf("%d",&num[i].data);
num[i].id = i;
}
sort(num+1 , num+n+1);
for(int i=1; i<=n; i++)
rank[num[i].id] = i;
欧拉函数
求0 ~ x范围 x的因子之和
int eular(int x) {
int ans = x;
for (int i = 2; i * i <= x; i++) {
if (x % i == 0) {
ans = ans / i * (i - 1);//multiply
while (x % i == 0) x /= i;
}
}
if (x > 1) ans = ans / x * (x - 1);
return ans;
}
求φ1 ~ φn
const int N = ...;
bool not_prime[N];
int prime[N],tot = 0,phi[N];
for(int i = 2;i <= n;i ++){
if(!not_prime[i]){
prime[++tot] = i;
phi[i] = i - 1;
}
for(int j = 1;j <= tot && i * prime[j] <= n;j ++){
not_prime[i * prime[j]] = 1;
if(i % prime == 0){
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
欧拉定理
如果(a,m) = 1,那么a ^ φ(m) ≡1(mod m);
那么a ^ [φ(m) - 1] 为a 的逆元
线性同余方程&中国剩余定理
#include <iostream>
#include<cstdio>
typedef long long int ll;
using namespace std;
long long gcd,ans,temp;
int temp1;
ll exgcd(ll a, ll b, ll &x, ll &y)
{
if (b == 0)
{
x = 1;
y = 0;
return a;
}
ans=edgcd(b, a%b, x, y);
temp = x;
x = y;
y = temp - (a / b)*y;
return ans;
}
ll cal(ll a, ll b, ll c)//扩展欧几里得算法，用来求二元一次方程解；
{
ll b0;
ll x, y; gcd=edgcd(a, b, x, y);
if (c%gcd) {temp1=0;
}//但c与最大公因数不互质时，方程无解。
else {
b0 = b / gcd;
return (x*(c / gcd) % b0 + b0) % b0;
}
}
int main()
{
int n;
ll a1,a2,r1,r2,s;
while(cin>>n)
{
temp1=1;
scanf("%lld%lld",&a1,&r1);
n--;
while(n--)
{
scanf("%lld%lld",&a2,&r2);
s=cal(a1,a2,r2-r1);
r1=s*a1+r1;
a1=a1*a2/gcd;//旧的a1和a2的最小公倍数
}
if(temp1==0)r1=-1;
printf("%lld\n",r1);
}
}
Min_25求1~n质数和
const int N=1000010;
namespace Min25 {
int prime[N], id1[N], id2[N], flag[N], ncnt, m;
ll g[N], sum[N], a[N], T, n;
inline int ID(ll x) {
return x <= T ? id1[x] : id2[n / x];
}
inline ll calc(ll x) {
return x * (x + 1) / 2 - 1;
}
inline ll f(ll x) {
return x;
}
inline void init() {
//for(int i=0;i<=N;i++) prime[i]=id1[i]=id2[i]=flag[i]=g[i]=sum[i]=a[i]=0,ncnt=0,m=0;
ncnt=m=0;
T = sqrt(n + 0.5);
for (int i = 2; i <= T; i++) {
if (!flag[i]) prime[++ncnt] = i, sum[ncnt] = sum[ncnt - 1] + i;
for (int j = 1; j <= ncnt && i * prime[j] <= T; j++) {
flag[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
}
}
for (ll l = 1; l <= n; l = n / (n / l) + 1) {
a[++m] = n / l;
if (a[m] <= T) id1[a[m]] = m; else id2[n / a[m]] = m;
g[m] = calc(a[m]);
}
for (int i = 1; i <= ncnt; i++)
for (int j = 1; j <= m && (ll)prime[i] * prime[i] <= a[j]; j++)
g[j] = g[j] - (ll)prime[i] * (g[ID(a[j] / prime[i])] - sum[i - 1]);
}
inline ll solve(ll x) {
if (x <= 1) return x;
return n = x, init(), g[ID(n)];
}
}
int main()
{
ll n;
cin >>n;
cout <<Min25::solve(n)<<endl;
}
Meisell-Lehmer区间筛
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 5e6 + 2;
bool np[N];
int prime[N], pi[N];
int getprime() {
int cnt = 0;
np[0] = np[1] = true;
pi[0] = pi[1] = 0;
for(int i = 2; i < N; ++i) {
if(!np[i]) prime[++cnt] = i;
pi[i] = cnt;
for(int j = 1; j <= cnt && i * prime[j] < N; ++j) {
np[i * prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
return cnt;
}
const int M = 7;
const int PM = 2 * 3 * 5 * 7 * 11 * 13 * 17;
int phi[PM + 1][M + 1], sz[M + 1];
void init() {
getprime();
sz[0] = 1;
for(int i = 0; i <= PM; ++i) phi[i][0] = i;
for(int i = 1; i <= M; ++i) {
sz[i] = prime[i] * sz[i - 1];
for(int j = 1; j <= PM; ++j) {
phi[j][i] = phi[j][i - 1] - phi[j / prime[i]][i - 1];
}
}
}
int sqrt2(LL x) {
LL r = (LL)sqrt(x - 0.1);
while(r * r <= x) ++r;
return int(r - 1);
}
int sqrt3(LL x) {
LL r = (LL)cbrt(x - 0.1);
while(r * r * r <= x) ++r;
return int(r - 1);
}
LL getphi(LL x, int s) {
if(s == 0) 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 + 1;
if(x <= prime[s]*prime[s]*prime[s] && x < N) {
int s2x = pi[sqrt2(x)];
LL ans = pi[x] - (s2x + s - 2) * (s2x - s + 1) / 2;
for(int i = s + 1; i <= s2x; ++i) {
ans += pi[x / prime[i]];
}
return ans;
}
return getphi(x, s - 1) - getphi(x / prime[s], s - 1);
}
LL getpi(LL x) {
if(x < N) return pi[x];
LL ans = getphi(x, pi[sqrt3(x)]) + pi[sqrt3(x)] - 1;
for(int i = pi[sqrt3(x)] + 1, ed = pi[sqrt2(x)]; i <= ed; ++i) {
ans -= getpi(x / prime[i]) - i + 1;
}
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 - 2) * (b - a + 1) / 2;
for (int i = a + 1; 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 - 1);
}
}
return sum;
}
int main() {
init();
LL n;
while(cin >> n) {
cout << lehmer_pi(n) << endl;
}
return 0;
}
Miller-Rabin
#include <cstdlib>
#include <iostream>
using namespace std;
typedef long long ll;
ll qul_mul(ll a, ll b, ll n) {
ll num = 0;
while (b) {
if (b & 1)
num = (num + a) % n;
a = (a + a) % n;
b >>= 1;
}
return num;
}
ll qul_pow(ll a, ll b, ll n) //快速幂计算快速求出a^m^的值{
ll sum = 1;
while (b) {
if (b & 1)
sum = sum * a % n;
a = a * a % n;
b >>= 1;
}
return sum;
}
bool Miller_Rabin(ll n) {
int m = n - 1;
int t = 0;
if (n == 2)
return true;
else if (n < 2 || !(n & 1))
return false;
//求出n-1 = m*2^t的m和t。
while (!(m & 1)) {
m >>= 1;
t++;
}
for (int i = 0; i < 20; i++) {
//随机数a
ll a = rand() % (n - 1) + 1;
// 计算a^m
ll x = qul_pow(a, m, n), y;
for (int j = 0; j < t; j++) {
y = qul_mul(x, x, n); //进行(x*x)%n操作。
//不满足二次探测定理，也就是y得1了但是x并不等于1或者n-1，那么n就一定不是质数。
if (y == 1 && x != 1 && x != n - 1) {
return false;
}
x = y;
}
//上面循环跑完了，如果最后x都不等于1，那么一定是一个合数了。
if (x != 1) {
return false;
}
}
return true;
}
int main() {
ll n;
while (cin >> n) {
cout << Miller_Rabin(n) << endl;
}
return 0;
}
大素数
#include <cstdio>
#include <cstdlib>
#include <cstring>
#define MAXL 4
#define M10 1000000000
#define Z10 9
const int zero[MAXL - 1] = {0};
struct bnum {
int data[MAXL]; // 断成每截9个长度
// 读取字符串并转存
void read() {
memset(data, 0, sizeof(data));
char buf[32];
scanf("%s", buf);
int len = (int)strlen(buf);
int i = 0, k;
while (len >= Z10) {
for (k = len - Z10; k < len; ++k) {
data[i] = data[i] * 10 + buf[k] - '0';
}
++i;
len -= Z10;
}
if (len > 0) {
for (k = 0; k < len; ++k) {
data[i] = data[i] * 10 + buf[k] - '0';
}
}
}
bool operator==(const bnum& x) {
return memcmp(data, x.data, sizeof(data)) == 0;
}
bnum& operator=(const int x) {
memset(data, 0, sizeof(data));
data[0] = x;
return *this;
}
bnum operator+(const bnum& x) {
int i, carry = 0;
bnum ans;
for (i = 0; i < MAXL; ++i) {
ans.data[i] = data[i] + x.data[i] + carry;
carry = ans.data[i] / M10;
ans.data[i] %= M10;
}
return ans;
}
bnum operator-(const bnum& x) {
int i, carry = 0;
bnum ans;
for (i = 0; i < MAXL; ++i) {
ans.data[i] = data[i] - x.data[i] - carry;
if (ans.data[i] < 0) {
ans.data[i] += M10;
carry = 1;
} else {
carry = 0;
}
}
return ans;
}
// assume *this < x * 2
bnum operator%(const bnum& x) {
int i;
for (i = MAXL - 1; i >= 0; --i) {
if (data[i] < x.data[i]) {
return *this;
} else if (data[i] > x.data[i]) {
break;
}
}
return ((*this) - x);
}
bnum& div2() {
int i, carry = 0, tmp;
for (i = MAXL - 1; i >= 0; --i) {
tmp = data[i] & 1;
data[i] = (data[i] + carry) >> 1;
carry = tmp * M10;
}
return *this;
}
bool is_odd() { return (data[0] & 1) == 1; }
bool is_zero() {
for (int i = 0; i < MAXL; ++i) {
if (data[i]) {
return false;
}
}
return true;
}
};
void mulmod(bnum& a0, bnum& b0, bnum& p, bnum& ans) {
bnum tmp = a0, b = b0;
ans = 0;
while (!b.is_zero()) {
if (b.is_odd()) {
ans = (ans + tmp) % p;
}
tmp = (tmp + tmp) % p;
b.div2();
}
}
void powmod(bnum& a0, bnum& b0, bnum& p, bnum& ans) {
bnum tmp = a0, b = b0;
ans = 1;
while (!b.is_zero()) {
if (b.is_odd()) {
mulmod(ans, tmp, p, ans);
}
mulmod(tmp, tmp, p, tmp);
b.div2();
}
}
bool MillerRabinTest(bnum& p, int iter) {
int i, small = 0, j, d = 0;
for (i = 1; i < MAXL; ++i) {
if (p.data[i]) {
break;
}
}
if (i == MAXL) {
// small integer test
if (p.data[0] < 2) {
return false;
}
if (p.data[0] == 2) {
return true;
}
small = 1;
}
if (!p.is_odd()) {
return false; // even number
}
bnum a, s, m, one, pd1;
one = 1;
s = pd1 = p - one;
while (!s.is_odd()) {
s.div2();
++d;
}
for (i = 0; i < iter; ++i) {
a = rand();
if (small) {
a.data[0] = a.data[0] % (p.data[0] - 1) + 1;
} else {
a.data[1] = a.data[0] / M10;
a.data[0] %= M10;
}
if (a == one) {
continue;
}
powmod(a, s, p, m);
for (j = 0; j < d && !(m == one) && !(m == pd1); ++j) {
mulmod(m, m, p, m);
}
if (!(m == pd1) && j > 0) {
return false;
}
}
return true;
}
int main() {
bnum x;
x.read();
puts(MillerRabinTest(x, 5) ? "Yes" : "No");
return 0;
}
佩尔方程
如果D是一个正整数且不是完全平方数， xx - Dy*y = 1 总是有正整数解
欧拉降幂
我们要求 $a^{b}% p $，如果 b 是个比较大的数（可能是个幂的形式），可以先让 b 对 (p−1)取模，是不影响最终答案的正确性的（前提是 a 和 p 互素）。这个取模就叫做欧拉降幂
随机数
#include <iostream>
#include <chrono>
#include <random>
using namespace std;
int main()
{
// 随机数种子
unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
mt19937 rand_num(seed); // 大随机数
cout << rand_num() << endl;
return 0;
}//无范围
#include <iostream>
#include <chrono>
#include <random>
using namespace std;
int main()
{
// 随机数种子
unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
mt19937 rand_num(seed); // 大随机数
uniform_int_distribution<long long> dist(0, 1000000000); // 给定范围
cout << dist(rand_num) << endl;
return 0;
}//有范围
