数值算法

拉格朗日插值（Lagrange 插值法）

P4781 【模板】拉格朗日插值

用于构造一个函数 \(f(x)\) 使其过点 \(P_1(x_1, y_1), P_2(x_2,y_2),\cdots,P_n(x_n,y_n)\).

首先设第 \(i\) 个点在 \(x\) 轴上的投影为 \(P_i^{\prime}(x_i,0)\).

考虑构造 \(n\) 个函数 \(f_1(x), f_2(x), \cdots, f_n(x)\)，使得对于第 \(i\) 个函数 \(f_i(x)\)，其图像过 \(\begin{cases}P_j^{\prime}(x_j,0),(j\neq i)\\P_i(x_i,y_i)\end{cases}\)，则可知题目所求的函数 \(f(x)=\sum\limits_{i=1}^nf_i(x)\).

那么可以设 \(f_i(x)=a\cdot\prod_{j\neq i}(x-x_j)\)，将点 \(P_i(x_i,y_i)\) 代入可以知道 \(a=\dfrac{y_i}{\prod_{j\neq i} (x_i-x_j)}\)，所以 \(f_i(x)=y_i\cdot\dfrac{\prod_{j\neq i} (x-x_j)}{\prod_{j\neq i} (x_i-x_j)}=y_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}.\)

那么我们就可以得出 Lagrange 插值的形式为：\(f(x)=\sum_{i=1}^ny_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}\)

朴素实现的时间复杂度为 \(O(n^2)\)，可以优化到 \(O(n\log^2 n)\)，需用多项式快速插值。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2e3+10,p=998244353;
ll x[N],y[N],ans;
ll ksm(ll a,ll b){
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}
int main(){
	int n;
	ll k;
	cin>>n>>k;
	for(int i=1;i<=n;i++) cin>>x[i]>>y[i];
	for(int i=1;i<=n;i++){
		ll sum=y[i]%p;
		for(int j=1;j<=n;j++){
			if(i==j) continue;
			sum=sum*(k-x[j]+p)%p*ksm((x[i]-x[j]+p)%p,p-2)%p;
		}
		ans=(ans+sum)%p;
	}
	cout<<ans;
	return 0;   
} 

FFT（快速傅里叶变换）

一、DFT（离散傅里叶变换）

二、FFT（快速傅里叶变换）

#include<bits/stdc++.h>
using namespace std;
const int N=4e6+10;
const double pi=acos(-1);//计算 π 的值
struct comp{
	double x,y;
	comp(double xx=0,double yy=0){//代入变量则为
		x=xx,y=yy;
	}
	comp(double rho,double k,double n){//算单位圆上坐标
		x=rho*cos(2*pi*k/n);
		y=rho*sin(2*pi*k/n);
	}
}a[N],b[N],c[N];
comp operator + (const comp &a,const comp &b){return comp(a.x+b.x,a.y+b.y);}//定义复数加减乘运算
comp operator - (const comp &a,const comp &b){return comp(a.x-b.x,a.y-b.y);}
comp operator * (const comp &a,const comp &b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int rev[N];
void FFT(comp *a,int len,int type){
	int wei=30-__builtin_clz(len);//为了快速变换，求与len值最接近的2的次幂
	for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<wei);//翻转
	for(int i=0;i<len;i++) if(rev[i]<i) swap(a[i],a[rev[i]]);
	for(int k=1;k*2<=len;k*=2){
		comp constomega=comp(1,type,2*k);//
		for(int j=0;j<len;j+=2*k){
			comp omega=comp(1,0);
			for(int i=0;i<k;i++){
				comp u=a[j+i],v=a[j+i+k];
				a[j+i]=u+omega*v;
				a[j+i+k]=u-omega*v;
				omega=omega*constomega;
			}
		}
	}
	if(type==-1) for(int i=0;i<len;i++) a[i].x/=len,a[i].y/=len;
}
void poly_mul(comp *a,int n,comp *b,int m,comp *c){//a*b
	int len=n+m;
	len=1<<(1+(int)log2(len-1));
	FFT(a,len,1);FFT(b,len,1);
	for(int i=0;i<len;i++) c[i]=a[i]*b[i];
	FFT(c,len,-1);
}
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	cout.tie(0);
	int n,m;
	cin>>n>>m;
	for(int i=0;i<=n;i++) cin>>a[i].x;
	for(int i=0;i<=m;i++) cin>>b[i].x;
	poly_mul(a,n,b,m,c);
	for(int i=0;i<=n+m;i++) cout<<int(c[i].x+0.5)<<' ';
	return 0;
} 

