I wrote a program which has a bug when the value of u is greater than 1.4. when I tried to find the bug for values of y greater than 1.4 I saw an strange point :
u=0.D0
do while (u.lt.20.D0) . . . print*,u, y if (u==1.4D0) print*,"ok" . . . u=u+0.1D0 end do ------------------------ the output: 0.00000000000000 -0.900316316157106 0.100000000000000 -0.896557866215700 0.200000000000000 -0.892752298842229 0.300000000000000 -0.889009961947742 0.400000000000000 -0.885300844120634 0.500000000000000 -0.881608066131983 0.600000000000000 -0.877943572939470 0.700000000000000 -0.874308949879819 0.800000000000000 -0.870674660655724 0.900000000000000 -0.867096152952160 1.00000000000000 -0.863546597445969 1.10000000000000 -0.860024078786788 1.20000000000000 -0.856526008002337 1.30000000000000 -0.853103843918795 1.40000000000000 -0.849649287329888
but the program cannot print "ok" while it prints u=1.40000000000000 . why???
> I wrote a program which has a bug when the value of u is greater than > 1.4. > when I tried to find the bug for values of y greater than 1.4 I saw an > strange point :
> but the program cannot print "ok" while it prints u=1.40000000000000 . > why???
> I have got confused. Is there anyone can help me?
> Thanks
The detailed immediate problem is that 1.4d0 and 14.0d0 * 0.1d0 are different numbers. They do not compare equal as you have tested then.
The real error is that you do not understand the technical issues of floating point computation, otherwise know as finite precision computation. The problem is that 0.1d0 is a notational fiction as there is no such number in binary floating point. Ask google for "What every computer scientist should know about floating point arithmetic" by Goldberg as supplied by Sun. Excellent paper.
0.1 would be a nonterminating repeating fraction in binary as so is only approximated in finite precision.
One quick fix is to multiply every thing by 10 to get integers if you are only interested in multiples of 0.1. Depending upon what you real intent is there are other fixes.
> I wrote a program which has a bug when the value of u is greater than > 1.4. > when I tried to find the bug for values of y greater than 1.4 I saw an > strange point :
> u=0.D0
> do while (u.lt.20.D0) > . > . > . > print*,u, y > if (u==1.4D0) print*,"ok"
> but the program cannot print "ok" while it prints u=1.40000000000000 . > why???
> I have got confused. Is there anyone can help me?
This is probably the single most frequently asked question concerning computing in general: the short answer is "because computer's floating- point arithmetics is NOT exact, and not the same as mathematics real number arithmetics".
There are many similar tutorials on the web, and even several books on the subject. Please get familiar with the concept, if you plan to do engineering computing for living.
-- Jugoslav www.xeffort.com Please reply to the newsgroup. You can find my real e-mail on my home page above.
> I wrote a program which has a bug when the value of u is greater than > 1.4. > when I tried to find the bug for values of y greater than 1.4 I saw an > strange point :
> but the program cannot print "ok" while it prints u=1.40000000000000 . > why???
> I have got confused. Is there anyone can help me?
> Thanks
The classic problem of floating-point data: a finite binary representation is not capable of exactly representing 1.4, just a finite decimal representation can not represent 1/3 exactly.
You need to define a small margin, like:
if (abs(u-1.4D0) < 2.0*epsilon(u)*1.4d0) print*,"ok"
where epsilon() is a standard function returning the smallest number E such that 1 and 1+E can be distinguished.
It has everything to do with the finite precision of computers.
> I wrote a program which has a bug when the value of u is greater > than 1.4. > when I tried to find the bug for values of y greater than 1.4 I > saw an strange point :
> u=0.D0
> do while (u.lt.20.D0) > . > . > . > print*,u, y > if (u==1.4D0) print*,"ok" > . > . > . > u=u+0.1D0 > end do > ------------------------ > the output: <--CUT--> > but the program cannot print "ok" while it prints > u=1.40000000000000 . why???
> I have got confused. Is there anyone can help me?
> Thanks
Try this simpler example:
program testFloat real*8 X character*3 C data C/'1.4'/ read(C,*)X if(X.eq.1.4)then write(*,*)' Equal' else write(*,*)' Not Equal' endif end
Run it and see what it prints and compare to what you expect. Change the '1.4' on line-6 to '1.4d0' and try again.Then, do enough reading to make your expections more often in agreement with what the computer gives.
> do while (u.lt.20.D0) > . > . > . > print*,u, y > if (u==1.4D0) print*,"ok" > . > . > . > u=u+0.1D0 > end do
Some other posters have explained your original problem. I wanted to point out something else related to the floating point precision problem.
Your above loop structure is probably intended to do 200 passes. but if you actually count the passes, it might do 199, or 200, or 201 (depending on exactly how to test "u" to terminate the loop). This is related to the fact that 0.1D0 is a decimal number that is not represented exactly as a floating point value. The actual floating point number might be either a little larger or a little smaller than the decimal number. So when the value of "u" gets incremented, it will always be incremented with a number that is a little too large or too small, and those errors accumulate for every addition. After the value of "u" starts to build up, the increment operation is also subject to roundoff error, making the result of the addition either a little larger or a little smaller than it should be. This cumulative error will affect not only your calculation (because "u" is used elsewhere), but also the number of passes through the loop. So, in the end, the actual number of passes might be right (200), or it could be off a little. In your case, it might be off only by 1, most likely by being a little too large, but in other situations it might be off by even more than that.
So, what should you do to eliminate this cumulative error? If "u" is changed only in that one place, and you know the number of passes that are supposed to occur (i.e. in exact arithmetic), then you can write the loop as, for example,
do iter = 0, 199 u = (iter) * 0.1D0 ... enddo
(I suggest also that you parameterize your double precision constant with a real type, 0.1_dp, but this is a separate issue.)
With this loop, there is still a little error associated with each value of u, but it is only a local error (i.e. probably no more than the last bit only), it does not accumulate with the number of passes through the loop. Furthermore, you will always get the right number of passes because the test on the integer "iter" is always correct. Finally, one other benefit is that compilers can optimize this loop structure better than the do while form (e.g. with loop unrolling, use of registers, etc.).
Some of the above comments are personal style issues, others are practical issues that affect the results.
> I wrote a program which has a bug when the value of u is greater than > 1.4. > when I tried to find the bug for values of y greater than 1.4 I saw an > strange point :
> u=0.D0
> do while (u.lt.20.D0) > . > . > . > print*,u, y > if (u==1.4D0) print*,"ok" > . > . > . > u=u+0.1D0 > end do > but the program cannot print "ok" while it prints u=1.40000000000000 . > why???
1. Others have explained why comparing floating point numbers for equality is almost always a problem.
2. It is traditional to use an integer variable in a DO loop, then generate U or X or whatever from that.
For example
integer i real u
do i=0,200 u=i/10.0 ... if(i == 14) print *,'OK' ... end do
3. There are several advanatage to this approach.
a. calculating U = U0 + i * delta_U or U = U0 + (i-1) * delta_U avoids accumulating errors unlike
U=U0 do
U=U + delta U if ( ) exit end do
b. It is easy to construct an array as a look-up table of values for later use instead of re-computing them.
real aqz(100) integer i real x,x0,delta_x
do i=1,100 x=x0 + (i-1)*delta_x aqz(i) = <some complicated function of x> end do ....
given an X, finding which aqz()[s] to use for interpolation is easy!
Gordon Sande wrote: > The problem is that 0.1d0 is a notational fiction
I really like that phrase to refer to the issue of finite precision computation. Succinct and imbibed with some mild shock value to make one think about the underlying mechanism(s).
>> The problem is that 0.1d0 is a notational fiction
> I really like that phrase to refer to the issue of finite precision > computation. Succinct and imbibed with some mild shock value to make one > think about the underlying mechanism(s).
Gib Bogle wrote: > Paul van Delst wrote: >> Gordon Sande wrote:
>>> The problem is that 0.1d0 is a notational fiction
>> I really like that phrase to refer to the issue of finite precision >> computation. Succinct and imbibed with some mild shock value to make >> one think about the underlying mechanism(s).
> Imbued, unless you've been drinking. ;-)
Oops, yes. A, uh, Lagerian slip? I could go for a cold one right now. :o)
> On 1 apr, 16:08, Fatemeh <fateme.mirj...@gmail.com> wrote:
> > Dear all ;
> > I wrote a program which has a bug when the value of u is greater than > > 1.4. > > when I tried to find the bug for values of y greater than 1.4 I saw an > > strange point :
> > but the program cannot print "ok" while it prints u=1.40000000000000 . > > why???
> > I have got confused. Is there anyone can help me?
> > Thanks
> The classic problem of floating-point data: a finite binary > representation > is not capable of exactly representing 1.4, just a finite decimal > representation > can not represent 1/3 exactly.
> You need to define a small margin, like:
> if (abs(u-1.4D0) < 2.0*epsilon(u)*1.4d0) print*,"ok"
> where epsilon() is a standard function returning the smallest number E > such that 1 and 1+E can be distinguished.
> It has everything to do with the finite precision of computers.
> Regards,
> Arjen- Tekst uit oorspronkelijk bericht niet weergeven -
> - Tekst uit oorspronkelijk bericht weergeven -
Actually, it is easier than what I wrote:
if ( abs(u-1.4D0) < 0.05D0 ) ...
where 0.05 is half the step (0.1).
u can not assume a value further off the intended value than 1/2 the step, unless you accumulate a lot of rounding errors. (It is a trick that is especially useful if you need to compare u to zero: epsilon would be tricky to use in that region!)
And read all the other postings for a better understanding too ;)
> In article > <bd511e22-1b33-4f3c-b2b6-71c2f42af...@j39g2000yqn.googlegroups.com>, > Fatemeh <fateme.mirj...@gmail.com> wrote:
> > u=0.D0
> > do while (u.lt.20.D0) > > . > > print*,u, y > > if (u==1.4D0) print*,"ok" > > . > > u=u+0.1D0 > > end do
> Some other posters have explained your original problem. I wanted > to point out something else related to the floating point precision > problem.
> Your above loop structure is probably intended to do 200 passes. > but if you actually count the passes, it might do 199, or 200, or > 201 (depending on exactly how to test "u" to terminate the loop). > This is related to the fact that 0.1D0 is a decimal number that is > not represented exactly as a floating point value. The actual > floating point number might be either a little larger or a little > smaller than the decimal number. So when the value of "u" gets > incremented, it will always be incremented with a number that is a > little too large or too small, and those errors accumulate for every > addition. After the value of "u" starts to build up, the increment > operation is also subject to roundoff error, making the result of > the addition either a little larger or a little smaller than it > should be. This cumulative error will affect not only your > calculation (because "u" is used elsewhere), but also the number of > passes through the loop. So, in the end, the actual number of > passes might be right (200), or it could be off a little. In your > case, it might be off only by 1, most likely by being a little too > large, but in other situations it might be off by even more than > that.
> So, what should you do to eliminate this cumulative error? If "u" > is changed only in that one place, and you know the number of passes > that are supposed to occur (i.e. in exact arithmetic), then you can > write the loop as, for example,
> do iter = 0, 199 > u = (iter) * 0.1D0
Definitely not. This multiplies the error in 0.1 by iter, and as iter approaches 199, you have multiplied this error by 199. A sure way to minimise the error is to divide by 10. Thus, u = iter/10.0d0 has a small error, while iter*0.1d0 has a large error.
> ... > enddo
> (I suggest also that you parameterize your double precision constant > with a real type, 0.1_dp, but this is a separate issue.)
> With this loop, there is still a little error associated with each > value of u, but it is only a local error (i.e. probably no more than > the last bit only),
> it does not accumulate with the number of passes > through the loop. Furthermore, you will always get the right number > of passes because the test on the integer "iter" is always correct. > Finally, one other benefit is that compilers can optimize this loop > structure better than the do while form (e.g. with loop unrolling, > use of registers, etc.).
In article <%k3Bl.397$y61....@news-server.bigpond.net.au>,
"robin" <robi...@bigpond.com> wrote: > > So, what should you do to eliminate this cumulative error? If "u" > > is changed only in that one place, and you know the number of passes > > that are supposed to occur (i.e. in exact arithmetic), then you can > > write the loop as, for example,
> > do iter = 0, 199 > > u = (iter) * 0.1D0
> Definitely not. This multiplies the error in 0.1 by iter, > and as iter approaches 199, you have multiplied this > error by 199.
Of course, but that is the best that can be done with floating point arithmetic.
> A sure way to minimise the error is to divide by 10. > Thus, > u = iter/10.0d0 > has a small error, while iter*0.1d0 has a large error.
This is just drivel.
When most compilers see the constant in the denominator, they will compute the reciprocal and use that in the expression. That is, instead of doing the slow division each pass, they would do the faster multiply. There can be differences in the two results, but it should always be just in the last binary digit. Even when iter is a multiple of 10, you cannot expect the result of (iter/10.0d0) to be accurate to all bits. It might be correct, but you cannot depend on it (and the fortran standard does not require it).
> > ... > > enddo
> > (I suggest also that you parameterize your double precision constant > > with a real type, 0.1_dp, but this is a separate issue.)
> > With this loop, there is still a little error associated with each > > value of u, but it is only a local error (i.e. probably no more than > > the last bit only),
> The error is magnified as u increases in value.
Of course, that is the way floating point arithmetic works. The error in the last bit is always the same relative error, but its value changes with the magnitude of the value. That error, relative to exact arithmetic, cannot be eliminated.
It is useful to experiment with something like the following program.
program divide implicit none integer, parameter :: wp = selected_real_kind(12) integer :: i, i1, i2 real(wp) :: r1, r2 write(*,*) 'input range:' read(5,*) i1, i2 do i = i1, i2 r1 = i * 0.1_wp r2 = i / 10._wp if ( r1 .ne. r2 ) write(*,'(i0,2(f30.20,z20))') i, r1, r1, r2, r2 enddo end program divide
You will get different results with different compilers, or even with different options on a given compiler, but it will be clear that what I said above is true. The differences in the two values, when they occur, will probably be only in the last binary digit (as seen in the z20 formats), even for large values of "i".
This is in contrast to the original code in which the value was the result of repeated accumulations. I that case, the error really does grow beyond the last bit.
Ron Shepard <ron-shep...@nospam.comcast.net> wrote:
(someone wrote)
>> Definitely not. This multiplies the error in 0.1 by iter, >> and as iter approaches 199, you have multiplied this >> error by 199. > Of course, but that is the best that can be done with > floating point arithmetic.
