数学知识啦啦啦

Bold!

已于 2022-08-16 15:24:54 修改

阅读量208

点赞数

文章标签：算法

于 2021-09-12 20:41:42 首次发布

本文链接：https://blog.csdn.net/srh20/article/details/119805304

版权

数学知识模板：https://www.acwing.com/blog/content/406/

Practice:

1.线性筛：质数距离

https://www.acwing.com/problem/content/description/198/

**考察：**线性筛+巧妙优化

思路：
因为l,r<=2^31,所以不能直接筛出1-r的所有素数。
因为任何一个合数n必定包含一个不超过sqrt(n)的质因子，又sqrt(2^31)=50000,所以先预处理，用筛法筛出2-50000之间的所有质数。
因为质数的倍数是合数，所以对于每个1-50000内的质数p，筛掉[L,R]范围内的合数（至少为2*p)，剩下的就是质数。

步骤
1、找出1 - 50000中的所有质因子
2、对于1 ~ 50000 中每个质数p，将[L,R]中所有p的倍数筛掉(至少2倍)
找到大于等于L的最小的p的倍数x，找下一个倍数时只需+= p

key:
ceil(L/P)=floor((L+P-1)/P)

AC代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e6+10;
int prime[N],num;
bool v[N];
void get(int n){//线性筛 
	memset(v,0,sizeof(v));
	num=0;
	for(int i=2;i<=n;i++){
		if(!v[i]) prime[num++]=i;
		for(int j=0;prime[j]*i<=n;j++){
			v[i*prime[j]]=true;
			if(i%prime[j]==0) break;
		}
	}
}
int main(){
	int l,r;
	while(cin>>l>>r){
		get(50000);//因为l,r<=2^31,所以sqrt(2^31)=50000,所以先预处理1- 50000的质数 
		memset(v,0,sizeof(v));//清false 
		for(int i=0;i<num;i++){
			ll x=prime[i];//质数 
			for(ll j=max(2*x,(l+x-1)/x*x);j<=r;j+=x){//筛掉范围内质数的倍数 ，因为l,r很小时，可能之间会包括质数x，所以至少从2*x开始 
				v[j-l]=true;
			}
		}
		num=0;
		for(int i=0;i<=r-l;i++){//得到l,r之间的质数 
			if(!v[i]&&i+l>1) prime[num++]=i+l;
		}
		if(num<2) puts("There are no adjacent primes.");
		else{
			int mi=0,ma=0;
			for(int i=0;i+1<num;i++){//因为要让i+1有意义 
				int d=prime[i+1]-prime[i];
				if(d<prime[mi+1]-prime[mi]) mi=i;
				if(d>prime[ma+1]-prime[ma]) ma=i;
			}
			printf("%d,%d are closest, %d,%d are most distant.\n",
			prime[mi],prime[mi+1],
			prime[ma],prime[ma+1]);
		
		}
	}
	
	return 0;
}

2.质数筛法：阶乘分解

https://www.acwing.com/problem/content/description/199/

考察：质数筛法，转化。

质数筛法模板：

void primes(int n){
	memset(v,0,sizeof(v));
	for(int i=2;i<=n;i++){
		if(!v[i]){
			prime[++num]=i;
			for(int j=i;j<=n/i;j++){
				v[i*j]=1;
			} 
		}
	}
}

思路：因为N！的每个质因子都不会超过N，所以使用素数筛找出1-N所有的质数，然后
再考虑N!中一共包含多少个质因子p,对于每一个p,N!所含p的数目为N/P+N/(P^2)+N/(P3)+``````+N(P^K)

代码如下：

#include <bits/stdc++.h>
using namespace std;
const int N = 1e6+10;
typedef long long ll;
ll prime[N],num,n;
bool v[N];
ll qk(ll a,ll b){//快速幂 
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a;
		a*=a;
		b>>=1;
	}
	return ans;
}
void get(ll n){//质数筛法 
	num=0;
	for(ll i=2;i<=n;i++){
		if(!v[i]) prime[++num]=i;
		for(ll j=i;j<=n/i;j++){
			v[i*j]=1;
		}
	}
}
int main(){
	cin>>n;
	get(n);
	for(ll i=2;i<=n;i++){
		if(!v[i]){
			ll ans=0;
			for(ll j=1;qk(i,j)<=n;j++){
				ans+=n/qk(i,j);
			}
			cout<<i<<" "<<ans<<endl;
		}
	}
	return 0;
}

