基于二次多项式的C语言快速指数实现：负输入结果异常问题

我基于二次多项式实现了指数函数e^x的快速近似，但当输入x<0.0时，程序无法输出合理结果，结果存在局部最小值等错误。

原实现代码 #

// Compile with: gcc -std=c99 fastexp.c -o fastexp

#include <stdio.h>
#include <math.h>

#define ALPHA 0.50617024f
#define ALPHAc 0.49382976f
#define INV_LN2 1.4426950408889634f

float fast2pxm1_m(float x) { // 0.0 <= x <= 1.0
  return x*(ALPHA+ALPHAc*x);
}

float fast2px(float x) {
  float m, mp;
  float e;

  m = modff(x, &e);
  mp = fast2pxm1_m(m) + 1.0f;
  
  return ldexpf(mp, e);
}

float fastexp(float x) {  // exportable
  return fast2px(x*INV_LN2);
}

int main(int argc, char *argv[]) {
  float x;

  for (int i = -10; i <= 10; i++) {
     x= (float) i/10.0*2;
     printf("%.8f => %.8f\n",x,fastexp(x));
  }
  return 0;
}

错误运行结果 #

-2.00000000 => 0.23474069
-1.79999995 => 0.21845233
-1.60000002 => 0.22272082
-1.39999998 => 0.24754614
-1.20000005 => 0.44696173
-1.00000000 => 0.43635058
-0.80000001 => 0.46685311
-0.60000002 => 0.93187356
-0.40000001 => 0.87235498
-0.20000000 => 0.89506382
0.00000000 => 1.00000000
0.20000000 => 1.18716359
0.40000001 => 1.45655441
0.60000002 => 1.80817270
0.80000001 => 2.17952919
1.00000000 => 2.64171839
1.20000005 => 3.26836252
1.39999998 => 4.04080629
1.60000002 => 4.81200027
1.79999995 => 5.91210270
2.00000000 => 7.34111547

问题分析 #

错误根源在于fast2px函数中使用modff拆分整数和小数部分的逻辑。modff返回的小数部分符号与原数一致，当x<0时，小数部分m为负数，超出了fast2pxm1_m函数仅处理[0,1)区间输入的设计范围，导致二次多项式计算出错误的2^m-1近似值，最终整个e^x的结果失真。

比如输入x=-0.6时，x*INV_LN2≈-0.8656，modff返回m=-0.8656、e=0，代入fast2pxm1_m得到错误的负值，加1后得到0.9319，再经ldexpf缩放后输出错误的0.9318，而正确的e^-0.6≈0.5488。

修复方案 #

修改fast2px函数的拆分逻辑，使用floorf将任意实数x拆分为整数部分k（可正可负）和小数部分f（始终落在[0,1)区间），确保fast2pxm1_m的输入符合要求：

float fast2px(float x) {
  float k = floorf(x);
  float f = x - k;
  float mp = fast2pxm1_m(f) + 1.0f;
  return ldexpf(mp, (int)k);
}

另外，替换为补充的优化常数可以进一步提升近似精度：

#define ALPHA 0.6563658616549711f               
#define ALPHAc 0.34363413834502887f             

修复后的程序运行结果会符合预期，负输入的e^x近似值将正确递减，不再出现局部最小值。

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