原根

P6091 【模板】原根

定义：原根是一些特殊元素，它的阶就等于所有模 \(m\) 简化剩余系的个数．

可用原根遍历完剩余系中的所有数，可用于 NTT .

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+10;
vector<int> num,ans;
int phi[N],vis[N],prime[N],cnt;
void Phi(int n){//求phi[1]-phi[n]
	phi[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]){
			phi[i]=i-1;
			prime[++cnt]=i;
		}
		for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0){
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			else phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}
}
bool check(int m){//先判断数是否有原根，根据证明可得：形似 2,4，p^k,(2p)^k（其中p为质数）这样的数有原根
	if(m==2||m==4) return true;
	if(m%2==0) m/=2;
	for(int i=2;prime[i]<=m;i++){
		if(m%prime[i]==0){
			while(m%prime[i]==0) m/=prime[i];
			return m==1;
		}
	}
	return false;
}
void div(int n){//求n的因数
	for(int i=2;i*i<=n;i++){
		if(n%i==0){
			num.push_back(i);
			while(n%i==0) n/=i;
		}
	}
	if(n>1) num.push_back(n);
}
int gcd(int a,int b){//手动求gcd
	if(b) return gcd(b,a%b);
	return a;
}
ll ksm(ll a,ll b,ll p){
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}
int main(){
	int T,n,d,x,siz;
	cin>>T;
	Phi(1000000);
	while(T--){
		cin>>n>>d;
		if(!check(n)){
			puts("0\n");
			continue;
		}
		num.clear(),ans.clear();
		div(phi[n]);
		for(int i=1;;i++){//暴力枚举i使其与n互质
			if(gcd(i,n)!=1) continue;
			int flag=1;
			for(int j=0;j<num.size();j++){
				if(ksm(i,phi[n]/num[j],n)==1){//找第一个原根：i^(任意phi[n]也就是群的阶（群的大小）的因数次)≠1的就是原根
					flag=0;
					break;
				}
			}
			if(flag){//找到第一个原根就退出循环
				x=i;
				break;
			}
		}
		ll y=1;
		siz=phi[phi[n]];//原根总个数
		for(int i=1;ans.size()<siz;i++){//求x^k是否为原根
			y=x*y%n;//计算x^i%n
			if(gcd(phi[n],i)==1) ans.push_back(y);//若gcd(phi[n],i)=1，则x^i为原根
		}
		sort(ans.begin(),ans.end());
		cout<<siz<<endl;
		for(int i=d-1;i<siz/d*d;i+=d) printf("%d ",ans[i]);//按题目要求输出
		puts("");
	}
	return 0;
}

NTT（快速数论变换）

P3803 【模板】多项式乘法（FFT）

\(e^{\frac{2\pi i}{n}\times k}=原根^{\frac{mod-1}{n}\times k}\)

#include<bits/stdc++.h>
using namespace std;
const int N=4e6+10,p=998244353;
const int G=3,Gi=332748118;//3是p的原根 
int a[N],b[N],c[N],d[N],ans[N];
int rev[N],ta[N],tb[N];
int ksm(int a,int b){
    int res=1;
    while(b){
        if(b&1) res=1ll*res*a%p;
        a=1ll*a*a%p;
        b>>=1;
    }
    return res;
}
void NTT(int *a,int len,int type){
    int wei=__lg(len);
    for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(wei-1));
    for(int i=0;i<len;i++) if(rev[i]<i) swap(a[i],a[rev[i]]);
    for(int k=1;k<len;k<<=1){
        int wn=ksm(type==1?G:Gi,(p-1)/(k<<1));//判断是用本身还是用逆元 
        for(int j=0;j<len;j+=k<<1){
            int w=1;
            for(int i=0;i<k;i++){
                int u=a[j+i],v=1ll*w*a[j+i+k]%p;
                a[j+i]=(u+v)%p;
                a[j+i+k]=(u-v+p)%p;
                w=1ll*w*wn%p;
            }
        }
    }
    if(type==-1){
        int inv_len=ksm(len,p-2);
        for(int i=0;i<len;i++) a[i]=1ll*a[i]*inv_len%p;
    }
}
void poly_mul(int *a,int n,int *b,int m,int *c){
    int len=1;
    while(len<=n+m) len<<=1;
    NTT(a,len,1),NTT(b,len,1);
    for(int i=0;i<len;i++) c[i]=1ll*a[i]*b[i]%p;
    NTT(c,len,-1);
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    int n,m;
    cin>>n>>m;
    for(int i=0;i<=n;i++) cin>>a[i];
    for(int i=0;i<=m;i++) cin>>b[i];
    poly_mul(a,n,b,m,c);
    for(int i=0;i<=n+m;i++) cout<<c[i]<<' ';
    return 0;
}