3.反素数

原题链接：https://www.acwing.com/problem/content/description/200/

引理：
1.1-N中最大的反质数，就是1-N中约数个数最多的数中的最小的一个。
2.1-N中任何数的质因子都不会超过10给，且所有质因子指数总和不超过30。
3.任意x属于1-N,x为反质数的必要条件是:x的质因子是连续的若干个最小的质数，并且质数单调递减。

思路：
采用dfs,从小到达枚举质因数，确定其个数，并且保证每一个的指数都小于上一个质因数的指数。

代码如下：

在这里插入代码片#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int INF=1e9;
ll zs[9]={2,3,5,7,11,13,17,19,23};
ll n,maxd,num;
void dfs(ll u,ll last,ll p,ll s){//当前质因数，上一质因数次数，当前数字，约数个数 
	if(s>maxd||(s==maxd&&p<num)){//如果约数个数最多或者约数个数相等但是数字更小，那么更新maxd和num 
		maxd=s;//约数最多为s 
		num=p;//最优解是数字p 
	}
	if(u==9) return;//根据引理2，质因子个数不会超过10个 
	for(ll i=1;i<=last;i++){//枚举次数，根据指数单调递减的性质，不能超过上一个质因数的次数last
		if(p*zs[u]>n) break;//如果当前数乘以该质因数>n ，直接break; 
		p*=zs[u];//否则，继续累乘 
		dfs(u+1,i,p,s*(i+1)); //往深处搜索，下一步的质因数为zs[u+1],上一质因数次数为i,当前数字为p,约数个数是s*(i+1) 
	}

}
int main(){
	cin>>n;
	dfs(0,30,1,1);//最开始的质因数是zs[0],次数最大是30，2^30,上一个数为1，约数个数为1（只有2） 
	cout<<num<<endl;
	return 0;
}

4.除法分块：AcWing199. 余数之和

原题链接：https://www.acwing.com/problem/content/description/201/

除法分块：
这个讲得挺好：https://blog.nowcoder.net/n/e21f70424dca4acdacc1adf783c38ff7
看的时候又发现了点其他知识，就顺便查了查，hh

1.莫比乌斯反演
2.杜教筛
3.

思路：
因为k%i=k-floor(k/i) * i，所以转化为求n*k -sum(floor(n/i)*i).难点就是求 floor(n/i).因为对于某一区间的i, floor(n/i)
的值是一样的，所以采用除法分块的思想。

除法分块模板：

for(L=ll;L<=rr;L=R+1){
	R=min(rr,x/(x/L));
	ans+=(R-L+1)*(x/L);//x/L出现的次数为R-L+1次 
}

代码如下：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll ans,n,k;
int main(){
	cin>>n>>k;
	ans=n*k;
	for(ll L=1,R;L<=n;L=R+1){//每次都从R+1开始 
		R=k/L?min(k/(k/L),n):n;//当n<k时，分块右边的k/(k/L)可能会大于n,但是分块右边边界最大为n,所以用min(n,k/(k/L)). 
		ans-=(k/L)*(R-L+1)*(L+R)/2;//与模板略微不同的是，由于每一次后面乘了个i,所以这里需要将一个范围内的i累加，额外乘以一个(R+L)/2 
	}
	cout<<ans<<endl;
	return 0;
}

质数：5.Hankson的趣味题

思路：先预处理1-100000内的所有质数，再对每个质数p，计算出x可能包含多少个p,结果累乘。

代码如下：
(1)书上解法：

