ICPC Regional Hanoi 2024 Editorial - A

Cho một đa giác đơn không tự cắt $P$ gồm $n$ đỉnh trên mặt phẳng 2D. Có $q$ truy vấn độc lập. Mỗi truy vấn mô tả một cơn bão là một hình tròn di chuyển từ tâm $S(x_s, y_s)$ với bán kính $r_s$ đến tâm $T(x_t, y_t)$ với bán kính $r_t$ theo đường thẳng:

1. Bán kính biến thiên tuyến tính theo khoảng cách di chuyển giữa $S$ và $T$.

2. Miền bị ảnh hưởng $H$ bởi cơn bão là tổng thể diện tích được bao phủ bởi hình tròn tại mọi thời điểm trong suốt hành trình di chuyển từ $S$ đến $T$ (về mặt hình học, đây là hợp của hai hình tròn tại $S, T$ và hình thang cân nối các tiếp tuyến chung ngoài của hai đường tròn này).

Với mỗi truy vấn, tính diện tích phần giao giữa đa giác đảo $P$ và miền ảnh hưởng $H$ của cơn bão.

Giới hạn: số test $t \le 10^4$, $3 \le n \le 3 \cdot 10^5$, tọa độ có trị tuyệt đối không quá $5000$, $1 \le r_s, r_t \le 100$, tổng $n \cdot q$ trên mọi test không vượt quá $10^6$. Đáp án được chấp nhận nếu sai số tuyệt đối hoặc tương đối không vượt quá $10^{-6}$.

Trong kì thi chính thức không có đội nào AC bài toán này.

Ý tưởng chính

Do tổng $n \cdot q$ trên mọi test không vượt quá $10^6$, với mỗi truy vấn ta có thể duyệt qua toàn bộ $n$ cạnh của đa giác. Bài toán có hai khó khăn chính: đa giác $P$ có thể không lồi, và biên của miền $H$ có những đoạn cong. Lời giải lần lượt xử lý hai khó khăn này qua bốn bước:

1. Ở Bước 1, ta biểu diễn đa giác $P$ thành tổng có dấu của $n$ tam giác, mỗi tam giác ứng với một cạnh của đa giác. Nhờ đó, bài toán quy về việc tính diện tích phần giao giữa một tam giác và miền $H$.

2. Ở Bước 2, ta chứng minh $H$ là bao lồi của hai hình tròn lúc đầu và lúc cuối, sau đó chia $H$ thành ba mảnh đơn giản: hai phần hình tròn và một hình thang.

3. Ở Bước 3, ta cắt mỗi tam giác theo ba mảnh nói trên. Phần nằm trong hình thang là một đa giác lồi, nên diện tích của phần này được tính bằng công thức shoelace.

4. Ở Bước 4, ta tính diện tích phần giao giữa một đa giác lồi và một hình tròn, áp dụng cho hai phần còn lại.

Trong lời giải, ta dùng các ký hiệu sau: $\mathcal{A}(X)$ là diện tích có hướng của hình $X$; $p \times q = p_x q_y - p_y q_x$ là tích có hướng của hai vector $p$ và $q$; $D_S$ là hình tròn tâm $S$ bán kính $r_s$ và $D_T$ là hình tròn tâm $T$ bán kính $r_t$, tức là hình tròn lúc cơn bão bắt đầu và lúc cơn bão kết thúc.

Bước 1: Biểu diễn đa giác thành tổng có dấu của các tam giác

Ta chọn một điểm $O$ bất kỳ làm điểm gốc (chẳng hạn trung điểm của $ST$). Với mỗi cạnh $v_i v_{i+1}$ của đa giác (quy ước $v_n = v_0$), ta xét tam giác $O v_i v_{i+1}$ cùng với dấu của tam giác đó: tam giác mang dấu dương nếu các đỉnh $O$, $v_i$, $v_{i+1}$ được sắp theo chiều ngược chiều kim đồng hồ, và mang dấu âm trong trường hợp ngược lại. Công thức shoelace quen thuộc chính là tổng diện tích có hướng của $n$ tam giác này. Tuy nhiên, để xử lý được phần giao với $H$, ta cần một tính chất mạnh hơn, phát biểu cho từng điểm của mặt phẳng.

