正弦餘弦模塊化擴展精度算術

Question

我在很多正弦/餘弦實現中見過所謂的擴展模塊化精度算法。但它是爲了什麼？ 例如在cephes implemetation中，在縮小到[0，pi/4]範圍後，他們正在進行模塊化精度算法以提高精度。

下文的代碼：

z = ((x - y * DP1) - y * DP2) - y * DP3; 

其中DP1，DP2和DP3有一些硬編碼的係數。 如何從數學上找到這些係數？我已經理解了大數字的「模塊化擴展算術」的目的，但是這裏的確切目的是什麼？

Answer 1

在減少三角函數的參數的情況下，您所看到的是Cody-Waite參數的減少，這是一本書中介紹的技巧：William J. Cody和William Waite，軟件手冊， Prentice-Hall，1980。儘管在中間計算中使用了subtractive cancellation，但其目的是爲了達到一定數量的參數而實現準確的減少的參數。爲此目的，通過使用多個遞減幅度的和（這裏：DP1,DP2,DP3），使得除了最不重要的中間產品之外的所有中間產品都可以通過以上的本地精度來表示相關常數沒有舍入誤差的情況下計算。

考慮以IEEE-754 binary32（單精度）計算sin（113）爲例。典型的參數減少將在概念上計算i=rintf(x/(π/2)); reduced_x = x-i*(π/2)。 binary32最接近π/ 2的數字是0x1.921fb6p+0。我們計算i=72，產品輪到0x1.c463acp+6，這與參數x=0x1.c40000p+6很接近。在減法期間，一些前導位取消，我們用reduced_x = -0x1.8eb000p-4結束。請注意重正化引入的尾隨零。這些零位沒有有用的信息。對簡化參數應用精確近似值sin(x) = -0x1.8e0eeap-4，而真實結果爲-0x1.8e0e9d39...p-4。我們結束了很大的相對誤差和大的ulp錯誤。

我們可以通過使用兩步Cody-Waite參數縮減來解決此問題。例如，我們可以使用pio2_hi = 0x1.921f00p+0和pio2_lo = 0x1.6a8886p-17。注意在pio2_hi的單精度表示中的八個尾隨零位，它允許我們乘以任何8位整數i並且仍然具有產品i * pio2_hi可表示爲正好作爲單精度數字。當我們計算((x - i * pio2_hi) - i * pio2_lo)時，我們得到reduced_x = -0x1.8eafb4p-4，因此sin(x) = -0x1.8e0e9ep-4，這是一個非常準確的結果。

將常量拆分爲總和的最佳方法取決於我們需要處理的i的大小，對於給定參數範圍的最大位數受到減法消除的影響（基於π的整數倍數/ 2可以達到整數）和性能考慮。典型的現實生活用例涉及2至4個階段的Cody-Waite簡化方案。融合多重加法（FMA）的可用性允許使用具有較少尾隨零位的組分常量。請參閱本文：Sylvie Boldo，Marc Daumas和Ren-Cang Li，「使用融合乘法加法正式驗證了參數的減少。」 IEEE Transactions on Computers,58：1139-1145,2009。對於使用fmaf()的工作示例，您可能想要查看one of my previous answers中的代碼。