#include <bits/stdc++.h>
using namespace std;
const int N= 1e6+10;
typedef long long ll;
ll n,a0,a1,b0,b1,x;
ll prime[N],cnt,ans;
bool v[N];
void primes(ll n){//线性筛 
	for(ll i=2;i<=n;i++){
		if(!v[i]){
			prime[++cnt]=i;
		}
		for(ll j=1;i*prime[j]<=n;j++){
			v[i*prime[j]]=1;
			if(i%prime[j]==0) break;
		}
	}
}
ll get(ll &x,ll &y){//采用引用，计算时边分解边缩小d 
	ll num=0;
	while(x%y==0){
		num++;
		x/=y;
	}
	return num;
}
ll work(ll p){
	if(p>b1) return 0;
	if(b1%p!=0) return 1;
	ll cntp=0;
	ll numa0=get(a0,p);
	ll numa1=get(a1,p);
	ll numb0=get(b0,p);
	ll numb1=get(b1,p);
	if(numa0>numa1&&numb0<numb1&&numa1==numb1) cntp=1;
	else if(numa0>numa1&&numb0==numb1&&numa1<=numb1) cntp=1;
	else if(numa0==numa1&&numb0<numb1&&numa1<=numb1) cntp=1;
	else if(numa0==numa1&&numb0==numb1&&numa1<=numb1) cntp=numb1-numa1+1;
	else cntp=0;
	ans*=cntp;//累乘结果 
	return 1;
}
int main(){
	//预处理出100000范围内的质数 
	primes(100000);//鉴于在2e9范围内，那么只需要筛sqrt()即可，即100000 
	cin>>n;
	while(n--){
		cin>>a0>>a1>>b0>>b1;
		ans=1;
		for(ll j=1;j<=cnt;j++){//对于每个质因子p,计算x可能包含多少个p 
			ll p=prime[j];
			if(work(p)) continue;
			break;
		}
		if(b1>1) work(b1);//如果最后b1>1,表明b1是某个大质数的倍数或者本身就是质数，因此再对b1进行分解 
		cout<<ans<<endl;
	}
}

（2）推公式结论：https://blog.csdn.net/nuclearsubmarines/article/details/77603154

这个思路值得借鉴学习，但是确实时间复杂度还是挺高，在洛谷上能过，在AcWing 上会超时。

5.筛素数，欧拉函数：可见的点

考察：筛选素数（埃氏筛法/线性筛），欧拉函数

思路：
可知除了（1，0），（0，1），（1，1）这三颗钉子外，当且仅当1<=x,y<=N,x!=y且gcd(x,y)=1时，一个坐标为(x,y)的钉子才能被看到。
并且可知这些钉子关于直线y=x对称，所以可以只考虑其中的一半，即对于每个2<=y<=N,计算有多少个x满足1<=x<=y并且gcd(x,y)=1,x的数量就是phi(y);
则ans=3+2*sum(phi(y))

代码如下：
（1）线性筛+欧拉函数