Nhận xét. Với mỗi điểm $X$ không nằm trên cạnh của tam giác nào, gọi $w(X)$ là tổng dấu của các tam giác chứa $X$: mỗi tam giác dương chứa $X$ được tính $+1$, mỗi tam giác âm chứa $X$ được tính $-1$. Nếu đa giác được cho theo chiều ngược chiều kim đồng hồ thì $w(X) = 1$ với mọi điểm $X$ nằm trong đa giác, và $w(X) = 0$ với mọi điểm $X$ nằm ngoài đa giác. Nếu đa giác được cho theo chiều kim đồng hồ thì giá trị ứng với các điểm nằm trong đa giác là $-1$.

Chứng minh. Gọi $\ell$ là tia xuất phát từ $X$ và có cùng hướng với vector $\overrightarrow{OX}$, tức là tia đi ra xa điểm $O$. Điểm $X$ nằm trong tam giác $O v_i v_{i+1}$ khi và chỉ khi đường thẳng $OX$ cắt cạnh $v_i v_{i+1}$ tại một điểm nằm xa $O$ hơn $X$, tức là khi và chỉ khi tia $\ell$ cắt cạnh $v_i v_{i+1}$. Ngoài ra, dấu của tam giác $O v_i v_{i+1}$ cho biết cạnh $v_i v_{i+1}$, đi theo chiều từ $v_i$ đến $v_{i+1}$, cắt tia $\ell$ từ bên phải sang bên trái hay từ bên trái sang bên phải. Do đó $w(X)$ bằng số lần biên đa giác cắt tia $\ell$ theo chiều này trừ đi số lần cắt theo chiều kia.

Bây giờ ta đi dọc tia $\ell$ từ $X$ ra xa vô cùng. Mỗi lần tia cắt biên, ta đi từ trong đa giác ra ngoài hoặc từ ngoài vào trong. Vì phần trong của đa giác luôn nằm về cùng một phía của các cạnh (bên trái nếu đa giác được cho ngược chiều kim đồng hồ), hai loại lần cắt này ứng với hai chiều cắt ngược nhau. Đầu xa của tia nằm ngoài đa giác. Nếu $X$ nằm ngoài đa giác, số lần đi vào bằng số lần đi ra, nên $w(X) = 0$. Nếu $X$ nằm trong đa giác, số lần đi ra nhiều hơn số lần đi vào đúng một lần, nên $w(X) = \pm 1$. Đây cũng chính là cơ sở của thuật toán kiểm tra một điểm có nằm trong đa giác hay không bằng cách bắn tia. $\blacksquare$

Theo nhận xét trên, khi cộng diện tích có hướng của các phần $O v_i v_{i+1} \cap H$, mỗi điểm $X$ thuộc $H$ được tính đúng $w(X)$ lần: các điểm nằm ngoài đa giác bị triệt tiêu, còn mỗi điểm nằm trong đa giác được tính đúng một lần, với dấu phụ thuộc vào chiều của đa giác. Do đó:

Vì đa giác có thể được cho theo chiều kim đồng hồ, đáp án của truy vấn là trị tuyệt đối của tổng trên. Như vậy, bài toán quy về việc tính diện tích phần giao giữa một tam giác và miền $H$.

Bước 2: Hình dạng của miền ảnh hưởng

Đặt $d = |ST|$ và $\Delta = r_t - r_s$. Khi tâm cơn bão đã đi được một tỉ lệ $\lambda \in [0, 1]$ của quãng đường $ST$, khoảng cách từ tâm đến $S$ và đến $T$ lần lượt là $d_s = \lambda d$ và $d_t = (1 - \lambda) d$. Theo công thức trong đề bài, tâm cơn bão lúc đó là điểm $(1-\lambda)S + \lambda T$, và bán kính là

