449. 巧克力糖衣(Chocolate Covered Candy)
糖果师傅菲尔正在制作一批裹有巧克力糖衣的新糖果。每块糖的中心是一个旋转椭球体，其方程为 \(b^2 x^2 + b^2 y^2 + a^2 z^2 = a^2 b^2\)。菲尔想知道，给一块这样的糖心裹上厚度均匀为 1 毫米的巧克力糖衣，需要多少巧克力。
若 \(a = 1\) mm 且 \(b = 1\) mm，所需巧克力为 \(\dfrac{28}{3}\pi\) mm\(^3\)。
若 \(a = 2\) mm 且 \(b = 1\) mm，所需巧克力约为 \(60.35475635\) mm\(^3\)。
求当 \(a = 3\) mm 且 \(b = 1\) mm 时所需巧克力的体积（单位 mm\(^3\)），并将结果四舍五入到小数点后 8 位。
分析：均匀厚度的糖衣对应三维凸体的“平行体”（Minkowski 和）。对光滑凸体，Steiner 公式给出外层体积的精确展开，从而把壳层体积化为表面积与平均曲率积分两项，再按扁/长旋转椭球分别代入闭式公式即可。
一、数学背景
1. 椭球参数化
将题面方程两边同除以 \(a^2 b^2\)，得到标准形式
$$
\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1.
$$
因此 \(x,y\) 方向的半轴长度均为 \(a\)（赤道半径），\(z\) 方向的半轴长度为 \(b\)（极半径）。这是一个绕 \(z\) 轴旋转而成的椭球体。
当 \(a=b\) 时为球体；当 \(a>b\) 时为扁椭球（oblate spheroid）；当 \(a<b\) 时为长椭球（prolate spheroid）。
2. 均匀糖衣 = 平行体
厚度为 \(t\) 的均匀涂层，几何上等于原椭球与半径为 \(t\) 的球做 Minkowski 和：
$$
K_t = K \oplus tB.
$$
对三维光滑凸体，Steiner 公式给出
$$
V(K_t)=V_0 + S\,t + M\,t^2 + \frac{4\pi}{3}t^3,
$$
其中 \(V_0\) 为原体积，\(S\) 为表面积，\(M=\int H\,\mathrm{d}A\) 为平均曲率在曲面上的积分。
糖衣体积（壳层）为
$$
V_{\text{choc}} = V(K_t)-V_0 = S\,t + M\,t^2 + \frac{4\pi}{3}t^3.
$$
本题 \(t=1\) mm。
3. 球体校验
当 \(a=b=R\) 时，\(S=4\pi R^2\)，\(M=4\pi R\)，代入得
$$
V_{\text{choc}} = 4\pi R^2 + 4\pi R + \frac{4\pi}{3} = \frac{28}{3}\pi \quad (R=1),
$$
与题面样例一致。
4. 旋转椭球闭式公式
记赤道半径 \(A=a\)，极半径 \(C=b\)。
扁椭球（\(A>C\)），离心率 \(e=\sqrt{1-C^2/A^2}\)：
$$
S = 2\pi A^2\left(1+\frac{C^2}{A^2}\cdot\frac{\operatorname{artanh}(e)}{e}\right),
$$
$$
M = 2\pi C + \frac{2\pi A}{e}\arctan\!\left(\frac{Ae}{C}\right).
$$
长椭球（\(A<C\)），离心率 \(e=\sqrt{1-A^2/C^2}\)：
$$
S = 2\pi A^2\left(1+\frac{C}{A}\cdot\frac{\arcsin(e)}{e}\right),
$$
$$
M = 2\pi C + \frac{2\pi A^2}{Ce}\operatorname{artanh}(e).
$$
目标情形 \(a=3,b=1\) 属于扁椭球分支。
二、算法设计
- 将题面参数 \((a,b)\) 映射为赤道半径 \(A=a\)、极半径 \(C=b\)。
- 判断 \(A\) 与 \(C\) 是否近似相等；若相等，走球体专用公式。
- 否则按 \(A>C\)（扁）或 \(A<C\)（长）选择对应分支，计算 \(S\) 与 \(M\)。
- 用 Steiner 公式计算 \(V_{\text{choc}}=S+M+\dfrac{4\pi}{3}\)（\(t=1\)）。
- 先用题面两组样例自测，再计算 \(a=3,b=1\) 并保留 8 位小数。
数值上，当 \(e\) 很小时，\(\operatorname{artanh}(e)/e\) 与 \(\arcsin(e)/e\) 用泰勒展开避免 \(0/0\) 型不稳定。
三、复杂度分析
- 时间：\(O(1)\)，仅含常数次超越函数调用。
- 空间：\(O(1)\)。
- 实测 3 次运行均在 \(10^{-4}\) 秒量级，远低于 60 秒限制。
四、代码实现与说明
"""Project Euler 449 - Chocolate Covered Candy."""
from __future__ import annotations
import math
from time import perf_counter
def _series_atanh_over_x(x: float) -> float:
"""Return atanh(x)/x with a stable Taylor expansion near zero."""
if abs(x) < 1e-8:
x2 = x * x
return 1.0 + x2 / 3.0 + (x2 * x2) / 5.0
return math.atanh(x) / x
def _series_asin_over_x(x: float) -> float:
"""Return asin(x)/x with a stable Taylor expansion near zero."""
if abs(x) < 1e-8:
x2 = x * x
return 1.0 + x2 / 6.0 + 3.0 * (x2 * x2) / 40.0
return math.asin(x) / x
def chocolate_volume(a: float, b: float, thickness: float = 1.0) -> float:
"""Return chocolate volume for a uniform shell around the candy ellipsoid.
The candy satisfies b^2 x^2 + b^2 y^2 + a^2 z^2 = a^2 b^2, which is
equivalent to x^2/a^2 + y^2/a^2 + z^2/b^2 = 1. Here `a` is the equatorial
semi-axis and `b` is the polar semi-axis. A shell of thickness `thickness`
is the outer parallel body minus the original ellipsoid volume.
"""
if a <= 0.0 or b <= 0.0 or thickness < 0.0:
raise ValueError("a and b must be positive; thickness must be non-negative")
pi = math.pi
t = thickness
if abs(a - b) < 1e-14:
surface_area = 4.0 * pi * a * a
mean_curvature_integral = 4.0 * pi * a
return surface_area * t + mean_curvature_integral * t * t + (4.0 * pi / 3.0) * t**3
equatorial = a
polar = b
if equatorial > polar:
eccentricity = math.sqrt(1.0 - (polar * polar) / (equatorial * equatorial))
surface_area = 2.0 * pi * equatorial * equatorial * (
1.0 + (polar * polar / (equatorial * equatorial)) * _series_atanh_over_x(eccentricity)
)
mean_curvature_integral = (
2.0 * pi * polar
+ (2.0 * pi * equatorial / eccentricity) * math.atan((equatorial * eccentricity) / polar)
)
else:
eccentricity = math.sqrt(1.0 - (equatorial * equatorial) / (polar * polar))
surface_area = 2.0 * pi * equatorial * equatorial * (
1.0 + (polar / equatorial) * _series_asin_over_x(eccentricity)
)
mean_curvature_integral = (
2.0 * pi * polar
+ (2.0 * pi * equatorial * equatorial / (polar * eccentricity)) * math.atanh(eccentricity)
)
return surface_area * t + mean_curvature_integral * t * t + (4.0 * pi / 3.0) * t**3
def verify_samples() -> None:
"""Check the two sample values given in the problem statement."""
sphere_case = chocolate_volume(1.0, 1.0)
expected_sphere = 28.0 * math.pi / 3.0
if abs(sphere_case - expected_sphere) > 1e-10:
raise AssertionError(f"sphere sample failed: {sphere_case} != {expected_sphere}")
oblate_case = chocolate_volume(2.0, 1.0)
if abs(oblate_case - 60.35475635) > 1e-8:
raise AssertionError(f"a=2,b=1 sample failed: {oblate_case}")
def solve() -> str:
"""Compute the answer for a=3 mm and b=1 mm, rounded to 8 decimal places."""
verify_samples()
return f"{chocolate_volume(3.0, 1.0):.8f}"
def main() -> None:
"""Run sample checks, solve the target case, and report timing."""
start = perf_counter()
verify_samples()
answer = solve()
elapsed = perf_counter() - start
print(answer)
print(f"elapsed_seconds={elapsed:.6f}")
if __name__ == "__main__":
main()
模块与导入
from __future__ import annotations：启用延迟注解，便于类型提示。import math：提供 \(\pi\)、反双曲函数、反三角函数等。from time import perf_counter：高精度计时，用于性能报告。
_series_atanh_over_x
- 计算 \(\operatorname{artanh}(x)/x\)。
- 当 \(|x|<10^{-8}\) 时用泰勒级数 \(1+x^2/3+x^4/5+\cdots\)，避免 \(x\to 0\) 时的数值不稳定。
- 否则直接调用
math.atanh(x)/x。
_series_asin_over_x
- 计算 \(\arcsin(x)/x\)。
- 小 \(|x|\) 时用级数 \(1+x^2/6+3x^4/40+\cdots\)。
- 否则调用
math.asin(x)/x。
chocolate_volume
- 参数
a,b对应题面方程中的 \(a,b\)；thickness为糖衣厚度（默认 1 mm）。 - 先校验参数为正、厚度非负。
- 球体分支：若 \(|a-b|<10^{-14}\)，用 \(S=4\pi a^2\)、\(M=4\pi a\)，代入 Steiner 公式。
- 一般分支：令赤道半径
equatorial=a，极半径polar=b。 - 扁椭球（
equatorial > polar）：算离心率 \(e\)，用扁椭球表面积与 \(M\) 公式，再组合 Steiner 三项。 - 长椭球（否则）：用长椭球对应公式。
- 返回值即糖衣体积 \(S t + M t^2 + \dfrac{4\pi}{3}t^3\)。
verify_samples
- 验证 \(a=b=1\) 时结果为 \(\dfrac{28}{3}\pi\)。
- 验证 \(a=2,b=1\) 时结果约为 \(60.35475635\)（容差 \(10^{-8}\)）。
solve
- 先调用
verify_samples()确保样例通过。 - 计算
chocolate_volume(3.0, 1.0)并格式化为 8 位小数字符串。
main
- 记录起始时间，执行样例校验与求解，打印答案与耗时。