我想在范围为-pi/4到pi/4的区间内，快速地将tan(x)近似到1 ULP以内。我已经找到了一个几乎足够好的解决方案，但最后的两倍因子仍然让我无法理解，即使我花了相当多的计算时间和模拟退火算法。
我寻求一个特定的数值近似tan(x)的方法，其形式为：
tan(x) = x*P(x^2)/Q(x^2)
其中P和Q都是x^2的三次方程。我需要对近似方法施加一些额外的约束条件。即在零点处以及理想情况下在pi/4处匹配函数值和梯度。目前，我已经决定放宽在pi/4处的梯度约束。我在这里应该说明的是，我寻求的近似方法在Hart的“计算机近似圣经”中是不可思议地缺失，这表明它是一个相当棘手的收敛问题。
通常他会给出整个序列，但是tan和pi/4的表格跳过了3,3，这很奇怪，因为它应该对于64位实数是完美的。精度是在P和Q的最高非零系数为N和M时获得的正确小数位数。
我所拥有的系数是通过模拟退火和一些本地技巧从Pade逼近中得出的。虽然还可以，但对于我的目的来说还不够好。有一个顽固的区间约为0.14，其中相对误差峰值保持在5.5e-16（大约60亿个测试案例中的每个）并且4.4e-16区间中的数值稍微过多。这是我在10^9个随机挑战测试中迄今为止所得到的相对误差直方图。
熵是比较直方图分布的一个快速评估指标。对于所有值都在一个区间的完美解决方案，熵将为零。
各种近似的系数如下所示：
Pade 0 (1) 1 (x^2) 2 (x^4) 3 (x^6)
Pn 135135 -17325 378 -1
Qn 135135 -62370 3150 28
Pade近似的最大相对误差为1.34e-4，发生在pi/4处。
具有小数系数的Remez在5.5e-16的精度下表现良好，足够快速查看，并且在Excel电子表格中可以正常工作，但由于小数舍入而不精确。
由于特定原因，我需要保持P0=Q0，以便对于非常小的值x，x == tan(x)。我还希望将x^6项的系数之一保持为1.0（节省一次乘法），但如果为了获得16位有效小数而必须做出妥协，那就是它。
这是一个示例实现（使用较短的、稍微不太准确的小数值）。升级为MRE 24/11/23。
优化之所以困难，是因为在最后阶段有一个平滑的函数，但上面有一片稀疏的尖钉，所以有无数个局部最优解，它们都差不多，只有很少几个好的解。密码学中的LLL方法被认为可以解决这个难题，但我个人的编程能力无法从零开始实现。问题归结为找到一组尾数值，它们的比率具有“神奇”属性，使得它们的比率比任何单个系数更精确。
如果有人能够使用Intel的LLL优化包、Sollya或ARM的Remez包，或者能找到稍微更好的近似解，我会非常感兴趣。我自己尝试过使用ARM的代码和Sollya，但效果不如预期。它离所需的精度非常接近，但还差一点。我也对其他启发式调整的建议很感兴趣，我有一些想法。
我寻求一个特定的数值近似tan(x)的方法，其形式为：
tan(x) = x*P(x^2)/Q(x^2)
其中P和Q都是x^2的三次方程。我需要对近似方法施加一些额外的约束条件。即在零点处以及理想情况下在pi/4处匹配函数值和梯度。目前，我已经决定放宽在pi/4处的梯度约束。我在这里应该说明的是，我寻求的近似方法在Hart的“计算机近似圣经”中是不可思议地缺失，这表明它是一个相当棘手的收敛问题。
通常他会给出整个序列，但是tan和pi/4的表格跳过了3,3，这很奇怪，因为它应该对于64位实数是完美的。精度是在P和Q的最高非零系数为N和M时获得的正确小数位数。
| Precision | N | M | Index |
|---|---|---|---|
| 10.66 | 2 | 2 | 4283 |
| 13.62 | 2 | 3 | 4284 |
| 19.74 | 4 | 3 | 4285 |
| 19.94 | 3 | 4 | 4286 |
-5.5511e-16 | 0
-4.4409e-16 | 704272
-3.3307e-16 | 0
-2.2204e-16 | 157067630
-1.1102e-16 | 0
0.0000e+00 | 485332167
1.1102e-16 | 240168741
2.2204e-16 | 99172390
3.3307e-16 | 16874505
4.4409e-16 | 680236
5.5511e-16 | 59
6.6613e-16 | 0
Mean 1.94185e-17 StdDev 1.31764e-16 Entropy 1.29231
熵是比较直方图分布的一个快速评估指标。对于所有值都在一个区间的完美解决方案，熵将为零。
各种近似的系数如下所示：
Pade 0 (1) 1 (x^2) 2 (x^4) 3 (x^6)
Pn 135135 -17325 378 -1
Qn 135135 -62370 3150 28
Pade近似的最大相对误差为1.34e-4，发生在pi/4处。
具有小数系数的Remez在5.5e-16的精度下表现良好，足够快速查看，并且在Excel电子表格中可以正常工作，但由于小数舍入而不精确。
| Remez | 0 (1) | 1 (x^2) | 2 (x^4) | 3 (x^6) |
|---|---|---|---|---|
| Pn | 133506.920100216 | -17127.5484486707 | 374.700589646508 | -1 |
| Qn | 133506.920100216 | -61629.8551487417 | 3117.06295917985 | -27.8423369154197 |
Remez完全准确的十六进制浮点表示（迄今为止最好的）
| Remez* | 0 (1) | 1 (x^2) | 2 (x^4) | 3 (x^6) |
|---|---|---|---|---|
| Pn | 0x1.04c175c5d80a5p+17 | -0x1.0b9e319c87430p+14 | 0x1.76b359d7d3ac0p+8 | -1 |
| Qn | 0x1.04c175c5d80a5p+17 | -0x1.e17bb5d60e518p+15 | 0x1.85a203c2f84f8p+11 | -0x1.bd7a3645ff105p+4 |
这是一个示例实现（使用较短的、稍微不太准确的小数值）。升级为MRE 24/11/23。
// Toy_Pade.cpp : Demonstrates the failing cases err>5.5e-16 and prints out failures
//
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <math.h>
double TanPQ(double x)
{
const double P[4] = { 133506.920100216, -17127.548448670765, 374.700589646508, -1.0 };
const double Q[4] = { 133506.920100216, -61629.8551487418, 3117.06295917985, -27.84233691541965 };
double sumP, sumQ, x2;
int i;
sumP = sumQ = 0;
x2 = x * x;
for (i = 3; i >= 0; i--)
{
sumP = sumP * x2 + P[i];
sumQ = sumQ * x2 + Q[i];
}
return x * sumP / sumQ;
}
double tan87_diff(double x)
// return x86 fptan x forced to the x87 coprocessor or tanl()
{
long double y,t;
if (x == 0) return 0;
t = TanPQ(x);
#ifdef _M_IX86
_asm {
fld qword ptr[x]
fptan
fxch st(1) // save the 1
fdivr qword ptr[t]
fsubp st(1), st
fstp qword ptr[y]
}
#else
y = 1 - t / tanl(x);
#endif
return y;
}
int main()
{
double x, dx, err;
dx = M_PI / 2000000000;
x = 0;
while (x < M_PI / 4)
{
err = tan87_diff(x);
if (fabs(err) > 5.0e-16) printf("%-26.18g %22.14a %g\n", x, x, err);
x = x + dx;
}
}
优化之所以困难，是因为在最后阶段有一个平滑的函数，但上面有一片稀疏的尖钉，所以有无数个局部最优解，它们都差不多，只有很少几个好的解。密码学中的LLL方法被认为可以解决这个难题，但我个人的编程能力无法从零开始实现。问题归结为找到一组尾数值，它们的比率具有“神奇”属性，使得它们的比率比任何单个系数更精确。
如果有人能够使用Intel的LLL优化包、Sollya或ARM的Remez包，或者能找到稍微更好的近似解，我会非常感兴趣。我自己尝试过使用ARM的代码和Sollya，但效果不如预期。它离所需的精度非常接近，但还差一点。我也对其他启发式调整的建议很感兴趣，我有一些想法。
d = 133506.920100216; printf("%.17g %a\n", d,d);-->133506.920100216 0x1.04c175c5d8086p+17，这与0x1.04c175c5d80a5p+17有所不同。 - undefinedfma()作为解决方案的可能性是否存在？ - undefined