Như vậy, cả tâm và bán kính đều thay đổi tuyến tính theo $\lambda$. Gọi $C_\lambda$ là hình tròn tại thời điểm đó. Sau khi tâm tới $T$, cơn bão đứng yên nên không phủ thêm điểm nào, do đó $H$ là hợp của các hình tròn $C_\lambda$ với $\lambda \in [0, 1]$.

Nhận xét. Miền $H$ là bao lồi của $D_S \cup D_T$.

Chứng minh. Trước hết, ta chứng minh mỗi hình tròn $C_\lambda$ gồm đúng các điểm có dạng $(1-\lambda)a + \lambda b$ với $a \in D_S$ và $b \in D_T$.

• Nếu $a \in D_S$ và $b \in D_T$, khoảng cách từ điểm $(1-\lambda)a + \lambda b$ đến tâm $(1-\lambda)S + \lambda T$ bằng $|(1-\lambda)(a - S) + \lambda(b - T)|$. Theo bất đẳng thức tam giác, đại lượng này không vượt quá $(1-\lambda)r_s + \lambda r_t$, nên điểm đang xét thuộc $C_\lambda$.

• Ngược lại, mỗi điểm của $C_\lambda$ có dạng $(1-\lambda)S + \lambda T + w$ với $|w| \le r$, trong đó $r = (1-\lambda) r_s + \lambda r_t$ là bán kính của $C_\lambda$. Đặt $a = S + \frac{r_s}{r} w$ và $b = T + \frac{r_t}{r} w$. Khi đó $a \in D_S$, $b \in D_T$, và $(1-\lambda)a + \lambda b = (1-\lambda)S + \lambda T + w$.

Lấy hợp theo mọi $\lambda \in [0, 1]$, ta được $H$ là tập các điểm $(1-\lambda)a + \lambda b$ với $a \in D_S$, $b \in D_T$ và $\lambda \in [0, 1]$. Đây chính là bao lồi của hợp hai tập lồi $D_S$ và $D_T$. $\blacksquare$

Nếu $d \le |\Delta|$, một hình tròn nằm trọn trong hình tròn kia, nên $H$ chính là hình tròn lớn hơn. Trường hợp này được xử lý riêng ở phần Cài đặt. Trong phần còn lại của Bước 2 và Bước 3, ta giả sử $d > |\Delta|$. Khi đó, bao lồi của hai hình tròn được giới hạn bởi hai tiếp tuyến chung ngoài và hai cung tròn ở hai đầu (Hình 2). Có thể hình dung trực quan như sau: khoảng cách có hướng từ một điểm đến một đường thẳng là hàm bậc nhất theo tọa độ của điểm, nên khoảng cách từ tâm của $C_\lambda$ đến mỗi tiếp tuyến chung ngoài thay đổi tuyến tính theo $\lambda$, bằng $r_s$ khi $\lambda = 0$ và bằng $r_t$ khi $\lambda = 1$. Vì vậy khoảng cách này luôn bằng bán kính của $C_\lambda$, tức là mọi hình tròn trung gian đều tiếp xúc với cả hai tiếp tuyến chung ngoài.

Để thuận tiện cho việc tính toán, ta đổi hệ tọa độ: gốc tọa độ là trung điểm $O$ của đoạn $ST$, trục $Ox$ cùng hướng với vector $\overrightarrow{ST}$. Trong hệ tọa độ mới, $S = \left(-\frac{d}{2}, 0\right)$ và $T = \left(\frac{d}{2}, 0\right)$.

