洛谷P3945 三体问题

题目背景

三体人所居住的星系由于三体运动的不确定性而导致三体星人生活动荡不安，善良的人类程序员（也就是你了！伟大的英雄！）决定帮助愚蠢得连程序都不会写的三体星人模拟天体的运动轨迹。这时，无聊的“歌者”文明决定戏弄一下你，于是给三体星系添加了一些新的星体。

题目描述

输入N个天体与他们在空间中的坐标（xi,yi,zi）、初速度（vx,vy,vz）与质量Mi，已知三体世界受到“歌者”影响时间的流动不是连续的（每0.01秒钟刷新一次），天体均视为质点，求t时刻所有天体的坐标。

#本题G取(6.67408*10^(-11))，在代码中可以写成：

#define G 6.67408e-11
###当你的答案与标准答案的相对误差不超过0.5%的时候，你在本测试点得到AC。

###也就是说，保留多少位小数你可以自行确定。

###标准答案将会保留12位小数。本题开启SPJ判断你的答案是否正确。

输入格式

第一行输入天体数N与时刻t，接下来逐行输入以空格分隔的各天体坐标、质量与初速度（三个方向上的分速度）。

输出格式

N行，每行为第i个天体在t时刻的坐标xi,yi,zi，空格隔开。\

题解

三维矢量

struct vector
{
	long double x,y,z;
	vector(long double xx=0,long double yy=0,long double zz=0)
	{
		x=xx,y=yy,z=zz;
	}
	void clear()
	{
		x=y=z=0;
	}
	void in()
	{
		cin>>x>>y>>z;
	}
	void out()
	{
		printf("%Lf %Lf %Lf",x,y,z);
	}
	vector operator+(const vector& b)const
	{
		return vector(x+b.x,y+b.y,z+b.z);
	}
	vector operator+=(const vector& b)
	{
		x+=b.x;
		y+=b.y;
		z+=b.z;
		return *this;
	}
	vector operator*(const long double& b)
	{
		return vector(x*b,y*b,z*b);
	}
	vector operator*=(const long double& b)
	{
		x*=b;
		y*=b;
		z*=b;
		return *this;
	}
	vector operator/(const long double& b)const
	{
		return vector(x/b,y/b,z/b);
	}
}F[100];

星球

struct Star
{
	vector pos,v;
	long double m;
	IL void in()
	{
		pos.in();
		cin>>m;
		v.in();
	}
	IL void out()
	{
		pos.out();
	}
	
}star[100];

计算

for(;T>0;T-=dt)
{
	for(int i=1;i<=n;i++) F[i].clear();
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			if(i==j) continue;
			r=(star[i].pos.x-star[j].pos.x)*(star[i].pos.x-star[j].pos.x)+(star[i].pos.y-star[j].pos.y)*(star[i].pos.y-star[j].pos.y)+(star[i].pos.z-star[j].pos.z)*(star[i].pos.z-star[j].pos.z);
			f=G*star[j].m/r;
			F[i]+=vector(f*(star[j].pos.x-star[i].pos.x)/sqrt(r),f*(star[j].pos.y-star[i].pos.y)/sqrt(r),f*(star[j].pos.z-star[i].pos.z)/sqrt(r));
		}
	}
	for(int i=1;i<=n;i++)
	{
		star[i].v+=F[i]*dt;
		star[i].pos+=(star[i].v)*dt;
	}
}

完整代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<limits.h>
#define IL inline
#define re register
#define LL long long
#define ULL unsigned long long
#ifdef TH
#define debug printf("Now is %d\n",__LINE__);
#else
#define debug
#endif
using namespace std;
const long double G = 6.67408e-11;
const long double dt = 1e-2;
int n;
long double T;
struct vector
{
	long double x,y,z;
	vector(long double xx=0,long double yy=0,long double zz=0)
	{
		x=xx,y=yy,z=zz;
	}
	void clear()
	{
		x=y=z=0;
	}
	void in()
	{
		cin>>x>>y>>z;
	}
	void out()
	{
		printf("%Lf %Lf %Lf",x,y,z);
	}
	vector operator+(const vector& b)const
	{
		return vector(x+b.x,y+b.y,z+b.z);
	}
	vector operator+=(const vector& b)
	{
		x+=b.x;
		y+=b.y;
		z+=b.z;
		return *this;
	}
	vector operator*(const long double& b)
	{
		return vector(x*b,y*b,z*b);
	}
	vector operator*=(const long double& b)
	{
		x*=b;
		y*=b;
		z*=b;
		return *this;
	}
	vector operator/(const long double& b)const
	{
		return vector(x/b,y/b,z/b);
	}
}F[100];
struct Star
{
	vector pos,v;
	long double m;
	IL void in()
	{
		pos.in();
		cin>>m;
		v.in();
	}
	IL void out()
	{
		pos.out();
	}
	
}star[100];
int main()
{
	cin>>n>>T;
	for(int i=1;i<=n;i++) star[i].in();
	long double f,r;
	for(;T>0;T-=dt)
	{
		for(int i=1;i<=n;i++) F[i].clear();
		for(int i=1;i<=n;i++)
		{
			for(int j=1;j<=n;j++)
			{
				if(i==j) continue;
				r=(star[i].pos.x-star[j].pos.x)*(star[i].pos.x-star[j].pos.x)+(star[i].pos.y-star[j].pos.y)*(star[i].pos.y-star[j].pos.y)+(star[i].pos.z-star[j].pos.z)*(star[i].pos.z-star[j].pos.z);
				f=G*star[j].m/r;
				F[i]+=vector(f*(star[j].pos.x-star[i].pos.x)/sqrt(r),f*(star[j].pos.y-star[i].pos.y)/sqrt(r),f*(star[j].pos.z-star[i].pos.z)/sqrt(r));
			}
		}
		for(int i=1;i<=n;i++)
		{
			star[i].v+=F[i]*dt;
			star[i].pos+=(star[i].v)*dt;
		}
	}
	for(int i=1;i<=n;i++)
	{
		star[i].out();
		cout<<endl;
	}
	return 0;
}

注意第9595行，r的定义是r^2 ，不先开方是为了精度。
以及玄学的物理计算在第102,103行