今天这场的出题人也比较仁慈，塞了几道 Div2A 难度的题（

A - An Experiment in Optics Lab

【题意】
从原点开始，xx 正半轴上依次摆放 n≤105n\le 10^5 面棱镜，棱镜 ii 的宽度为 li∈[1,1000]l_i\in [1,1000]、折射系数为 ri∈[104,1.6×104]r_i\in [10^4, 1.6\times 10^4]. 考察分界处，光的折射与左右两面棱镜 i−1,ii-1, i 满足：

ri−1sin⁡θi−1=risin⁡θir_{i-1}\sin \theta_{i-1}=r_i\sin\theta_i

其中 θi−1\theta_{i-1} 是入射角，θi\theta_i 是折射角。

然后，你需要维护 q≤105q\le 10^5 次操作，每次操作会是两种中的一个：

  1. 1 id l r: 用一块新棱镜替换原来的第 idid 块棱镜，其宽度为 ll，折射系数为 rr
  2. 2 p angle: 查询操作。假设有一束光从 (p,0)(p,0) 的位置、与 xx 正半轴成 angle3600\frac{angle}{3600} 度夹角（在 xx 轴上方）发射出去，当这束光从最后一面棱镜出来的时候，折射点的 yy 坐标是多少。

查询的精度在 10−610^{-6} 以内，时间 4s4s.

我们先梳理一下我们可以推出来的信息。考虑在一块棱镜 ii 内部，一束光从其左侧以 θi\theta_i 的角度到达右侧，则这束光的折射点升高了 Δy=litan⁡θi\Delta y=l_i\tan\theta_i，也就是说，我们最后的答案是

∑ilitan⁡θi \sum_i l_i\tan\theta_i

注意到第 ii 次折射的折射角等于第 i+1i+1 次折射的入射角，所以，我们可以在 sin⁡θi\sin\theta_i 之间转移。现在的问题就是如何从 sin⁡θi\sin\theta_i 的关系推理出 tan⁡θi\tan\theta_i 之间的关系

sin⁡θi\sin\theta_i 之间可以以 O(1)O(1) 的方式转移，因为根据折射公式，sin⁡θi+1=riri+1sin⁡θi\sin\theta_{i+1}=\frac{r_i}{r_{i+1}}\sin\theta_i，但是我们从中却推不出 tan⁡θi+1=f(tan⁡θi)\tan\theta_{i+1}=f(\tan\theta_i)，这个 f()f() 怎么样都得带个开根号，我们就很难用数据结构快速维护了，线段树只能快速维护一堆东西求和。

所以，一个解决办法是，把 tan⁡θi\tan\theta_i 对 sin⁡θi\sin\theta_i 进行泰勒展开！这样 tan⁡θi=∑jajsin⁡jθi\tan\theta_i=\sum_j a_j\sin^j\theta_i，而 aja_j 是固定的，这意味着我们就可以用多棵线段树分别维护 sin⁡jθi\sin^j\theta_i

这就是算法的 backbone: 我们考虑这样一个向量乘法，我们假设先用泰勒展开的前三项对 tan⁡(x)=f(sin⁡(x))\tan(x)=f(\sin(x)) 近似

[yiaisin⁡θia2sin⁡2θia3sin⁡3θi]T[1000liriri+100li0(riri+1)20li00(riri+1)3]=[yi+li(aisin⁡θi+a2sin⁡2θi+a3sin⁡3θi)a1riri+1sin⁡θia2(riri+1)2sin⁡2θia3(riri+1)3sin⁡3θi]T=[yi+1aisin⁡θi+1a2sin⁡2θi+1a3sin⁡3θi+1]T \begin{bmatrix} y_i\\a_i\sin\theta_i\\a_2\sin^2\theta_i\\a_3\sin^3\theta_i \end{bmatrix}^T \begin{bmatrix} 1&0&0&0\\ l_i&\frac{r_i}{r_{i+1}}&0&0\\ l_i&0&(\frac{r_i}{r_{i+1}})^2&0\\ l_i&0&0&(\frac{r_i}{r_{i+1}})^3\\ \end{bmatrix}\\ =\begin{bmatrix} y_i+l_i(a_i\sin\theta_i+a_2\sin^2\theta_i+a_3\sin^3\theta_i)\\a_1\frac{r_i}{r_{i+1}}\sin\theta_i\\a_2(\frac{r_i}{r_{i+1}})^2\sin^2\theta_i\\a_3(\frac{r_i}{r_{i+1}})^3\sin^3\theta_i \end{bmatrix}^T\\ =\begin{bmatrix} y_{i+1}\\a_i\sin\theta_{i+1}\\a_2\sin^2\theta_{i+1}\\a_3\sin^3\theta_{i+1} \end{bmatrix}^T