Ta tìm tiếp tuyến chung ngoài phía trên dưới dạng $\mathbf{n} \cdot p = c$, trong đó $\mathbf{n}$ là vector pháp tuyến đơn vị hướng ra phía ngoài $H$. Với một điểm $p$ nằm cùng phía với $H$ so với đường thẳng này, khoảng cách từ $p$ đến đường thẳng bằng $c - \mathbf{n} \cdot p$. Đường thẳng tiếp xúc với $D_S$ và $D_T$ khi và chỉ khi khoảng cách từ $S$ và $T$ đến đường thẳng lần lượt bằng $r_s$ và $r_t$:

Trừ vế theo vế hai đẳng thức, ta được $\mathbf{n} \cdot (T - S) = r_s - r_t$. Vì $T - S = (d, 0)$, đẳng thức này cho $n_x = -\frac{\Delta}{d}$. Do $\mathbf{n}$ là vector đơn vị hướng lên trên, $n_y = \sqrt{1 - n_x^2} = \frac{h}{d}$ với $h = \sqrt{d^2 - \Delta^2}$. Thay vào đẳng thức thứ nhất, ta tìm được hằng số $c$:

Do tính đối xứng qua trục $Ox$, tiếp tuyến chung ngoài phía dưới có phương trình $\mathbf{n}' \cdot p = c$, với $\mathbf{n}' = \left(-\frac{\Delta}{d}, -\frac{h}{d}\right)$ và cùng hằng số $c$.

Tiếp điểm của một đường tròn với một tiếp tuyến là chân đường vuông góc hạ từ tâm xuống tiếp tuyến đó. Vì vậy, hai tiếp điểm trên $D_S$ là $S + r_s \mathbf{n}$ và $S + r_s \mathbf{n}'$, còn hai tiếp điểm trên $D_T$ là $T + r_t \mathbf{n}$ và $T + r_t \mathbf{n}'$. Hai vector $\mathbf{n}$ và $\mathbf{n}'$ có cùng hoành độ, nên hai tiếp điểm trên cùng một đường tròn nằm trên cùng một đường thẳng đứng. Hoành độ của hai đường thẳng đứng này lần lượt là

và ta có $x_R - x_L = \frac{h^2}{d} > 0$. Hai đường thẳng đứng $x = x_L$ và $x = x_R$ chia $H$ thành ba mảnh (Hình 2):

1. Trong nửa mặt phẳng $x \le x_L$, biên của $H$ là cung tròn của $D_S$ nối hai tiếp điểm. Do đó, phần của $H$ nằm trong nửa mặt phẳng này trùng với phần của $D_S$ nằm trong nửa mặt phẳng đó.

2. Trong dải $x_L \le x \le x_R$, biên của $H$ gồm hai đoạn tiếp tuyến. Do đó, phần của $H$ nằm trong dải này là một hình thang cân, chính là giao của bốn nửa mặt phẳng $x \ge x_L$, $x \le x_R$, $\mathbf{n} \cdot p \le c$ và $\mathbf{n}' \cdot p \le c$.

3. Tương tự, trong nửa mặt phẳng $x \ge x_R$, phần của $H$ trùng với phần của $D_T$.

Bước 3: Cắt một tam giác theo ba mảnh của miền ảnh hưởng

Xét một tam giác $\triangle = O v_i v_{i+1}$ ở Bước 1. Ta ký hiệu:

• $\triangle_L$ là phần của tam giác nằm trong nửa mặt phẳng $x \le x_L$;

• $\triangle_M$ là phần của tam giác nằm trong hình thang ở Bước 2;

• $\triangle_R$ là phần của tam giác nằm trong nửa mặt phẳng $x \ge x_R$.

Theo cách chia ở Bước 2, phần chung của $\triangle_L$ và $H$ chính là phần chung của $\triangle_L$ và $D_S$; tương tự, phần chung của $\triangle_R$ và $H$ là phần chung của $\triangle_R$ và $D_T$; còn $\triangle_M$ nằm trọn trong $H$. Do đó:

Số hạng ở giữa là diện tích của một đa giác lồi, được tính bằng công thức shoelace. Hai số hạng còn lại là diện tích phần giao giữa một đa giác lồi và một hình tròn, sẽ được tính ở Bước 4.

