<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

265. Wizards
time limit per test: 0.5 sec.
memory limit per test: 65536 KB
input: standard
output: standard



The magical school of nature exists for thousands years. The wizards of the school are mastered in controlling animate nature. But it is still extremely difficult for them to manipulate matter itself. The best mages of the school were trying to learn the ways of manipulating the world around them for a long time. They started with the space distortion. Ancient books mentioned three basic space distortion spells: Translation, Scaling and Rotation. But all the experiments failed because it was extremely hard to predict the way the spell would affect the space. Nowadays the experiments arose again and you are to write the program which would help wizards in predicting of the results of some complex distortion spells. Complex spell consists of N consecutive basic distortion spells. And you are to predict its effect on M control vertices. For your convenience wizards measure the space in some 3-dimensional Cartesian system.
While casting Translation spell wizard specifies Translation vector and all the vertices are moved by the given x, y and z values.
While casting Scaling spell wizard specifies three Scaling parameters which show the Scaling factor for each of three coordinates (each x, y, and z coordinate of every vertex is multiplied by the corresponding Scaling factor).
While casting Rotation spell wizard specifies a point and Rotation angle and all the vertices are rotated in a clockwise direction about the ray from the origin through the specified point. The angle parameter specifies the angle of Rotation in degrees.

Input
The first line of the input file contains integer number N (1 <= N <= 1000). The next N lines describe complex spell. Each line contains one basic spell description in order they are applied. The first letter of each line specifies basic spell type. `T' stands for Translation spell, `S' stands for Scaling spell and `R' stands for Rotation spell. For Translation spell Translation vector is specified on the same line (three numbers: x, y, and z). For Scaling spell three Scaling factors are specified on the same line (three numbers: x, y, and z, each not greater than 1 by absolute value). For Rotation spell the point and the angle are specified on the same line (x, y, z and a, x^2+y^2+z^2>0). After that vertices description follows. First an integer number of vertices M is specified (1 <= M <= 100000). Then descriptions of M vertices follow. Each vertex description consists of 3 numbers - Cartesian coordinates of the vertex.
All the numbers in input are real and do not exceed 1000 by their absolute value.

Output
For each vertex in the input output its coordinates after the entire spell takes effect.

Sample test(s)

Input
Test #1
3
T 1.0 -2.4 0.7
S 0.2 4.0 2.0
R 0.0 0.0 1.0 2.0
3
0.0 0.0 1.0
0.0 1.0 0.0
1.0 0.0 0.0

Test #2
3
S 1.0 2.0 3.0
R 1.0 0.0 0.0 90.0
R 0.0 1.0 0.0 90.0
3
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0

Output
Test #1
0.53 -9.59 3.40
0.40 -5.59 1.40
0.73 -9.58 1.40

Test #2
0.00 0.00 -1.00
2.00 0.00 0.00
0.00 -3.00 0.00
Author:	Igor A. Kulkin
Resource:	Saratov SU Contest: Golden Fall 2004
Date:	October 2, 2004

<|response|>
1. Abridged Problem Statement

We are given N 3D transformations applied in order:

- T x y z: translate every point by vector (x, y, z).
- S x y z: scale coordinates independently by factors (x, y, z).
- R x y z a: rotate clockwise by a degrees around the ray from the origin through point (x, y, z).

After these transformations, we are given M vertices. For each vertex, output its final coordinates.

Constraints:
- 1 <= N <= 1000
- 1 <= M <= 100000

A direct simulation of every transformation on every point would be too slow: O(NM).

2. Key Observations

Observation 1: Every operation is affine
Translation, scaling, and rotation in 3D can all be represented as affine transformations. An affine transformation has the form p' = A * p + b where A is a 3 x 3 matrix and b is a translation vector.

Observation 2: Use homogeneous coordinates
To combine all transformations uniformly, represent a point (x, y, z) as the column vector (x, y, z, 1). Then every transformation becomes a 4 x 4 matrix, and applying it is p' = M * p.

Observation 3: Compose all transformations first
If transformations are applied in this order A1, A2, A3, ..., AN, then the final matrix is Final = AN * ... * A3 * A2 * A1. The newest transformation is multiplied on the left. So while reading transformations, if cur is the current transformation matrix, Final = cur * Final. Then each point is transformed once.

