bupt244

source: http://acm.bupt.edu.cn/onlinejudge/newoj/showProblem/show_problem.php?problem_id=244

title:     A Letter  to Programmers 

 

/*
坐标变换的公式+矩阵快速幂  
*/
#include <iostream>
#include <cmath>
using namespace std;
const int N = 105;
const double pi = acos(-1.0);
const double eps = 1e-4;
struct point{
	double x, y, z;
	point(double _x=0, double _y=0, double _z=0):x(_x), y(_y), z(_z){};
	void read(){
		scanf("%lf%lf%lf", &x, &y, &z);
	}
};
//三维点的变换矩阵,点做变换的时候是乘在矩阵的右边的，所以，后来的变换矩阵应该乘在
//原有矩阵的左边。
struct  mat{
	static const int N = 4;
	double v[N][N];
	int n;
	mat(){  n =  N; }
	void  setv(double _v[N][N]){
		int i, j;
		for(i = 0;  i < N; i++){
			for(j = 0; j < N; j++){
				v[i][j] = _v[i][j];
			}
		}
	}
	//点绕轴(0,0,0)->(x,y,z)旋转(旋转方向为从 点(x,y,z)往(0,0,0))
	//看是逆时针方向.如果所绕的轴的一个端点不在原点，则可以移到端
	//点再移回来
	void rotate(double x, double y, double z, double ang){
		double cosa = cos(ang), sina = sin(ang), len;
		len = sqrt(x*x+y*y+z*z);
		x/=len; y/=len; z/=len;
		double a[N][N] = {
			cosa+(1-cosa)*x*x, (1-cosa)*x*y-z*sina, (1-cosa)*x*z+y*sina, 0,
			(1-cosa)*x*y+z*sina, cosa+(1-cosa)*y*y, (1-cosa)*y*z-x*sina, 0,
			(1-cosa)*x*z-y*sina, (1-cosa)*y*z+x*sina, cosa+(1-cosa)*z*z, 0,
			0, 0, 0, 1
		};
		setv(a);
	}
	//移动矩阵，点(x0,y0,z0)变为(x0+x,y0+y,z0+z)
	void translate(double x, double y, double  z){
		double a[N][N] = {
			1, 0, 0, x,
			0, 1, 0, y,
			0, 0, 1, z,
			0, 0, 0, 1
		};
		setv(a);
	}
	//伸缩矩阵,点(x0,y0,z0)变为(x0*x,y0*y,z0*y)
	void scale(double x, double y, double z){
		double a[N][N] = {
			x, 0, 0, 0,
			0, y, 0, 0,
			0, 0, z, 0,
			0, 0, 0, 1
		};
		setv(a);
	}
	//对点st应用该矩阵
	point change(point st){
		double ds[N] = {st.x, st.y, st.z, 1}, ans[N-1];
		int i, j;
		//点做矩阵变换的时候是点是乘在矩阵的右边的
		for(i = 0; i < N-1; i++){
			ans[i] = 0;
			for(j = 0; j < N; j++){
				ans[i] += v[i][j] * ds[j];
			}
		}
		return point(ans[0], ans[1], ans[2]);
	}

	//矩阵乘法
	mat operator*(mat tm){
		mat ans;
		int i, j, k;
		for(i = 0; i < n; i++){
			for(j = 0; j < n; j++){
				ans.v[i][j] = 0;
				for(k = 0; k < n; k++){
					ans.v[i][j] += v[i][k] * tm.v[k][j];
				}
			}
		}
		return ans;
	}
	//矩阵幂
	mat operator^(int k){
		mat ans, org = *this;
		ans.scale(1., 1., 1.); //初始化为单位矩阵
		for(; k; k >>= 1, org = org * org){
			if(k&1){
				ans = ans * org;
//				ans = org * ans;
			}
		}
		return ans;
	}
};
char ts[N];  //ts[i]为'm'表示矩阵，为'r'表示重复
mat ms[N];
int  ks[N];

int main() {
//	freopen("in.txt", "r", stdin);

	int n, top, i, k;
	char op[20];
	double tx, ty, tz, d;
	mat m;
	point tp;
	while(~scanf("%d", &n) && n > 0){
		top=-1;
		while(true){
			scanf("%s", op);
			if(strcmp(op, "translate") == 0){
				scanf("%lf%lf%lf", &tx, &ty, &tz);
				ms[++top].translate(tx, ty, tz);
				ts[top] = 'm';
			}else if(strcmp(op, "scale") == 0){
				scanf("%lf%lf%lf", &tx, &ty, &tz);
				ms[++top].scale(tx, ty, tz);
				ts[top] = 'm';
			}else if(strcmp(op, "rotate") == 0){
				scanf("%lf%lf%lf%lf", &tx, &ty, &tz, &d);
				d = (d * pi) / 180.0;
				ms[++top].rotate(tx, ty, tz, d);
				ts[top] = 'm';
			}else if(strcmp(op, "repeat") == 0){
				scanf("%d", &k);
				ks[++top] = k;
				ts[top] = 'r';
			}else if(strcmp(op, "end") == 0){
				m.scale(1, 1, 1);
				while(top >= 0 && ts[top] != 'r'){
					m = m * ms[top];
					top--;
				}
				if(top < 0) break;
				m = m ^ ks[top];
				ms[top] = m;
				ts[top] = 'm';
			}
		}
		for(i = 0; i < n; i++){
			tp.read();
			tp = m.change(tp);
			printf("%.2f %.2f %.2f\n", tp.x + eps, tp.y + eps, tp.z + eps);
		}
		printf("\n");
	}
	return 0;
}