Cả ba đa giác $\triangle_L$, $\triangle_M$ và $\triangle_R$ đều thu được bằng cách cắt một đa giác lồi bởi một nửa mặt phẳng, lặp lại một số lần: $\triangle_L$ và $\triangle_R$ cần một lần cắt, còn $\triangle_M$ cần bốn lần cắt, ứng với bốn nửa mặt phẳng tạo nên hình thang. Thao tác cắt này được thực hiện bằng thuật toán Sutherland–Hodgman (hàm clip_halfplane trong code), được trình bày chi tiết trong khung bên dưới. Mỗi lần cắt cho ra một đa giác lồi có cùng chiều với đa giác ban đầu, nên dấu của diện tích có hướng được giữ nguyên; đồng thời, số đỉnh tăng thêm không quá một. Vì vậy, $\triangle_M$ có không quá $3 + 4 = 7$ đỉnh.

Thuật toán Sutherland–Hodgman

Bài toán. Cho đa giác lồi $Q = q_0 q_1 \dots q_{m-1}$, với các đỉnh được liệt kê theo thứ tự trên biên, và nửa mặt phẳng $\{p : \mathbf{n} \cdot p \le c\}$. Ta cần tìm đa giác là phần chung của $Q$ và nửa mặt phẳng này.

Phân loại điểm. Đặt $f(p) = \mathbf{n} \cdot p - c$. Điểm $p$ nằm trong nửa mặt phẳng nếu $f(p) \le 0$, và nằm ngoài nếu $f(p) > 0$. Vì $f$ là hàm bậc nhất, trên một cạnh $xy$ giá trị của $f$ thay đổi tuyến tính từ $f(x)$ đến $f(y)$. Do đó, cạnh $xy$ cắt đường thẳng biên $f(p) = 0$ khi và chỉ khi hai đầu mút $x$ và $y$ nằm ở hai phía khác nhau.

Thuật toán. Ta đi một vòng quanh biên của $Q$, lần lượt xét các cạnh $q_0 \to q_1$, $q_1 \to q_2$, …, $q_{m-1} \to q_0$. Với mỗi cạnh $x \to y$, ta ghi các điểm sau vào đa giác kết quả, theo đúng thứ tự:

1. Nếu $x$ nằm trong, ta ghi điểm $x$.

2. Nếu $x$ và $y$ nằm ở hai phía khác nhau, ta ghi giao điểm của cạnh $xy$ với đường thẳng biên.

Giao điểm có dạng $x + t(y - x)$ với $t \in [0, 1]$. Vì $f(x + t(y - x)) = f(x) + t\big(f(y) - f(x)\big)$, giao điểm ứng với giá trị

$$t = \frac{f(x)}{f(x) - f(y)}$$

Ví dụ. Hình 4 minh họa việc cắt tứ giác $ABCD$, trong đó $A$ và $B$ nằm trong nửa mặt phẳng, còn $C$ và $D$ nằm ngoài. Thuật toán xét bốn cạnh như sau:

Cạnh	Đầu cạnh	Cuối cạnh	Các điểm được ghi
$A \to B$	trong	trong	$A$
$B \to C$	trong	ngoài	$B$, sau đó giao điểm $I_1$
$C \to D$	ngoài	ngoài	không có
$D \to A$	ngoài	trong	giao điểm $I_2$

Kết quả là tứ giác $A B I_1 I_2$. Các đỉnh nằm trong ($A$ và $B$) được giữ nguyên theo thứ tự cũ, còn phần biên nằm ngoài ($I_1 \to C \to D \to I_2$) được thay bằng đoạn thẳng $I_1 I_2$ nằm trên đường thẳng biên.

