[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 
#include 
#include 
#include 
#include 

#define PI M_PI

const 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]);
    }
}