[jzoj]3498. 【NOIP2013模拟联考14】图形变换(transform) （计算几何+矩阵乘法）

Problem

四种操作：

Trans(dx,dy) 表示平移图形，即把图形上所有的点的横纵坐标分别加上dx和dy；

Scale(sx,sy) 表示缩放图形，即把图形上所有点的横纵坐标分别乘以sx和sy；

Rotate(θ,x0,y0) 表示旋转图形，即把图形上所有点的坐标绕(x0,y0)顺时针旋转θ角度

Loop(m)

…

End

表示把Loop和对应End之间的操作循环执行m次，循环可以嵌套。

Data constraint

保证操作中坐标值不会超过double范围，输出不会超过int范围；

指令总共不超过1000行；

对于所有的数据,所有循环指令中m<=1000000；

对于60%的数据,所有循环指令中m<=1000；

对于30%的数据不含嵌套循环。

Solution

矩阵乘法

• 矩阵乘法有许多优美的性质，如：

  • 结合律.（非常重要，这也是矩阵乘法的用处）
  • 分配率.
  • 却不满足分配率
    
    • 如：
    • 

• 矩阵A乘矩阵B，那么矩阵A的列数必须等于矩阵B的行数.

• 如果矩阵A是m∗n，矩阵B是n∗r，那么乘出来的矩阵C应该是m∗r的.

• 结合律：

  • 因为满足结合律，所以我们还可以再求n个相同矩阵相乘时运用快速幂.

回到正题

• 为什么可以用矩阵乘法？

• 因为此题其实只要两种操作，加法，乘法，而这样的操作却相同的重复了很多次.

• 于是，我们可以把一个循环看作一个矩阵，也就是在一个循环中，每一次循环，答案矩阵都会乘上一个对应矩阵.

• 这个对应矩阵是相同的，所以我们可以运用结合律，来进行快速幂，使时间复杂度降为O(log2n)

• 因为题目中的旋转，是相当于乘法+加法，所以我们同样可以处理.

• 我们可以设一个3∗3的矩阵.

• f[1][1/2/3]表示横坐标，加了多少，乘了多少个x[i]，乘了多少个y[i].

• f[2][1/2/3]表示纵坐标，加了多少，乘了多少个x[i]，乘了多少个y[i].

• 最后我们再x[i],y[i]带进去，就可以求得答案了.

• 我们来看一下3*3的矩阵，是如何进行乘法的.

• code
  

#include <iostream>#include <cstdio>#include <cstring>#include <cstdlib>#include <cmath>#define PI M_PIconst int Maxn = 110, Maxm = 1010;using namespace std;typedef double matrix[4][4];double x[Maxn],y[Maxn];char ch;int n,len;matrix g,h;struct set{    int bz,time;    double x,y,z;} a[Maxm];void multip(matrix &a, matrix b){    matrix c;    memset(c,0,sizeof(c));    for (int i=1;i<=3;i++)        for (int j=1;j<=3;j++)            for (int k=1;k<=3;k++)                c[i][j] += a[i][k] * b[k][j];    memcpy(a,c,sizeof(c));}int dfs(int time,int first){    int i = first-1;    matrix f;    memset(f,0,sizeof(f)), f[1][1] = f[2][2] = f[3][3] = 1;    while (++i<=len)    {        memset(g,0,sizeof(g)), g[1][1] = g[2][2] = g[3][3] = 1;        if (a[i].bz == 5) break;        if (a[i].bz == 1) g[3][1] = a[i].x, g[3][2] = a[i].y; else        if (a[i].bz == 2) g[1][1] = a[i].x, g[2][2] = a[i].y; else        if (a[i].bz == 3)        {            double p=((360-a[i].z)*PI)/180, w1 = cos(p), w2 = sin(p);            g[1][1] = w1, g[1][2] = w2, g[2][1] =-w2, g[2][2] = w1;            g[3][1] =(-a[i].x)*w1+a[i].y*w2+a[i].x, g[3][2] = (-a[i].y)*w1-a[i].x*w2+a[i].y;        }; //注意一下这里的旋转操作，必须先转到与原点相关，然后再放回去.        if (a[i].bz != 4) multip(f,g); else        {            int k = a[i].time;            i = dfs(a[i].time,i+1);            while (k)            {                if (k & 1) multip(f,h);                multip(h,h), k >>= 1;            }        }    }    memcpy(h,f,sizeof(f));    return i;}int main(){    freopen("transform.in","r",stdin);    freopen("transform.out","w",stdout);    scanf("%d",&n);    for (int i=1;i<=n;i++)        scanf("%lf%lf\n",&x[i],&y[i]);    while ( scanf("\n%c",&ch) != EOF)    {        if (ch=='T') a[++len].bz = 1, scanf("rans(%lf,%lf)", &a[len].x,&a[len].y); else        if (ch=='S') a[++len].bz = 2, scanf("cale(%lf,%lf)", &a[len].x,&a[len].y); else        if (ch=='R') a[++len].bz = 3, scanf("otate(%lf,%lf,%lf)",&a[len].z,&a[len].x,&a[len].y); else        if (ch=='L') a[++len].bz = 4, scanf("oop(%d)",&a[len].time); else                     a[++len].bz = 5, scanf("nd");    }    dfs(1,1);    for (int i=1;i<=n;i++)    {        matrix f;        memset(f,0,sizeof(f)), f[1][1] = x[i], f[1][2] = y[i], f[1][3] = 1;        multip(f,h);        printf("%.4lf %.4lf\n",f[1][1],f[1][2]);    }}