寻求|x|<π/4范围内tan(x)的双精度1 ULP逼近优化方案

需求说明 #

• 目标逼近形式：tan(x) = x * P(x²) / Q(x²)，其中P、Q均为x²的三次多项式
• 约束条件：
  • 必须满足x=0处函数值与梯度匹配（即P₀=Q₀，保证极小x时tan(x)≈x）
  • 理想情况满足x=π/4处匹配，目前可放宽此处梯度约束
  • 优先保留某一x⁶项系数为1.0以节省乘法运算，若为达到16位有效十进制精度可放宽此要求

当前误差表现 #

通过Pade逼近结合模拟退火得到的系数，在x≈0.14处存在相对误差尖峰（约5.5e-16），十亿次随机测试中约60个样本超出阈值，误差直方图如下：

-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逼近系数 #

十进制系数Remez逼近系数 #

高精度十六进制Remez*逼近系数 #

示例实现代码 #

// 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;
  }
}

优化难点与求助方向 #

• 瓶颈：最终优化阶段存在大量局部最优解，难以收敛到1 ULP精度
• 潜在方案：LLL方法可能有效，但暂无实现能力
• 寻求支持：
  • Intel LLL优化包、Sollya、ARM Remez包的使用指导
  • 更优的逼近方案设计
  • 启发式优化建议

内容的提问来源于stack exchange，提问作者Martin Brown