你遇到的核心问题是灾难性数值抵消：当|x|非常小时，cosx和cos2x都极其接近1，直接相减会丢失几乎所有有效精度，同时因为z趋近于0，直接计算条件数时确实会出现分母趋近于0导致数值爆炸的情况。解决这个问题的关键是通过三角恒等变换，把原式转换成没有抵消风险的形式，同时能让条件数保持在可控范围。

第一步：用三角恒等式重构表达式 #

回忆余弦差的恒等式：

cosA - cosB = -2sin[(A+B)/2]sin[(A-B)/2]

把A=x、B=2x代入，化简后得到：

cosx - cos2x = 2*sin(3x/2)*sin(x/2)

这个变换是精确的，没有任何近似，完全等价于原式。

第二步：分析变换后的稳定性与条件数 #

当|x|≪1时，sinθ≈θ - θ³/6 + θ⁵/120（泰勒展开），所以变换后的表达式可以近似为：

z ≈ 2*(3x/2)*(x/2) = (3/2)x²

这时候我们来看相对条件数（更能反映数值稳定性的指标）：
相对条件数定义为 K_rel = |(dz/z)/(dx/x)| = |(x/z)*(dz/dx)|

对变换后的z求导：

dz/dx = 3*cos(3x/2)*sin(x/2) + sin(3x/2)*cos(x/2)

当x→0时，cosθ≈1，sinθ≈θ，代入后：

dz/dx ≈ 3*(x/2) + (3x/2)*1 = 3x

再代入相对条件数公式：

K_rel ≈ |(x/(3x²/2))*3x| = 2

看到了吗？这个条件数是一个常数，不会随着x趋近于0而变大，完美解决了你担心的条件数爆炸问题。

第三步：实现稳定算法的具体步骤 #

当|x|≪1时，推荐这样的计算流程：

• 计算两个角度：theta1 = 3*x/2，theta2 = x/2
• 计算这两个角度的正弦值：可以直接调用编程语言的标准库sin()函数（现代标准库的sin函数在0附近的计算已经做了优化，精度很高）；如果x极小到需要更高精度，也可以用泰勒展开近似：sin(theta) ≈ theta - theta^3/6 + theta^5/120
• 最终结果：z = 2 * sin(theta1) * sin(theta2)

为什么这个算法稳定？ #

直接计算cosx - cos2x时，两个接近1的浮点数相减，会丢失大部分有效位（比如x=1e-4时，cosx≈1-5e-9，cos2x≈1-2e-8，单精度浮点数下这两个值都会被存储为1.0，相减结果为0，完全错误）。而变换后的表达式是两个小量的乘积，不会出现抵消问题，计算出来的z能保留足够的有效精度。

备注：内容来源于stack exchange，提问作者J3ck_Budl7y