Tính đúng đắn. Phần chung của $Q$ và nửa mặt phẳng là giao của hai tập lồi, nên cũng là một đa giác lồi. Nếu mọi đỉnh của $Q$ nằm trong nửa mặt phẳng, kết quả là toàn bộ $Q$; nếu mọi đỉnh nằm ngoài, kết quả là tập rỗng; thuật toán cho đúng kết quả trong cả hai trường hợp này. Trong trường hợp còn lại, đường thẳng biên cắt biên của $Q$ tại đúng hai điểm và chia biên của $Q$ thành hai đường gấp khúc: một đường nằm trong nửa mặt phẳng và một đường nằm ngoài. Biên của đa giác kết quả gồm đường gấp khúc nằm trong, được khép lại bởi đoạn thẳng nối hai giao điểm. Khi đi một vòng quanh $Q$, thuật toán ghi lại đúng các đỉnh của đường gấp khúc nằm trong, cùng với hai giao điểm ở hai đầu, theo thứ tự xuất hiện trên biên. Vì vậy, đa giác kết quả có cùng chiều với $Q$. Hơn nữa, có ít nhất một đỉnh của $Q$ bị loại và chỉ có hai giao điểm được thêm vào, nên đa giác kết quả có không quá $(m - 1) + 2 = m + 1$ đỉnh.

Để cắt một đa giác bởi nhiều nửa mặt phẳng, ta áp dụng thuật toán lần lượt với từng nửa mặt phẳng. Hàm clip_halfplane trong code cài đặt đúng thuật toán trên, kèm theo hai chi tiết để tránh sai số: một điểm nằm ngoài nhưng rất gần đường thẳng biên (sai lệch cỡ $10^{-16}$, tính tương đối theo độ lớn của các giá trị) vẫn được coi là nằm trong, và tham số $t$ được giới hạn trong đoạn $[0, 1]$.

Bước 4: Diện tích giao giữa đa giác lồi và hình tròn

Cho đa giác lồi $Q = q_0 q_1 \dots q_{m-1}$ và hình tròn $D$ có tâm $C$, bán kính $r$. Nhận xét ở Bước 1 đúng với mọi điểm gốc và mọi hình được đem giao, nên ta có thể áp dụng nhận xét đó với điểm gốc là tâm $C$ và hình $D$ thay cho $H$:

Như vậy, ta chỉ cần tính diện tích có hướng của phần giao giữa hình tròn $D$ và một tam giác $CAB$ có một đỉnh là tâm $C$. Ta đặt gốc tọa độ tại $C$. Đoạn $AB$ gồm các điểm $A + t(B - A)$ với $t \in [0, 1]$, và các giao điểm của đoạn $AB$ với đường tròn ứng với các nghiệm $t \in (0, 1)$ của phương trình bậc hai

Phương trình này có nhiều nhất hai nghiệm, nên các giao điểm chia đoạn $AB$ thành nhiều nhất ba đoạn con. Mỗi đoạn con $pq$ không chứa giao điểm nào ở bên trong, nên nằm hoàn toàn bên trong hoặc hoàn toàn bên ngoài hình tròn. Tương ứng, tam giác $CAB$ được chia thành các tam giác nhỏ $Cpq$ (Hình 5), và ta tính phần giao của từng tam giác nhỏ với hình tròn:

• Nếu đoạn $pq$ nằm trong hình tròn, cả ba điểm $C$, $p$, $q$ đều thuộc hình tròn. Vì hình tròn là một tập lồi, toàn bộ tam giác $Cpq$ nằm trong hình tròn, nên phần giao có diện tích có hướng bằng $\frac{1}{2}(p \times q)$.

• Nếu đoạn $pq$ nằm ngoài hình tròn, mỗi tia xuất phát từ $C$ và đi qua một điểm của đoạn $pq$ sẽ ra khỏi hình tròn trước khi chạm tới đoạn $pq$. Do đó, phần giao chính là hình quạt giới hạn bởi hai tia $Cp$ và $Cq$, có diện tích có hướng bằng $\frac{1}{2} r^2 \theta$, trong đó $\theta = \operatorname{atan2}(p \times q,\ p \cdot q)$ là góc có hướng từ vector $p$ đến vector $q$.

