Bug in evaluation of sin(x) for large x

FTN95 fails to do argument reduction when sin(x) is called with values of x >= 1E19 and returns the argument as the value of sin(x)!

program sinerr
implicit none
real theta,s
integer i
theta=1e10
print *,'          i       theta          sin(theta)'
do i=1,15
   s = sin(theta)
   print *,i,theta,s
   theta=theta*10
end do
end program

The output:

           i       theta          sin(theta)
            1    1.000000E+10   -0.487506
            2    1.000000E+11    0.998109
            3    1.000000E+12   -2.080298E-02
            4    1.000000E+13    0.968876
            5    1.000000E+14    0.967149
            6    1.000000E+15    0.994434
            7    1.000000E+16   -0.496563
            8    1.000000E+17   -0.570077
            9    1.000000E+18   -0.215481
           10    1.000000E+19    1.000000E+19
           11    1.000000E+20    1.000000E+20
           12    1.000000E+21    1.000000E+21
           13    1.000000E+22    1.000000E+22
           14    1.000000E+23    1.000000E+23
           15    1.000000E+24    1.000000E+24

Is argument reduction required by the standard? I remember being surprised by Microsoft Fortran losing accuracy in trig functions at high arguments when I was evaluating a Fourier series over 20 years ago. I can't remember where I read up on it, probably '8087 applications and programming for the IBM PC and other PCs' by Richard Startz, althoughI don't have a copy to hand to check. The hard-wired routines in an x87 have a limited (although still huge) range.

I formed the view, correctly or not, that parameter reduction was the job of the programmer. I suspect I am about to learn differently!

Eddie

The Fourier solution turned out to take a great deal more computation than solving the same problem to the same accuracy using a small number of 8-noded numerically-integrated isoparametric finite elements, and turned me off spending time on those mathematical solutions.

The standard puts the details of the calculation of mathematical functions under the 'processor-dependent' category.

The X87 had the FPREM (and, in later models, the FPREM1) instruction for the purpose of range reduction. What is more, after range reduction two flag bits were set so that one could find the signs of the sine and cosine of the argument.

Providing any answer for SIN(1.e8) or SIN(1.d18) is debatable. Doesn't your example use real*4 values ? At least the value supplied could indicate a gross error.

Thank you. I have made a note of this.

Quoted from JohnCampbell Providing any answer for SIN(1.e8) or SIN(1.d18) is debatable.

Why? The number (1E8, i.e., 100 million) is exactly representable in 32 or 64 bit IEEE format. Old Lady Sine traverses a graceful path between two straight insurmountable walls, -1 on the left and +1 on the right.

If any answer is given all, that answer should lie between -1 and +1, limits included.

mecej4,

I may not have explained my comment well enough, but the value of sin(1.e8) or sin(1.d18) is going to be rubbish, as you can not define the angle to any reasonable precision.

You say 'The number (1E8, i.e., 100 million) is exactly representable in 32 or 64 bit IEEE format' but I don't agree as I don't think it can be used to that accuracy. You can't claim that the computation will be based on a number to any accuracy that could be exact, say =/- .0000001.

The following program demonstrates (to me) that 1.e8 is accurate to +/- 8 and 1.d18 is accurate to +/- 128. I would suggest that the sine of the angle being calculated is an unknown angle to +/- 8 or 128 radians, which would return a random value of sin (x). What good is getting a sin value in the range -1:1, if the angle has no accuracy to a fraction of 2_pi.

John

!  test of precision using spacing and nearest intrinsic
!
      real*4 :: x4 = 1.e8
      real*8 :: x8 = 1.e18
!
      write (*,*) 'x       =', x4
      write (*,*) 'sin     =', sin(x4)
      write (*,*) 'spacing =', spacing (x4)       
      write (*,*) 'nearest =', nearest (x4, -1.0)
      write (*,*) nearest (x4, 1.0)-nearest (x4, -1.0)
!
      write (*,*) ' '
      write (*,*) 'x       =', x8
      write (*,*) 'sin     =', sin(x8)
      write (*,*) 'spacing =', spacing (x8)       
      write (*,*) 'nearest =', nearest (x8, -1.0d0)      
      write (*,*) nearest (x8, 1.0d0)-nearest (x8, -1.0d0)
