
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
Author 
Message 
mecej4
Joined: 31 Oct 2006 Posts: 1072

Posted: Sun Apr 02, 2017 9:55 pm Post subject: 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)!
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.080298E02
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 


LitusSaxonicum
Joined: 23 Aug 2005 Posts: 1946 Location: Yateley, Hants, UK

Posted: Sun Apr 02, 2017 11:40 pm Post subject: 


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 hardwired 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 8noded numericallyintegrated isoparametric finite elements, and turned me off spending time on those mathematical solutions. 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1072

Posted: Mon Apr 03, 2017 1:17 am Post subject: 


The standard puts the details of the calculation of mathematical functions under the "processordependent" 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 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2053 Location: Sydney

Posted: Mon Apr 03, 2017 7:36 am Post subject: 


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 


PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 5760 Location: Salford, UK

Posted: Mon Apr 03, 2017 8:00 am Post subject: 


Thank you. I have made a note of this. 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 1072

Posted: Mon Apr 10, 2017 8:24 pm Post subject: Re: 


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 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2053 Location: Sydney

Posted: Tue Apr 11, 2017 1:36 pm Post subject: 


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 


LitusSaxonicum
Joined: 23 Aug 2005 Posts: 1946 Location: Yateley, Hants, UK

Posted: Tue Apr 11, 2017 4:28 pm Post subject: 


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 


mecej4
Joined: 31 Oct 2006 Posts: 1072

Posted: Wed Apr 12, 2017 4:41 am Post subject: 


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 highprecision 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(1E18int(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 inmemory representations of REAL*4 and REAL*8 variables. In the X87, with its 80bit 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 


mecej4
Joined: 31 Oct 2006 Posts: 1072

Posted: Wed Apr 12, 2017 4:53 am Post subject: Re: 


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 


LitusSaxonicum
Joined: 23 Aug 2005 Posts: 1946 Location: Yateley, Hants, UK

Posted: Wed Apr 12, 2017 12:48 pm Post subject: 


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 


mecej4
Joined: 31 Oct 2006 Posts: 1072

Posted: Wed Apr 12, 2017 1:01 pm Post subject: 


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 coprocessor (the brain of the user). A slide rule is still faster for solving oneoff 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 


LitusSaxonicum
Joined: 23 Aug 2005 Posts: 1946 Location: Yateley, Hants, UK

Posted: Wed Apr 12, 2017 1:51 pm Post subject: 


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 threedimensional finite element model employing viscoelastoplastic 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 prewar Austin 7s, which at one time could be had for £3 to £5.) 

Back to top 




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