观察结果向量，我们从中发现了可以递推的关系，并且乘的矩阵 MiM_i 本身只和第 i,i+1i,i+1 个棱镜有关，而与入射角无关，经过多面棱镜，就相当于接着乘 Mi+1,Mi+2,…M_{i+1},M_{i+2},\dots，所以我们可以线段树维护 MiM_i，修改棱镜就相当于在线段树上单点修改两次。

时间复杂度 O(kn+kqlog⁡n)O(kn+kq\log n). 其中 kk 是泰勒展开的项数，越高越精确，也越慢。

AC Code

需要注意几点：

不要真的写矩阵乘法，这样的话矩阵乘就是 O(K3)O(K^3) 的。实际上除了第一列和主对角线，MiM_i 剩下的位置都是 00，而且这样两个矩阵乘起来 MiMi+1M_iM_{i+1} 也依然满足除了第一列和主对角线剩下元素都是零. 所以我们可以 O(K)O(K) 记录、计算矩阵乘法。详见 combine(Matrix, Matrix) 函数。

处理查询的时候，首先要判断出初始光线在那个棱镜内部，然后先让这条光线在棱镜内部运动，直到碰到右边界。后面的交给矩阵乘法处理。这里我用二分+树状数组找光源在哪个棱镜内部，所以我的实现其实是 O(kn+qlog⁡2n+kqlog⁡n)O(kn+q\log^2 n+kq\log n) 的，O(kn)O(kn) 是线段树的初始化，O(qlog⁡2n)O(q\log^2 n) 是二分+树状数组找棱镜，O(kqlog⁡n)O(kq\log n) 则是从线段树里提取出矩阵并 combine()。

其次注意精度问题，一是用 long double，二是我这里用了泰勒展开前 K=40K=40 项近似 f(x)=x1−x2f(x)=\frac{x}{\sqrt{1-x^2}}. 这个式子的泰勒展开比较有规律，只有 sin⁡jx\sin^j x 指数为奇数的项，而且通项为

f(x)=x+∑n=1∞1⋅3⋅5⋅⋯⋅(2n−1)2nn!x2n−1 f(x)=x+\sum_{n=1}^{\infin}\frac{1\cdot 3\cdot 5\cdot \dots\cdot (2n-1)}{2^n n!}x^{2n-1}

最后，为了提速，Row, Matrix 最好都用 std::array<ld, K> 而不是 std::vector<ld>.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast,unroll-loops")
#include <bits/stdc++.h>
using ld        = long double;
constexpr int N = 1e5 + 10, K = 40;
constexpr ld pi = 3.14159265358979323846;
std::vector<ld> C;

int n, q;
struct Row {
    ld y;
    std::array<ld, K> sines;
};
struct Mirror {
    ld l, r{1};
} m[N];
struct Matrix {
    std::array<ld, K> col, diag;

    void set(int i) {
        auto &x = m[i];
        for (int i = 0; i < K; i++) col[i] = x.l;

        ld v = m[i].r / m[i + 1].r, t = v;
        for (int i = 0; i < K; i++, t = t * v * v) diag[i] = t;
    }