Để biết một đoạn con nằm trong hay nằm ngoài hình tròn, ta chỉ cần kiểm tra một điểm bất kỳ ở bên trong đoạn đó. Code kiểm tra hai điểm, tại vị trí $\frac{1}{3}$ và $\frac{2}{3}$ của đoạn, để hạn chế ảnh hưởng của sai số.

Cài đặt

Với mỗi truy vấn, thuật toán thực hiện như sau:

1. Nếu $d \le |\Delta|$, miền $H$ là hình tròn lớn hơn trong hai hình tròn. Ta chọn $O$ là tâm của hình tròn này và cộng các đại lượng $\mathcal{A}(O v_i v_{i+1} \cap H)$, trong đó mỗi đại lượng được tính bằng Bước 4.

2. Trong trường hợp còn lại, ta đổi sang hệ tọa độ ở Bước 2 và tính $x_L$, $x_R$, $\mathbf{n}$, $\mathbf{n}'$ và $c$. Với mỗi $i$, ta cắt tam giác $O v_i v_{i+1}$ thành ba phần như ở Bước 3, rồi cộng ba đại lượng tương ứng vào tổng.

3. Đáp án là trị tuyệt đối của tổng.

Mỗi tam giác chỉ cần một số lượng phép tính không đổi, nên mỗi truy vấn có độ phức tạp $\mathcal{O}(n)$, và độ phức tạp của toàn bộ thuật toán là $\mathcal{O}\left(\sum n \cdot q\right)$.

Độ chính xác: Mỗi số hạng $\mathcal{A}(O v_i v_{i+1} \cap H)$ có thể lớn tới cỡ diện tích của $H$ (khoảng $3 \cdot 10^6$), và các số hạng dương, âm triệt tiêu lẫn nhau, trong khi đáp án có thể rất nhỏ nhưng vẫn cần sai số tuyệt đối không quá $10^{-6}$. Code xử lý vấn đề này bằng ba cách:

1. Dùng kiểu long double thay cho double. Trên máy chấm Linux x86, long double có khoảng 19 chữ số có nghĩa, so với khoảng 16 chữ số của double.

2. Cộng dồn đáp án bằng tổng Kahan (xem khung bên dưới), nhờ đó sai số gần như không tích lũy theo số lượng số hạng.

3. Giải phương trình bậc hai theo cách ổn định số. Xét phương trình $\alpha t^2 + \beta t + \gamma = 0$. Khi $\beta^2$ lớn hơn nhiều so với $|4\alpha\gamma|$, ta có $\sqrt{\beta^2 - 4\alpha\gamma} \approx |\beta|$, nên trong công thức nghiệm $\frac{-\beta \pm \sqrt{\beta^2 - 4\alpha\gamma}}{2\alpha}$, một trong hai nghiệm là hiệu của hai số gần bằng nhau, và phép trừ này làm mất gần hết các chữ số có nghĩa. Để tránh điều đó, code chỉ tính trực tiếp nghiệm ứng với phép cộng hai số cùng dấu, rồi suy ra nghiệm còn lại từ định lý Viète $t_1 t_2 = \frac{\gamma}{\alpha}$:

   $$u = -\frac{1}{2}\left(\beta + \operatorname{sgn}(\beta)\sqrt{\beta^2 - 4\alpha\gamma}\right), \qquad t_1 = \frac{u}{\alpha}, \qquad t_2 = \frac{\gamma}{u}$$

Tổng Kahan

Vấn đề. Máy tính lưu mỗi số thực với một số lượng chữ số có nghĩa cố định (với double là 53 bit, tương đương khoảng 16 chữ số thập phân). Vì vậy, kết quả của mỗi phép cộng đều bị làm tròn. Khi cộng một số nhỏ vào một tổng lớn, các chữ số ở hàng thấp của số nhỏ không được lưu lại. Khi cộng liên tiếp $n$ số, các sai số làm tròn này có thể tích lũy, và trong trường hợp xấu nhất, sai số tăng tỉ lệ với $n$.