#include <bits/stdc++.h>
using namespace std;
#define N 50000
int n,c,phi[N],prime[N],m;
bool v[N];
void euler(int n){//线性筛，求出2-n中每个数的欧拉函数 
	memset(v,0,sizeof(v));//最小质因子
	m=0;
	for(int i=2;i<=n;i++){
		if(v[i]==0){
			prime[m++]=i;
			phi[i]=i-1; 
		}
		for(int j=0;i*prime[j]<=n;j++){
			v[i*prime[j]]=1;
		    if(i%prime[j]==0){//体现线性 ，由性质4和性质5 若p|n且p^2|n，则phi(n)=phi(n/p)*p;若 p|n但p^2不能整除n，则phi(n)=phi(n/p)*(p-1)
		    	phi[i*prime[j]]=phi[i]*prime[j];
		    	break;
			}
			phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}
}

int main(){
	cin>>c;
	for(int i=1;i<=c;i++){
		cin>>n;
		euler(n);
		int ans=3;
		cout<<i<<" "<<n<<" ";
		for(int j=2;j<=n;j++){//对于每个2<=y<=N, 
			ans+=2*phi[j];
		}
		if(n==0) cout<<0<<endl;//如果初始是原点 
		else cout<<ans<<endl;
	}
	return 0;
	
}

（2）埃氏筛法+欧拉函数

#include <bits/stdc++.h>
using namespace std;
#define N 50000
int n,c,phi[N]; 
void euler(int n){//埃氏筛法，求出2-n中每个数的欧拉函数 
	for(int i=2;i<=n;i++) phi[i]=i;
	for(int i=2;i<=n;i++){
		if(phi[i]==i){
			for(int j=i;j<=n;j+=i){
				phi[j]=phi[j]/i*(i-1);
			}
		}
	}
}

int main(){
	cin>>c;
	for(int i=1;i<=c;i++){
		cin>>n;
		euler(n);
		int ans=3;
		cout<<i<<" "<<n<<" ";
		for(int j=2;j<=n;j++){
			ans+=2*phi[j];
		}
		if(n==0) cout<<0<<endl;
		else cout<<ans<<endl;
	}
	return 0;
	
}

6.欧拉函数，快速幂，龟速乘：最幸运的数字

思路：第一次遇到快速幂会爆的情况，www，这个题真的太棒了！！！强推！
在这里插入图片描述 快速幂（龟速乘强化版）

龟速乘就是求两个数的乘积的:ab%mod。
这个版本就是原来的ans=ansa%mod改成了ans=mul(ans,a,mod)%mod,
a=a*a%mod改成了a=mul(a,a,mod)%mod

ll mul(ll a,ll b,ll mod){//龟速乘 
	ll ans=0;
	while(b){
		if(b&1) ans=(ans+a)%mod;
		a=a*2%mod;
		b>>=1;
	}
	return ans;
}
ll qk(ll a,ll b,ll mod){//快速幂（龟速乘强化版） 
	ll ans=1;
	while(b){
		if(b&1) ans=mul(ans,a,mod)%mod;//here
		a=mul(a,a,mod)%mod;//here
		b>>=1;
	}
	return ans;
}

代码如下：


#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

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

ll mul(ll a,ll b,ll mod){//龟速乘 
	ll ans=0;
	while(b){
		if(b&1) ans=(ans+a)%mod;
		a=a*2%mod;
		b>>=1;
	}
	return ans;
}
ll qk(ll a,ll b,ll mod){//快速幂 （龟速乘强化版）
	ll ans=1;
	while(b){
		if(b&1) ans=mul(ans,a,mod)%mod;
		a=mul(a,a,mod)%mod;
		b>>=1;
	}
	return ans;
}

ll phi(ll n){//求欧拉函数 
	ll ans=n;
	for(ll i=2;i*i<=n;i++){
		if(n%i==0){
			ans=ans/i*(i-1);
			while(n%i==0) n/=i;
		}
	}
	if(n>1) ans=ans/n*(n-1);
	return ans;
}

int main(){
	ll L,num=0;
	while(cin>>L&&L){
		ll d=gcd(L,8);
		ll mod=L/d*9;
		ll p=phi(mod);//得到mod的欧拉函数 
		ll ans=1e18;
		for(ll x=1;x*x<=p;x++){//在1-sqrt(p)范围内枚举p的约数 
			if(p%x!=0) continue;//如果不是p的约数
			//如果x是p的约数， 10^x和1模mod同余，也就是10^x%mod=1 
			if(qk(10,x,mod)==1) ans=min(ans,x);
			if(qk(10,p/x,mod)==1) ans=min(ans,p/x);//根据约数成对存在，那么另一个约数就是p/x 
		}
		printf("Case %lld: ",++num);
		if(ans==1e18) puts("0");
		else printf("%lld\n",ans);
	}
	return 0;
}

先把A分解质因数，那么A^B的约数之和可以表示为
在这里插入图片描述
，每一项都是一个等比数列，用快速幂求出分子的值，将分母用费马小引理转化为乘法逆元。由于费马小引理的条件是p和mod互质，所以当p[i]-1是9901的倍数时，乘法逆元不存在，但是此时1+p+p^{2+```+p^(b*c[i])和1+1+1^2+```+1}(bc[i]) mod9901同余，那么第一项就等于bc[1]+1。对于每一项的值进行累乘，得到最终答案。

代码如下：


#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll a,b,m,p[20],c[20];
ll ans=1,mod=9901;
void divide(ll n){//求出n的所有约数以及约数个数 
	m=0;
	for(ll i=2;i*i<=n;i++){
		if(n%i==0){
			p[++m]=i,c[m]=0; //第m个约数对应的个数为c[m] 
			while(n%i==0){
				n/=i;
				c[m]++;
			}
		}
	}
	if(n>1) {
		p[++m]=n;
		c[m]=1;
	}
}
ll qk(ll a,ll b){//快速幂 
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}
int main(){
	cin>>a>>b;
	if(a==0) {
		puts("0");
		return 0;
	} 
	divide(a);//先求出a的所有约数 
	for(ll i=1;i<=m;i++){
		if((p[i]-1)%mod==0){//没有逆元时 
			ans=ans%mod*(b*c[i]+1)%mod;
			continue;
		}
		ll x=qk(p[i],b*c[i]+1)%mod;
		x=(x-1+mod)%mod;//分子 
		ll y=p[i]-1;
		y=qk(y,mod-2);//分母转化为乘法逆元 
		ans=ans*x%mod*y%mod;//费马小引理求乘法逆元 
	}
	cout<<ans<<endl;
	return 0;
}

8.同余：同余方程

思路：
ax≡1(mod b) 等价于 ax-1是b的倍数不妨设为-y倍，那么ax+by=1。
先用扩展欧几里得算法求出一组特解x0,y0,则x0就是原方程的一个解，通解为
所有模b与x0同余的整数，再通过取模操作把解的范围移动到1-b之间，就得到了最小整数解

代码如下：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll a,b,x,y;
ll exgcd(ll a,ll b,ll &x,ll &y){//用扩展欧几里得算法求出ax+by=1的一组解 
	if(!b){
		x=1,y=0;
		return a;
	}
	ll d=exgcd(b,a%b,x,y);
	ll z=x;x=y;y=z-y*(a/b);
	return d;
}
int main(){
	cin>>a>>b;
	exgcd(a,b,x,y);
	cout<<(x%b+b)%b<<endl;//通解为所有模b与x0同余的整数，再通过取模操作把解的范围移动到1-b之间，就得到了最小整数解 
	return 0;
}

(2)买不到的数目

考察：线性同余方程

思路：
x≡ma(mod b) 1<=m<=b-1;
即x=ma+n*b;
如果n>=0,那么x必然能用a,b表示出来，不合题意。
并且要x最大，所以n=-1。所以x=ma-b;
要使x最大，则因为 1<=m<=b-1 ，所以m=(b-1)

所以x=(b-1)*a-b=ab-a-b;

代码如下：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll a,b;

void rd(long long &x){
    char c;
    for (x=0,c=getchar();c<48;c=getchar());
    for (;c>47;c=getchar()) x=(x<<1)+(x<<3)+(c^48);
}



int main(){
	rd(a),rd(b);
	cout<<a*b-a-b<<endl;
	return 0;
}

9.中国剩余定理:表达整数的奇怪方式

考察：
中国剩余定理，扩展中国剩余定理

**思路：**啊啊啊，还不是很懂，懂了再来更吧！

代码如下：

//扩展中国剩余定理模板题 
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll R[50],M[50],n;
ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){
		x=1,y=0;
		return a;
	}
	ll d=exgcd(b,a%b,x,y);
	ll z=x;
	x=y;
	y=z-y*(a/b);
	return d;
} 
ll exgct(){
	ll m=M[1],r=R[1],x,y,gcd;
	for(ll i=2;i<=n;i++){
		gcd=exgcd(m,M[i],x,y);
		if((r-R[i])%gcd!=0) return -1;
		x=(r-R[i])/gcd*x%M[i];
		r-=m*x;
		m=m/gcd*M[i];
		r%=m;
	}
	return (r%m+m)%m;
}
int main(){
	cin>>n;
	for(ll i=1;i<=n;i++){
		cin>>M[i]>>R[i];
	}
	cout<<exgct()<<endl;
	return 0;
}