Observation 4: Rotation uses Rodrigues' formula
Rotation is around an arbitrary axis from the origin through (x, y, z). First normalize the axis: len = sqrt(x^2 + y^2 + z^2), ux = x/len, uy = y/len, uz = z/len. Let c = cos(angle), s = sin(angle). Then the rotation matrix is:
  [ c+ux²(1-c)      uxuy(1-c)-uzs   uxuz(1-c)+uys ]
  [ uyux(1-c)+uzs   c+uy²(1-c)      uyuz(1-c)-uxs ]
  [ uzux(1-c)-uys   uzuy(1-c)+uxs   c+uz²(1-c)    ]
This is placed in the upper-left 3 x 3 block of the 4 x 4 matrix.

3. Full Solution Approach

Step 1: Initialize the final matrix
Start with the identity matrix (which represents doing nothing).

Step 2: Build a matrix for each spell
Translation by (tx, ty, tz):
  [1 0 0 tx]
  [0 1 0 ty]
  [0 0 1 tz]
  [0 0 0  1]
Scaling by (sx, sy, sz):
  [sx 0  0  0]
  [0  sy 0  0]
  [0  0  sz 0]
  [0  0  0  1]
Rotation around axis (x, y, z) by angle a degrees: normalize (x, y, z), convert angle to radians, apply Rodrigues' formula to fill the upper-left 3 x 3 block. The last row and column stay as in the identity matrix.

Step 3: Compose transformations
If Final stores all previous transformations, and cur is the current spell matrix, then Final = cur * Final. This order is important because points are treated as column vectors.

Step 4: Apply final matrix to each point
For each point (x, y, z):
  x' = Final[0][0]*x + Final[0][1]*y + Final[0][2]*z + Final[0][3]
  y' = Final[1][0]*x + Final[1][1]*y + Final[1][2]*z + Final[1][3]
  z' = Final[2][0]*x + Final[2][1]*y + Final[2][2]*z + Final[2][3]
Print the result with two digits after the decimal point.

Complexity
Matrix size is constant, so each matrix multiplication costs O(1). Total complexity O(N + M); memory O(1) apart from input/output buffers.

4. C++ Implementation