Giả sử máy tính chỉ lưu được 4 chữ số có nghĩa (theo hệ thập phân), và ta cần tính $1000 + 1{,}4 + 1{,}4 + 1{,}4$. Kết quả chính xác là $1004{,}2$; khi làm tròn tới 4 chữ số có nghĩa, kết quả tốt nhất có thể đạt được là $1004$. Ở phép cộng đầu tiên, $1000 + 1{,}4 = 1001{,}4$, nhưng máy chỉ lưu được $1001$, nên phần $0{,}4$ bị mất (Hình 6). Nếu cộng theo cách thông thường, cả ba phép cộng đều mất $0{,}4$, và kết quả cuối cùng là $1003$, lệch $1{,}2$ so với kết quả chính xác.

Ý tưởng. Sau mỗi phép cộng, thuật toán Kahan tính lại phần vừa bị làm tròn mất và bù phần đó vào phép cộng tiếp theo. Biến correction lưu phần bị mất của phép cộng gần nhất:

C++

1
void add(ld x) {
2
  ld y = x - correction;
3
  ld next = sum + y;
4
  correction = (next - sum) - y;
5
  sum = next;
6
}

Ý nghĩa của từng dòng như sau:

• y = x - correction là số hạng mới, đã được bù phần bị mất ở phép cộng trước.

• next = sum + y là tổng mới; giá trị này đã bị làm tròn.

• next - sum là phần thực sự đã được cộng vào tổng. Vì next và sum có giá trị gần nhau, phép trừ này thường không bị làm tròn.

• correction = (next - sum) - y là phần thực sự đã được cộng trừ đi phần cần cộng, tức là sai số của phép cộng vừa thực hiện. Đại lượng này âm nếu kết quả bị làm tròn xuống, và dương nếu kết quả bị làm tròn lên. Ở phép cộng tiếp theo, đại lượng này được trừ khỏi số hạng mới, nghĩa là phần bị mất được cộng bù trở lại.

Bảng sau chạy lại ví dụ trên với máy tính chỉ lưu 4 chữ số có nghĩa:

Lần cộng	Cộng thông thường	y	next	correction	Tổng Kahan
1	$1001$	$1{,}4$	$1001$	$-0{,}4$	$1001$
2	$1002$	$1{,}4 + 0{,}4 = 1{,}8$	$1003$	$0{,}2$	$1003$
3	$1003$	$1{,}4 - 0{,}2 = 1{,}2$	$1004$	$-0{,}2$	$1004$

Tổng Kahan cho kết quả $1004$, đúng bằng kết quả chính xác sau khi làm tròn, trong khi cách cộng thông thường cho $1003$. Tổng quát, gọi $\varepsilon$ là sai số tương đối của một phép làm tròn. Khi cộng $n$ số $x_1, x_2, \dots, x_n$, sai số của cách cộng thông thường bị chặn bởi khoảng $n \varepsilon \sum |x_i|$, còn sai số của tổng Kahan bị chặn bởi khoảng $2 \varepsilon \sum |x_i|$, gần như không phụ thuộc vào $n$. Trong bài này, mỗi truy vấn cộng tới $9 \cdot 10^5$ số hạng (ba số hạng cho mỗi tam giác), và các số hạng có giá trị lớn, triệt tiêu lẫn nhau, nên tổng Kahan giúp đáp án chính xác hơn đáng kể.

Lưu ý. Không biên dịch chương trình với cờ -ffast-math, cũng như không dùng #pragma GCC optimize("Ofast"). Với các tùy chọn này, trình biên dịch được phép biến đổi biểu thức theo các quy tắc đại số, chẳng hạn rút gọn (next - sum) - y thành $0$ vì next = sum + y, làm mất hoàn toàn phần bù.