求组合数：10.计算系数

快速求组合数方法

考察：二项式，组合数

思路：
由二项式展开可知，x^n*ym的系数为c(k,n)*a^n*bn
用快速幂求出a^n和bn，利用lucas定理算出c(k,n)

代码如下：


#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll a,b,k,n,m,mod=10007;
ll qk(ll a,ll b,ll mod){//快速幂 
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
} 
ll comp(ll n,ll m,ll mod){//组合数 
	ll ans=1,up=1,down=1;
	if(n<m||m<0) return 0;
	if(n==m) return 1;
	for(ll i=0;i<m;i++){
		up=up*(n-i)%mod;
		down=down*(m-i)%mod;
	}
	ans=up*qk(down,mod-2,mod)%mod;
	return ans;
}
ll lucas(ll n,ll m,ll mod){//lucas定理 
	ll ans=1;
	while(n&&m&&ans){
		ans=ans*comp(n%mod,m%mod,mod)%mod;
		n/=mod;
		m/=mod;
	}
	return ans;
}
int main(){
	cin>>a>>b>>k>>n>>m;
	ll ans=1;
	ans*=qk(a,n,mod)%mod;
	ans=ans*qk(b,m,mod)%mod;
	ans=ans*lucas(k,n,mod)%mod;
	cout<<ans<<endl;
	return 0;
}

11.容斥原理：

1.三重集：
A∪B∪C=A+B+C - A∩B - A∩C - B∩C + A∩B∩C

2.四重集：
A∪B∪C∪D=A+B+C+D - A∩B - B∩C - C∩A - A∩D - B∩D - C∩D + A∩B∩C + A∩B∩D + A∩C∩D + B∩C∩D - A∩B∩C∩D （规律集合数 奇加偶减）