    friend Matrix combine(const Matrix &lhs, const Matrix &rhs) {
        Matrix res;
        int k = lhs.col.size();

        for (int i = 0; i < lhs.col.size(); i++) {
            res.col[i]  = lhs.col[i] + lhs.diag[i] * rhs.col[i];
            res.diag[i] = lhs.diag[i] * rhs.diag[i];
        }
        return res;
    }
    friend Row combine(const Row &lhs, const Matrix &rhs) {
        Row res;
        res.y = lhs.y;
        for (int i = 0; i < K; i++) res.y += lhs.sines[i] * rhs.col[i];
        for (int i = 0; i < K; i++) res.sines[i] = lhs.sines[i] * rhs.diag[i];
        return res;
    }
};
Matrix T[N * 4];

void build(int p = 1, int l = 1, int r = n) {
    if (l == r) {
        T[p].set(l);
        return;
    }
    int mid = l + ((r - l) >> 1);
    build(p * 2, l, mid), build(p * 2 + 1, mid + 1, r);
    T[p] = combine(T[p * 2], T[p * 2 + 1]);
}
void modify(int q, int p = 1, int l = 1, int r = n) {
    if (l == r) {
        T[p].set(l);
        return;
    }
    int mid = l + ((r - l) >> 1);
    if (q <= mid) modify(q, p * 2, l, mid);
    else modify(q, p * 2 + 1, mid + 1, r);
    T[p] = combine(T[p * 2], T[p * 2 + 1]);
}
Matrix range(int ql, int qr, int p = 1, int l = 1, int r = n) {
    if (ql <= l && r <= qr) return T[p];
    int mid = l + ((r - l) >> 1);
    if (qr <= mid) return range(ql, qr, p * 2, l, mid);
    else if (ql > mid) return range(ql, qr, p * 2 + 1, mid + 1, r);
    else {
        auto L = range(ql, qr, p * 2, l, mid);
        auto R = range(ql, qr, p * 2 + 1, mid + 1, r);
        return combine(L, R);
    }
}

std::array<int, N> b;

int sum(int p) {
    int res = 0;
    for (; p > 0; p -= p & -p) res += b[p];
    return res;
}
int ask(int l, int r) { return l > r ? 0 : sum(r) - sum(l - 1); }
int at(int p) { return ask(p, p); }
void add(int p, int v) {
    for (; p <= n; p += p & -p) b[p] += v;
}
void set(int p, int v) { add(p, -at(p)), add(p, v); }

int start_from(int p) {
    int l = 0, r = n + 1;
    while (l + 1 < r) {
        int mid = l + ((r - l) >> 1);
        if (sum(mid) > p) r = mid;
        else l = mid;
    }
    return r;
}

int main() {
#ifdef ONLINE_JUDGE
    std::cin.tie(0)->sync_with_stdio(0);
    std::cout.tie(0);
#endif
    //! initialize C
    C.push_back(1), C.push_back(0.5);
    for (int i = 2; i < K; i++) C.push_back(C.back() * (2 * i - 1) / 2 / i);

    std::cin >> n;

    for (int i = 1, L; i <= n; i++) {
        std::cin >> L >> m[i].r;
        add(i, L);
        m[i].l = L;
    }
    build(1, 1, n);

    std::cin >> q;
    while (q--) {
        int op, id, p, deg;
        std::cin >> op;
        if (op == 2) {
            std::cin >> p >> deg;

            int r    = start_from(p);
            ld theta = pi * deg / 3600 / 180;

            int dist = sum(r) - p;

            ld ans = (ld)dist * std::tan(theta);
            if (r + 1 <= n) {
                ld ntheta = m[r].r * std::sin(theta) / m[r + 1].r;

                Row initial;
                initial.y = ans;
                ld v = ntheta, t = v;
                for (int i = 0; i < K; i++, t = t * v * v) initial.sines[i] = t * C[i];

                auto mat = range(r + 1, n);
                initial  = combine(initial, mat);
                ans      = initial.y;
            }
            std::cout << std::setprecision(20) << std::fixed << ans << "\n";
        } else {
            std::cin >> id;
            std::cin >> m[id].l >> m[id].r;
            modify(id);
            if (id > 1) modify(id - 1);
            set(id, int(m[id].l));
        }
    }
}