Fortran中MacLaurin展开实现正弦函数的微小精度差异问题

首先得说，你的程序逻辑是完全正确的——递推计算每一项的方式没问题，项数也卡对了：输入30项时，确实覆盖到了x^59/59!这一项，循环次数和初始项的累加也没毛病。那结果的细微差别主要来自这几个点：

• 级数近似的固有误差：MacLaurin级数是sin(x)的无穷级数展开，取前30项只是个近似值，而教授给出的输出应该是Fortran内置sin()函数的结果。内置函数用的是更高效、精度拉满的算法（比如CORDIC或者优化的多项式逼近），能达到double precision的理论精度极限，和级数近似的结果自然会有最后几位的差异。

• 浮点运算的累积舍入误差：30次递推计算里，每次乘法、除法都会引入一点点舍入误差，累积起来就可能影响到结果的最后几位。你可以试试调整递推的计算顺序，比如把除法提前：

  temp = temp * (-x*x) / (i*(i-1))
  

  虽然数学上和你当前的写法等价，但浮点运算的顺序变化可能会让舍入误差更小一点。

• 编译器的浮点模式差异：不同编译器或者编译选项会影响浮点计算的精度，比如有些编译器默认用80位扩展精度寄存器算中间结果，再截断到64位double存储，这和严格64位计算的结果会有细微差别。比如用gfortran的话，可以试试加-ffloat-store选项强制用64位精度计算。

如果教授要求输出必须和内置函数完全一致，那直接调用sin(x)是最直接的办法；如果必须用MacLaurin级数，你可以试试把项数加到40左右，这样近似结果会更接近精确值，差异会缩小到可以忽略的程度。

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