forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Bug in evaluation of sin(x) for large x

 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support
View previous topic :: View next topic  
Author Message
mecej4



Joined: 31 Oct 2006
Posts: 1884

PostPosted: Sun Apr 02, 2017 9:55 pm    Post subject: Bug in evaluation of sin(x) for large x Reply with quote

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)!

Code:

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:
Code:

           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
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Sun Apr 02, 2017 11:40 pm    Post subject: Reply with quote

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.
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1884

PostPosted: Mon Apr 03, 2017 1:17 am    Post subject: Reply with quote

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.
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2551
Location: Sydney

PostPosted: Mon Apr 03, 2017 7:36 am    Post subject: Reply with quote

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.
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 7914
Location: Salford, UK

PostPosted: Mon Apr 03, 2017 8:00 am    Post subject: Reply with quote

Thank you. I have made a note of this.
Back to top
View user's profile Send private message AIM Address
mecej4



Joined: 31 Oct 2006
Posts: 1884

PostPosted: Mon Apr 10, 2017 8:24 pm    Post subject: Re: Reply with quote

JohnCampbell wrote:
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.
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2551
Location: Sydney

PostPosted: Tue Apr 11, 2017 1:36 pm    Post subject: Reply with quote

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

Code:
!  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
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Tue Apr 11, 2017 4:28 pm    Post subject: Reply with quote

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.
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1884

PostPosted: Wed Apr 12, 2017 4:41 am    Post subject: Reply with quote

Quote:
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))*2*pi), 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 REAL*4 and REAL*8 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.


Last edited by mecej4 on Wed Apr 12, 2017 12:49 pm; edited 2 times in total
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1884

PostPosted: Wed Apr 12, 2017 4:53 am    Post subject: Re: Reply with quote

LitusSaxonicum wrote:
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.


Last edited by mecej4 on Mon Sep 18, 2017 1:15 am; edited 1 time in total
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Wed Apr 12, 2017 12:48 pm    Post subject: Reply with quote

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
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1884

PostPosted: Wed Apr 12, 2017 1:01 pm    Post subject: Reply with quote

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!


Last edited by mecej4 on Wed Apr 12, 2017 2:02 pm; edited 2 times in total
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Wed Apr 12, 2017 1:51 pm    Post subject: Reply with quote

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.)
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group