```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

const double PI = acos(-1.0);

namespace fastio {
constexpr int BUF_SIZE = 1 << 16;
char buf[BUF_SIZE];
int buf_pos = 0, buf_len = 0;

int get_char() {
    if(buf_pos == buf_len) {
        buf_len = (int)fread(buf, 1, BUF_SIZE, stdin);
        buf_pos = 0;
    }
    return buf_pos < buf_len ? (unsigned char)buf[buf_pos++] : EOF;
}

int skip_ws() {
    int c = get_char();
    while(c == ' ' || c == '\n' || c == '\r' || c == '\t') {
        c = get_char();
    }
    return c;
}

int read_int() {
    int c = skip_ws();
    int sign = 1;
    if(c == '-') {
        sign = -1;
        c = get_char();
    } else if(c == '+') {
        c = get_char();
    }
    int x = 0;
    while(c >= '0' && c <= '9') {
        x = x * 10 + (c - '0');
        c = get_char();
    }
    return sign * x;
}

double read_double() {
    int c = skip_ws();
    int sign = 1;
    if(c == '-') {
        sign = -1;
        c = get_char();
    } else if(c == '+') {
        c = get_char();
    }
    double x = 0.0;
    while(c >= '0' && c <= '9') {
        x = x * 10.0 + (c - '0');
        c = get_char();
    }
    if(c == '.') {
        c = get_char();
        double f = 0.1;
        while(c >= '0' && c <= '9') {
            x += (c - '0') * f;
            f *= 0.1;
            c = get_char();
        }
    }
    return sign * x;
}

char read_char() { return (char)skip_ws(); }
}  // namespace fastio

int n, m;
vector<tuple<char, array<double, 4>>> spells;
vector<array<double, 3>> points;

void read() {
    n = fastio::read_int();
    spells.assign(n, {});
    for(auto& [c, params]: spells) {
        c = fastio::read_char();
        int k = (c == 'R') ? 4 : 3;
        for(int i = 0; i < k; i++) {
            params[i] = fastio::read_double();
        }
    }

    m = fastio::read_int();
    points.assign(m, {});
    for(auto& p: points) {
        p[0] = fastio::read_double();
        p[1] = fastio::read_double();
        p[2] = fastio::read_double();
    }
}

void solve() {
    // Each basic spell is an affine map on R^3, so on the homogeneous vector
    // (x, y, z, 1) it becomes a plain linear map - a 4x4 matrix. Concretely:
    //
    //   Translation(tx,ty,tz): identity 4x4 with (tx,ty,tz) in the last column.
    //   Scaling(sx,sy,sz):     diag(sx, sy, sz, 1).
    //   Rotation(axis, angle): Rodrigues' formula in the upper-left 3x3, with
    //                          the last row and column kept as (0,0,0,1).
    //
    // Composing the whole spell list is one big 4x4 product M = T_n * ... * T_1
    // (each new spell multiplies on the left because it is applied after all
    // earlier ones), and producing the final coordinates for every vertex is
    // just M * (x, y, z, 1)^T. So we pay O(N) for the composition and O(M)
    // for the application, independently of each other.
    //
    // Translation and scaling are immediate from the definitions above - a
    // translation adds the vector to (x,y,z) and a scaling multiplies each
    // coordinate by its factor, which is exactly what the matrices encode.
    // Rotation is the interesting one because the axis is arbitrary, not one
    // of the coordinate axes. We first normalize the axis to a unit vector
    // u = (ux, uy, uz) (the input only guarantees x^2 + y^2 + z^2 > 0). Then
    // Rodrigues' formula gives the 3x3 rotation matrix R (what we drop into
    // the upper-left of t) by angle a around u as
    //
    //   R = I + sin(a) * K + (1 - cos(a)) * K^2,
    //
    // where I is the 3x3 identity and K is the cross-product matrix of u
    // (so K*v = u x v for every vector v):
    //
    //          [  0  -uz   uy ]
    //      K = [ uz    0  -ux ].
    //          [-uy   ux    0 ]
    //
    // Geometrically: any vector v splits into a part parallel to u (which is
    // fixed by the rotation) and a part perpendicular to u (which lives in
    // the plane orthogonal to u and is rotated there by angle a using the
    // 2D rotation written in the orthonormal pair (v_perp, u x v_perp));
    // expanding that gives exactly the R above. We hardcode the resulting
    // nine entries of R directly into the 3x3 block of t below to avoid
    // building K and K^2 explicitly.
    //
    // The problem says rotation is "clockwise about the ray from origin
    // through the point". Standard Rodrigues with a right-handed axis already
    // matches this when you look along the ray from the origin outward (it
    // matches the sample with R about +y by 90 sending (1,0,0) to (0,0,-1)).

    auto id = []() {
        array<array<double, 4>, 4> r{};
        for(int i = 0; i < 4; i++) {
            r[i][i] = 1.0;
        }
        return r;
    };

    auto mul = [](const array<array<double, 4>, 4>& a,
                  const array<array<double, 4>, 4>& b) {
        array<array<double, 4>, 4> r{};
        for(int i = 0; i < 4; i++) {
            for(int j = 0; j < 4; j++) {
                double s = 0.0;
                for(int k = 0; k < 4; k++) {
                    s += a[i][k] * b[k][j];
                }
                r[i][j] = s;
            }
        }
        return r;
    };

    auto M = id();
    for(auto& [c, p]: spells) {
        auto t = id();
        if(c == 'T') {
            t[0][3] = p[0];
            t[1][3] = p[1];
            t[2][3] = p[2];
        } else if(c == 'S') {
            t[0][0] = p[0];
            t[1][1] = p[1];
            t[2][2] = p[2];
        } else {
            double x = p[0], y = p[1], z = p[2];
            double len = sqrt(x * x + y * y + z * z);
            x /= len;
            y /= len;
            z /= len;

            double rad = p[3] * PI / 180.0;
            double co = cos(rad), si = sin(rad);
            t[0][0] = co + x * x * (1 - co);
            t[0][1] = x * y * (1 - co) - z * si;
            t[0][2] = x * z * (1 - co) + y * si;
            t[1][0] = y * x * (1 - co) + z * si;
            t[1][1] = co + y * y * (1 - co);
            t[1][2] = y * z * (1 - co) - x * si;
            t[2][0] = z * x * (1 - co) - y * si;
            t[2][1] = z * y * (1 - co) + x * si;
            t[2][2] = co + z * z * (1 - co);
        }

        M = mul(t, M);
    }

    cout << fixed << setprecision(2);
    for(auto& v: points) {
        double x = M[0][0] * v[0] + M[0][1] * v[1] + M[0][2] * v[2] + M[0][3];
        double y = M[1][0] * v[0] + M[1][1] * v[1] + M[1][2] * v[2] + M[1][3];
        double z = M[2][0] * v[0] + M[2][1] * v[1] + M[2][2] * v[2] + M[2][3];
        cout << x << ' ' << y << ' ' << z << '\n';
    }
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        solve();
    }

    return 0;
}
```

