floating-point - 使用 MacLaurin 展开的 Fortran Sine 函数的微小差异

Question

我正在 Fortran 中创建一个程序，它以弧度表示 sin(x) 的 x，然后是要计算的项数。

这是我的程序:

! Sine value using MacLaurin series 

program SineApprox
    implicit none
    integer :: numTerms, denom, i
    double precision :: x, temp, summ

    ! Read Angle in Radians, and numTerms
    read(*,*) x, numTerms

    ! Computing MacLaurin series
    denom = (2*numTerms)-1  !Gives denominator
    temp = x                !Temp calculates each term
    summ = x

    do i = 3, denom, 2
        temp = (((-temp)*(x**2))/(i*(i-1)))
        summ = summ + temp 
    end do
    print *, summ

end program SineApprox

但是，我没有得到我的教授要求输入的相同值:5 30

我的代码的输出是:

-0.95892427466314001  

但是，所需的输出是:

-0.95892427466313568
                ^^^^

我不知道错误在哪里。

Answer 1

I can't figure out where the error is.

一个高精度的sine(5.0)是-0.95892427466313846889...

OP 的结果比 Prof 的要好。

OP 的结果与最佳答案相差 14 ULP，而 Prof 的结果与最佳答案相差 25 ULP。

所以对 OP 来说没有问题。要获得与教授的答案完全匹配的结果，您必须编写一种劣等方法。

教授较差答案的一个简单可能原因是，如果教授的代码仅在 i<=30(15 个术语，而不是 30 个术语)时循环 - 这几乎可以解释差异。尝试使用迭代次数较少的代码，看看哪些迭代次数最接近教授的答案。

---

//  sine(5.0)             Delta from correct  Source
//  
//-0.95892427466313846889...    0 000000      calc
//-0.95892427466313823        -23 889...      chux (via some C code)
//-0.95892427466314001       +154 111...      OP
//-0.95892427466313568       -278 889...      Prof
// 0.00000000000000089      +/-89 ...         ulp(5.0)
// 0.00000000000000011      +/-11 ...         ulp(-0.9589)
//   12345678901234567(digit count)

---

注意事项:
在 x == 5.0 附近，大约在第 17 项之后，temp 项非常小，不会显着影响结果。所以肯定使用了足够多的术语。

对于普通的FP表示，5的最后一位的单位是~89e-17。如果是 ~11e-17，则 -0.9589 的 ULP。