!
      end

Aren't the trig functions computed with essentially a finite number of terms in a series, so that for beyond a certain value of x, sin(x) (for example) rapidly gets less accurate.

I think one should beware when stepping outside -2pi to + 2pi, and that means doing the reduction yourself.

John Campbell: I may not have explained my comment well enough, but the value of sin(1.e8) or sin(1.d18) is going to be rubbish, as you can not define the angle to any reasonable precision.

John, if you correct your program by writing '1D18' instead of '1E18', compile and run it, all the digits that it prints for sin(1E8) and sin(1D18) are correct, as you can verify by using a high-precision calculator/calculator program. That is certainly not rubbish.

A useful online calculator is http://keisan.casio.com/calculator . Set the mode to Real Rad, set the number of digits to 30, paste in the expression sin(1E18-int(1E18/(2*pi))2pi), and after a couple of seconds you will see -0.9929693207404, which is identical to the output of your program.

The SPACING, NEAREST, etc. that you used to reach your conclusions pertain to the in-memory representations of REAL4 and REAL8 variables. In the X87, with its 80-bit registers with 64 bits mantissae, if variables are stored in the FPU registers as much as possible, better results can be obtained, as your test program proves.

It is certainly possible to lose precision if care is not taken or the numbers are far larger in magnitude than, say, 1D18. Sometimes, it is sufficient to do only a part of a calculation in higher precision; for example, using a double precision register to accumulate the sum of squares of single precision vectors, which is a remedy that you are quite familiar with.

Quoted from LitusSaxonicum I think one should beware when stepping outside -2pi to + 2pi, and that means doing the reduction yourself.

Eddie, as you can see at http://x86.renejeschke.de/html/file_module_x86_id_114.html, you do not need to do range reduction yourself unless the argument of sin() falls outside the range [-2^{63, +2}63], rather than [-2 Pi, +2 Pi]. It is quite likely that range reduction is performed within the microcode of the X87 unit before a truncated power series or a Tchebycheff approximation is used.

Hi Mecej4,

I agree that you can go outside of the - to + 2pi range, my suggestion was to be careful when doing so.

If you keep all your calculations in terms of meaningful and visualizable quantities during a calculation, you can usually forget about all sorts of silliness that computers do. The more elaborate the calculation, the more often that silliness comes to the fore. Matrix solution is one of the areas where the intermediate results are not identifiable in terms of their physical meaning, and where things like roundoff catch you out. It's the underlying reason why aircraft and suspension bridges designed with slide rules worked.

As angles in the -2pi to +2pi range can be drawn with a protractor, then my simple rule of thumb tells me that computer functions ought to work as reliably as (say) the numbers in my log table book!

Eddie

Eddie, I remember my slide rule with fondness. I still have one. In contrast to doing calculating using tables, using the slide rule allowed one to develop and retain a sense of reasonable and unreasonable magnitudes, because the slide rule dealt only with normalized numbers, the exponents being processed in the co-processor (the brain of the user). A slide rule is still faster for solving one-off triangles using the sine rule than using any computer.

I remember a humorous situation in a laboratory. A young engineer was puzzled by an incorrect estimate. I pointed out that the plate thickness that she had measured in mm with a micrometer should have been expressed in meters for use in the formula. She fished out her calculator, entered the measured value, then 1000, and then pressed the / button!

It's worse than that. I work in a field where in some parts of it, the only numerical skills required were to multiply by 2, or divide by 10, the latter skill being changed after metrication to be able to multiply by 10.

Answers on a postcard as to what I am talking about, if anyone knows. I will post the answer in due course.

The sad thing is that the slide rule, then the pocket calculator and latterly the smartphone are commonly deployed to do the calculation.

There are people who think the answer from the above calculation is improved by doing a three-dimensional finite element model employing visco-elasto-plastic behaviour.

Eddie

(My slide rule is an early Thornton PIC, bought at the cost of 4 guineas. I had friends who bought cars for that or less, but admittedly they were rustbucket pre-war Austin 7s, which at one time could be had for £3 to £5.)