积分

定积分的定义
简单来说，函数 \(f(x)\) 在区间 \([l,r]\) 上的定积分 \(\int_{l}^{r}f(x)\mathrm{d}x\) 指的是 \(f(x)\) 在区间 \([l,r]\) 中与 \(x\) 轴围成的区域的面积（其中 \(x\) 轴上方的部分为正值，\(x\) 轴下方的部分为负值）．

很多情况下，我们需要高效，准确地求出一个积分的近似值．下面介绍的 辛普森法，就是这样一种求数值积分的方法．

辛普森法

P4525 【模板】自适应辛普森法 1

这个方法的思想是将被积区间分为若干小段，每段套用二次函数的积分公式进行计算．

对于一个二次函数 \(f(x)=ax^2+bx+c\)，有：

\[\int_l^r f(x) {\mathrm d}x = \frac{(r-l)(f(l)+f(r)+4 f(\frac{l+r}{2}))}{6} \]

推导过程： 对于一个二次函数 \(f(x)=ax^2+bx+c\)；求积分可得 \(F(x)=\int_0^x f(x) {\mathrm d}x = \frac{a}{3}x^3+\frac{b}{2}x^2+cx+D\) 在这里 \(D\) 是一个常数，那么

\[\begin{aligned} \int_l^r f(x) {\mathrm d}x &= F(r)-F(l) \\ &= \frac{a}{3}(r^3-l^3)+\frac{b}{2}(r^2-l^2)+c(r-l) \\ &=(r-l)(\frac{a}{3}(l^2+r^2+lr)+\frac{b}{2}(l+r)+c) \\ &=\frac{r-l}{6}(2al^2+2ar^2+2alr+3bl+3br+6c)\\ &=\frac{r-l}{6}((al^2+bl+c)+(ar^2+br+c)+4(a(\frac{l+r}{2})^2+b(\frac{l+r}{2})+c)) \\ &=\frac{r-l}{6}(f(l)+f(r)+4f(\frac{l+r}{2})) \end{aligned}\]

自适应辛普森法

由于普通的方法为保证精度在时间方面会受到 \(n\) 的限制，所以我们用一种更加合适的方法．

若有一段图像已经很接近二次函数的话，直接带入公式求积分，得到的值精度就很高了，不需要再继续分割这一段．

于是采用这样的分割方法：每次判断当前段和二次函数的相似程度，如果足够相似的话就直接代入公式计算，否则将当前段分割成左右两段递归求解．

如何判断每一段和二次函数是否相似

我们把当前段代入公式求积分，再将当前段从中点分割成两段，把这两段再直接代入公式求积分．

如果当前段的积分和分割成两段后的积分之和相差很小的话，就可以认为当前段和二次函数很相似，不用再递归分割．

上面就是自适应辛普森法的思想．在分治判断的时候，除了判断精度是否正确，一般还要强制执行最少的迭代次数．

代码如下：

#include<bits/stdc++.h>
using namespace std;
double a,b,c,d;
double f(double x){
	return (c*x+d)/(a*x+b);
}
double simpson(double l,double r){
	double mid=(l+r)/2;
	return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double eps,double ans){
	double mid=(l+r)/2;
	double x=simpson(l,mid),y=simpson(mid,r);
	if(fabs(x+y-ans)<=15*eps) return x+y+(x+y-ans)/15;
	return asr(l,mid,eps/2,x)+asr(mid,r,eps/2,y);
}
int main(){
	double l,r;
	cin>>a>>b>>c>>d>>l>>r;
	double eps=1e-7;
	double ans=asr(l,r,eps,simpson(l,r));
	printf("%.6lf",ans);
	return 0;
}