算法中数学相关内容的整理和补齐
数学算法
此专栏持续更新，可能要等到猴年马月更新结束
1.快速幂
算法出现背景：
给出 a,b, 要求出 ab mod p,并且这时的 a,b 往往又很大，直接求会爆炸，所以快速幂应运而生，事实
上，快速幂诞生于上个世纪六七十年代（？）
核心思想：
1.我们在计算随意一个数的幂次方时，计算过程中会产生许多对我们后续计算有用的数，我们可以将它们
在过程中保存下来，减少总共的计算次数。
比方说，我们在计算 44 时，你肯定不会连续把 4 乘以 4 次，而是在算出 42 等于 16 时，直接将 16 进行
平方得到256，想一想，这样是不是比直接计算 44 更快更有效？
2.对于取余运算，我们又有如下性质：
证明以后补上
因此我们可以在计算过程中不断取模，而不是最后对一个很大的数进行取模，这便是快速幂算法另一
方面的优化。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll a,b,p;
int quickpower(ll a,ll b)
{
int ans=1,base=a;
while(b)
{
if(b&1)
{
ans*=base;
ans%=p;
}
base*=base;
base%=p;
b>>=1;
}
return ans;
}
int main()
{
cin>>a>>b>>p;
int ans;
ans=quickpower(a,b);
cout<<ans<<endl;
}
我们可以这样理解代码，base 是一个准备好的基，一开始为 a ，然后不断升级，变为 a2 , a4, ......
而 ans 通过二进制来询问每一位是否为 1 ，如果为 1 这时准备好的 base 便可以排上用场。
2.乘法逆元
算法出现背景：
a 模 p 的乘法逆元定义为 ax \(\equiv\) 1 (mod p) 的解，因为在加减乘除四种运算中，只有除法不满足一般的
性质，因此人们想出了逆元来解决这一问题。
求解逆元的方式以及各自的核心思想：
首先说一下存在逆元的充要条件：a 与 p 互质
1.拓展欧几里得
这一思想就是把求逆元的问题转化为线性同余方程 a * x \(\equiv\) c (mod b)在 c = 1 时的解，因此我们可以
先考虑方程 ax + by =1 。
对于一般情况下 ax + by = c，求整数解（x ， y )，这时我们要用到朴素欧几里得。
对于上述方程，存在整数解的充要条件是 gcd(a,b) | c
对于求解gcd(a,b),我们应用朴素欧几里得，朴素欧几里得如下:
gcd(a,b)=gcd(b,a mod b)
由此我们得到了计算 gcd 的函数
点击查看代码
int gcd(int a, int b) {
return !b ? a : gcd(b, a % b);
}
扩展欧几里得代码：
点击查看代码
void Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) x = 1, y = 0;
else Exgcd(b, a % b, y, x), y -= a / b * x;
}
完整代码：
点击查看代码
void Exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) x = 1, y = 0;
else Exgcd(b, a % b, y, x), y -= a / b * x;
}
int main() {
ll x, y;
Exgcd (a, p, x, y);
x = (x % p + p) % p;
printf ("%d\n", x); //x是a在mod p下的逆元
}
特点：速度比较快，但只在单个查找时效率比较高。
2.快速幂
这一算法利用了费马小定理。
代码如下：
点击查看代码
ll fpm(ll x, ll power, ll mod) {
x %= mod;
ll ans = 1;
for (; power; power >>= 1, (x *= x) %= mod)
if(power & 1) (ans *= x) %= mod;
return ans;
}
int main() {
ll x = fpm(a, p - 2, p); //x为a在mod p意义下的逆元
}
3.递推算法
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,p;
long long inv[20000529];
int main()
{
scanf("%d%d",&n,&p);
inv[1]=1;
inv[0]=0;
printf("1\n");
for(long long i=2;i<=n;++i)
{
inv[i]=p-(p/i)*inv[p%i]%p;
printf("%d\n",inv[i]);
}
}
特点：对于连续查找多个值非常好用。
3.多项式乘法(FFT)
背景：
用来解决多项式乘以多项式的系数问题，如果直接把两个多项式的系数取出来对应相乘，我们会发现
复杂度是 \(O(n^2)\) 的，那么还有没有办法优化呢？
事实上，多项式是有两种表现方式的，除了我们常见的系数表示法之外，还有点表示法。
点表示法：
最简单地，我们都知道两点确定一条直线，三点确定一条抛物线，对于任意的 n 次多项式，只要我们
知道它上面 n + 1 个点的坐标，就可以通过这些点来表示它。
也许到这里，你仍然不知道它对优化时间复杂度有什么帮助，但至少，你能感受到，它相对于系数表示
法显得更加灵活，毕竟系数是确定的，但这些点却是可以任选的。
优化的实现：
伟大的人类在研究此问题时想到了取 n 个复数单位根，把 0 到 n-1 次幂的单位根带入，
将多项式 \(P(x)\) 按照下标的奇偶性分类：
所以\(A(x)=A_{1}(x^2)+xA_{2}(x^2)\)
我们将 \(\omega_{n}^k\) 和 \(\omega_{n}^{\frac{n}{2}+k}\) 分别代入，会发现两个式子相差的只有常数项。因此我们就将问题缩小了一半。
递归地想，我们可以将问题的规模不算缩小，直到多项式只剩下常数项。
这样的复杂度为 \(O(n logn)\)
然而，还有最后一点，点值表示法在我们生活中很少被用到，所以我们还需要利用傅里叶逆变换将
点值表示法变回系数表示法。
可得两者之间的关系为：
所以总体思路就是：
系数表示法 -> 点值表示法 -> 系数表示法
利用二进制的数位反转，还可以简化奇偶分类这一过程。
4.中国剩余定理
伟大的中国古人发现的算法耶！
因为在之前学习基础数论的时候已经对相关内容做过了解，并且大多时候也直接用板，所以暂时不作
具体理解的介绍。（虽然我其实也还没完全明白）
点击查看代码
#include<bits/stdc++.h>
using namespace std;
long long a[15],mul=1,b[15],Mi[15],X;
int n;
void exgcd(long long a,long long b,long long&x,long long&y)
{
if(b==0){
x=1;
y=0;
return;
}
exgcd(b,a%b,x,y);
int z=x;
x=y;
y=z-y*(a/b);
}
int main()
{
cin>>n;
for(int i=1;i<=n;i++)
{
cin>>a[i];
mul*=a[i];
cin>>b[i];
}
for(int i=1;i<=n;i++)
{
Mi[i]=mul/a[i];
long long x=0,y=0;
exgcd(Mi[i],a[i],x,y);
X+=b[i]*Mi[i]*(x<0?x+a[i]:x);
}
cout<<X%mul<<endl;
return 0;
}
过程中需要用到扩展欧几里得求逆元。
5.裴蜀定理
裴蜀定理内容：
\(ax+by=c,x\in Z，y\in Z\)成立的充分必要条件是 \(gcd(a,b)|c\).
这个定理是可以推广到 n 个数的。
所以对于洛谷的这道模板题，我们只需要求出这些数的 gcd (注意对于负数要先变为正数)
然后最小值便就是 gcd 了。
代码：
点击查看代码
#include<bits/stdc++.h>
using namespace std;
inline int gcd(int x,int y)
{
return y?gcd(y,x%y):x;
}
int main()
{
int n;
cin>>n;
int ans=0,tmp;
for(int i=1;i<=n;i++)
{
cin>>tmp;
if(tmp<0)tmp=-tmp;
ans=gcd(ans,tmp);
}
cout<<ans<<endl;
return 0;
}
6.自适应辛普森法
原理：
辛普森公式是在积分区间 \([a,b]\) 上选取 \(a,b,mid=\frac{a+b}{2}\) 这三个点，然后通过这三个点构建出一条
抛物线来拟合原函数 \(f(x)\)。
代码：
点击查看代码
inline double simpson(double l,double r) {
double mid=(l+r)/2;
return (f(l)+4*f(mid)+f(r))*(r-l)/6;
}
自适应辛普森法的“自适应”就体现能根据精度需求，调整区间分割的大小，分割的过程就是不断二分，
满足精度要求后就终止。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
double a,b,c,d,l,r;
inline double f(double x)
{
return (c*x+d)/(a*x+b);
}
inline double simpson(double l,double r)
{
double mid=(l+r)/2;
return(f(l)+4*f(mid)+f(r))*(r-l)/6;
}
double asr(double l,double r,double eps,double ans)
{
double mid=(l+r)/2;
double l1=simpson(l,mid),r1=simpson(mid,r);
if(fabs(l1+r1-ans)<=15*eps)return l1+r1+(l1+r1-ans)/15;
return asr(l,mid,eps/2,l1)+asr(mid,r,eps/2,r1);
}
inline double asr(double l,double r,double eps)
{
return asr(l,r,eps,simpson(l,r));
}
int main()
{
cin>>a>>b>>c>>d>>l>>r;
cout<<fixed<<setprecision(6)<<asr(l,r,1e-6);
return 0;
}
7.卢卡斯定理
应用：
解决组合数取模问题。
代码：
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int p;
ll a[100010];
ll pow(ll y,int z,int p)
{
y%=p;
ll ans=1;
for(int i=z;i;i>>=1,y=y*y%p)
{
if(i&1)ans=ans*y%p;
}
return ans;
}
ll C(ll n,ll m)
{
if(m>n)return 0;
return ((a[n]*pow(a[m],p-2,p))%p*pow(a[n-m],p-2,p)%p);
}
ll lucas(ll n,ll m)
{
if(!m)return 1;
return C(n%p,m%p)*lucas(n/p,m/p)%p;
}
int main()
{
int t;
cin>>t;
while(t--)
{
int n,m;
cin>>n>>m>>p;
a[0]=1;
for(int i=1;i<=p;i++)
{
a[i]=(a[i-1]*i)%p;
}
cout<<lucas(n+m,n)<<endl;
}
}
浙公网安备 33010602011771号