5. Python Implementation

```python
import sys
import math


def identity_matrix():
    """
    Return a 4x4 identity matrix.

    The identity matrix represents a transformation that does nothing.
    """
    mat = [[0.0] * 4 for _ in range(4)]

    for i in range(4):
        mat[i][i] = 1.0

    return mat


def multiply(a, b):
    """
    Multiply two 4x4 matrices.

    Returns:
        a * b
    """
    res = [[0.0] * 4 for _ in range(4)]

    for i in range(4):
        for j in range(4):
            total = 0.0

            for k in range(4):
                total += a[i][k] * b[k][j]

            res[i][j] = total

    return res


def fix_zero(x):
    """
    Avoid printing -0.00 caused by floating-point round-off errors.
    """
    if abs(x) < 0.0005:
        return 0.0
    return x


def main():
    data = sys.stdin.buffer.read().split()
    idx = 0

    n = int(data[idx])
    idx += 1

    final_matrix = identity_matrix()

    for _ in range(n):
        spell_type = data[idx].decode()
        idx += 1

        cur = identity_matrix()

        if spell_type == 'T':
            """
            Translation by vector (tx, ty, tz).

            Matrix:

            [1 0 0 tx]
            [0 1 0 ty]
            [0 0 1 tz]
            [0 0 0  1]
            """
            tx = float(data[idx])
            ty = float(data[idx + 1])
            tz = float(data[idx + 2])
            idx += 3

            cur[0][3] = tx
            cur[1][3] = ty
            cur[2][3] = tz

        elif spell_type == 'S':
            """
            Scaling by factors (sx, sy, sz).

            Matrix:

            [sx 0  0  0]
            [0  sy 0  0]
            [0  0  sz 0]
            [0  0  0  1]
            """
            sx = float(data[idx])
            sy = float(data[idx + 1])
            sz = float(data[idx + 2])
            idx += 3

            cur[0][0] = sx
            cur[1][1] = sy
            cur[2][2] = sz

        else:
            """
            Rotation around the ray from the origin through (x, y, z).

            We normalize the axis and use Rodrigues' rotation formula.
            """
            x = float(data[idx])
            y = float(data[idx + 1])
            z = float(data[idx + 2])
            angle_degrees = float(data[idx + 3])
            idx += 4

            length = math.sqrt(x * x + y * y + z * z)

            ux = x / length
            uy = y / length
            uz = z / length

            angle = angle_degrees * math.pi / 180.0

            c = math.cos(angle)
            s = math.sin(angle)
            one_minus_c = 1.0 - c

            """
            Fill the upper-left 3x3 block using Rodrigues' formula.
            """
            cur[0][0] = c + ux * ux * one_minus_c
            cur[0][1] = ux * uy * one_minus_c - uz * s
            cur[0][2] = ux * uz * one_minus_c + uy * s

            cur[1][0] = uy * ux * one_minus_c + uz * s
            cur[1][1] = c + uy * uy * one_minus_c
            cur[1][2] = uy * uz * one_minus_c - ux * s

            cur[2][0] = uz * ux * one_minus_c - uy * s
            cur[2][1] = uz * uy * one_minus_c + ux * s
            cur[2][2] = c + uz * uz * one_minus_c

        """
        Compose transformations.

        If previous transformations are represented by F
        and current transformation is C, then:

            new_F = C * F

        because the current transformation is applied after all previous ones.
        """
        final_matrix = multiply(cur, final_matrix)

    m = int(data[idx])
    idx += 1

    output = []

    for _ in range(m):
        x = float(data[idx])
        y = float(data[idx + 1])
        z = float(data[idx + 2])
        idx += 3

        """
        Apply final_matrix to the homogeneous vector [x, y, z, 1]^T.
        """
        nx = (
            final_matrix[0][0] * x +
            final_matrix[0][1] * y +
            final_matrix[0][2] * z +
            final_matrix[0][3]
        )

        ny = (
            final_matrix[1][0] * x +
            final_matrix[1][1] * y +
            final_matrix[1][2] * z +
            final_matrix[1][3]
        )

        nz = (
            final_matrix[2][0] * x +
            final_matrix[2][1] * y +
            final_matrix[2][2] * z +
            final_matrix[2][3]
        )

        nx = fix_zero(nx)
        ny = fix_zero(ny)
        nz = fix_zero(nz)

        output.append(f"{nx:.2f} {ny:.2f} {nz:.2f}")

    sys.stdout.write("\n".join(output))


if __name__ == "__main__":
    main()
```