eg:给出一个数N，求1至N中，有多少个数不是2 3 5 7的倍数。例如N = 10，只有1不是2 3 5 7的倍数。

输入
输入1个数N(1 <= N <= 10^18)。
输出
输出不是2 3 5 7的倍数的数共有多少。
输入样例
10
输出样例
1

代码如下：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

void solve(){
	int n;
	cin>>n;
	ll ans=0;
	ans=n/2+n/3+n/5+n/7;
	ans=ans-n/6-n/10-n/14-n/15-n/21-n/35+n/30+n/42+n/105+n/70-n/210;
	cout<<n-ans<<endl;
	
}


int main(){
	int t;
	t=1;
	while(t--){
		solve();
	}
	return 0;
}

3.容斥原理：
容斥原理1

在这里插入图片描述容斥原理2

(2)硬币购物

思路：
反面思考背包+容斥 +差分
容斥：奇加偶减，反面思考，先预处理，用完全背包求出没有限制的方案数
总答案为：没有限制的方案数-有一个限制的方案数+ 有两个限制的方案数-有三个限制的方案数+有四个限制的方案数

考察：
背包+容斥 +差分

代码如下：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5+10;
int n;
ll dp[N],c[5],d[5],ans;

//dfs 
void dfs(int index,ll sum,int f){//f的初值为1 
	if(sum<0) return;
	if(index>4){
		ans+=f*dp[sum];
		return;
	}
	dfs(index+1,sum,f);//合法 
	//不合法 
	dfs(index+1,sum-c[index]*(d[index]+1),-f);//奇加偶减 
}
int main(){
	for(int i=1;i<=4;i++) cin>>c[i];
	cin>>n;
	dp[0]=1;
	//预处理：完全背包求出1-N之间的无限制方案数 
	for(int i=1;i<=4;i++){
		for(int j=c[i];j<N;j++){
			dp[j]+=dp[j-c[i]];
		}
	}
	while(n--){
		for(int i=1;i<=4;i++) cin>>d[i];
		ans=0;//每次记得清0 
		ll s;
		cin>>s;
		dfs(1,s,1);//进行dfs 
		cout<<ans<<endl;
	}
}

12.高斯消元

传送门
例题：
1.球形空间产生器

思路：一个球体上的所有点到球心的距离相等，所以只需要求出一个点(x1,x2,x3,x4```xn)，使得sum(ai,j-xj)^2=C,C为常数，通过相邻两方程作差，可以得到：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 20;
int n;
double a[N][N],b[N],c[N][N];

int main(){
	cin>>n;
	for(int i=1;i<=n+1;i++){
		for(int j=1;j<=n;j++){
			cin>>a[i][j];
		}
	}
	//c是系数矩阵，b是常数，二者一起构成增广矩阵 
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			c[i][j]=2*(a[i][j]-a[i+1][j]);//系数 
			b[i]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];//常数 
	
		}
	}
	//高斯消元 
	for(int i=1;i<=n;i++){
		for(int j=i;j<=n;j++){
			if(fabs(c[j][i])>1e-8){//先找到系数不为0的一个方程 
				for(int k=1;k<=n;k++) swap(c[i][k],c[j][k]);
				swap(b[i],b[j]);
			}
		}
		//消去其他方程的x[i]的系数 
		for(int j=1;j<=n;j++){
			if(i==j) continue;
			double rate = c[j][i]/c[i][i];//计算出倍数 
			for(int k=i;k<=n;k++) c[j][k]-=c[i][k]*rate;//消去 
			b[j]-=b[i]*rate;
		}
	}
	for(int i=1;i<n;i++) printf("%.3lf ",b[i]/c[i][i]);
	printf("%.3lf\n",b[n]/c[n][n]);
	return 0;
}

Bold!

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
数学知识啦啦啦

数学知识模板：https://www.acwing.com/blog/content/406/Practice:1.质数距离https://www.acwing.com/problem/content/description/198/**考察：**线性筛+巧妙优化思路：因为l,r<=2^31,所以不能直接筛出1-r的所有素数。因为任何一个合数n必定包含一个不超过sqrt(n)的质因子，又sqrt(2^31)=50000,所以先预处理，用筛法筛出2-50000之间的所有质数。因为质数的倍
复制链接

扫一扫