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

我正在用 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
                ^^^^

我无法弄清楚错误出在哪里。

我无法弄清楚错误出在哪里。

精度高sine(5.0) is -0.95892427466313846889...

OP的结果比Prof的好。

OP的结果在14以内ULP https://en.wikipedia.org/wiki/Unit_in_the_last_place最佳答案，而教授的则有 25 ULP 折扣。

所以OP没有问题。为了得到一个exact与教授的答案相匹配，您将不得不编写一种较差的方法。

教授的回答较差的一个简单的可能原因是教授的代码是否仅在 while 时循环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)

Notes:
Near x == 5.0，大约第 17 个任期后，temp项很小，不会对结果产生显着影响。所以当然使用了足够多的术语。

使用常见的 FP 表示法，5 的最后一位单位是 ~89e-17。如果是 ~11e-17，则 ULP 为 -0.9589。