




Fortran Programming Language









Epsilon, Precision or Tiny ???
Hello; In evaluating the Fortran statement: x = dacos(1.0D0*(x1a)/x2) with: xh = 0.20 x1 = 0.50*(1.+xh) x2 = 0.50*(1. xh) a = 0.20 I get : x = NaN instead of 0.0!! Q: Is it possible to use the Epsilon, Precision or Tiny intrinsic functions to avoid this difficulty ?? How ?? Thank you kindly. Monir
"monir" <mon @mondenet.com> wrote in message news:1178635745.415154.54050@w5g2000hsg.googlegroups.com... > Hello; > In evaluating the Fortran statement: > x = dacos(1.0D0*(x1a)/x2) > with: xh = 0.20 > x1 = 0.50*(1.+xh) > x2 = 0.50*(1. xh) > a = 0.20 > I get : x = NaN instead of 0.0!! > Q: Is it possible to use the Epsilon, Precision or Tiny intrinsic > functions to avoid this difficulty ?? How ??
Using double precision consistently gets you closer: double precision xh,x1,x2,a,x xh = 0.20d0 x1 = 0.50d0*(1.d0+xh) x2 = 0.50d0*(1.d0 xh) a = 0.20d0 x = acos(1.0D0*(x1a)/x2) Regards, Mike Metcalf
Michael Metcalf wrote: > "monir" <mon @mondenet.com> wrote in message > news:1178635745.415154.54050@w5g2000hsg.googlegroups.com... >> Hello; >> In evaluating the Fortran statement: >> x = dacos(1.0D0*(x1a)/x2) >> with: xh = 0.20 >> x1 = 0.50*(1.+xh) >> x2 = 0.50*(1. xh) >> a = 0.20 >> I get : x = NaN instead of 0.0!! >> Q: Is it possible to use the Epsilon, Precision or Tiny intrinsic >> functions to avoid this difficulty ?? How ?? > Using double precision consistently gets you closer: > double precision xh,x1,x2,a,x > xh = 0.20d0 > x1 = 0.50d0*(1.d0+xh) > x2 = 0.50d0*(1.d0 xh) > a = 0.20d0 > x = acos(1.0D0*(x1a)/x2) > Regards, > Mike Metcalf
An easy way to get NaN is for the argument of acos to be greater than one or less than minus one. Your values evaluated exactly in decimal notation yields exactly 1.0D0 but a computer works in binary which can lead to roundoff error. Even one bit above 1.0D0 will give NaN for acos. You may want to consider restricting the argument with something like: double precision xc,xc2 xc = 1.0D0*(x1a)/x2) xc2 = max(1.0,min(1.0,xc)) x = acos(xc2) There is also a discussion not opened so far about why the use of dacos instead of acos. Does it have something to do with the Fortran you are using? What kind? Fortran77? Fortran 90? Fortran 95? Bob Stryk
On May 8, 11:31 am, Bob Stryk <rast@tcq.nospam.net> wrote:
> Michael Metcalf wrote: > > "monir" <mon @mondenet.com> wrote in message > > news:1178635745.415154.54050@w5g2000hsg.googlegroups.com... > >> Hello; > >> In evaluating the Fortran statement: > >> x = dacos(1.0D0*(x1a)/x2) > >> with: xh = 0.20 > >> x1 = 0.50*(1.+xh) > >> x2 = 0.50*(1. xh) > >> a = 0.20 > >> I get : x = NaN instead of 0.0!! > >> Q: Is it possible to use the Epsilon, Precision or Tiny intrinsic > >> functions to avoid this difficulty ?? How ?? > > Using double precision consistently gets you closer: > > double precision xh,x1,x2,a,x > > xh = 0.20d0 > > x1 = 0.50d0*(1.d0+xh) > > x2 = 0.50d0*(1.d0 xh) > > a = 0.20d0 > > x = acos(1.0D0*(x1a)/x2) > > Regards, > > Mike Metcalf > An easy way to get NaN is for the argument of acos to be greater than > one or less than minus one. Your values evaluated exactly in decimal > notation yields exactly 1.0D0 but a computer works in binary which can > lead to roundoff error. Even one bit above 1.0D0 will give NaN for > acos. You may want to consider restricting the argument with something > like: > double precision xc,xc2 > xc = 1.0D0*(x1a)/x2) > xc2 = max(1.0,min(1.0,xc)) > x = acos(xc2) > There is also a discussion not opened so far about why the use of dacos > instead of acos. Does it have something to do with the Fortran you are > using? What kind? Fortran77? Fortran 90? Fortran 95? > Bob Stryk Hide quoted text  >
Double Precision alone would not solve the problem of possible round off errors in the limiting values of the inverse cosine argument. The following works fine (F77 & g95); ensuring: 0.0 <= arccos(x) <= PI [ 1.0 <= x <= 1.0] Double Precision x, xh, x1, x2, a xh = 0.20d0 x1 = 0.50d0*(1.0d0 + xh) x2 = 0.50d0*(1.0d0  xh) a = 0.20d0 x = 1.0d0*(x1a)/x2 x = dmax1(1.0d0, dmin1(1.0d0, x)) x = dacos(x) The result is x = 0.0 And, if a = 1.0d0, the result is x = PI ... Perfect! Thank you very much, Mike Metcalf and Bob Stryk, for your help. Monir
On 8 mei, 20:41, monir <mon@mondenet.com> wrote:
> On May 8, 11:31 am, Bob Stryk <rast @tcq.nospam.net> wrote: > > Michael Metcalf wrote: > > > "monir" <mon@mondenet.com> wrote in message > > >news:1178635745.415154.54050@w5g2000hsg.googlegroups.com... > > >> Hello; > > >> In evaluating the Fortran statement: > > >> x = dacos(1.0D0*(x1a)/x2) > > >> with: xh = 0.20 > > >> x1 = 0.50*(1.+xh) > > >> x2 = 0.50*(1. xh) > > >> a = 0.20 > > >> I get : x = NaN instead of 0.0!! > > >> Q: Is it possible to use the Epsilon, Precision or Tiny intrinsic > > >> functions to avoid this difficulty ?? How ?? > > > Using double precision consistently gets you closer: > > > double precision xh,x1,x2,a,x > > > xh = 0.20d0 > > > x1 = 0.50d0*(1.d0+xh) > > > x2 = 0.50d0*(1.d0 xh) > > > a = 0.20d0 > > > x = acos(1.0D0*(x1a)/x2) > > > Regards, > > > Mike Metcalf > > An easy way to get NaN is for the argument of acos to be greater than > > one or less than minus one. Your values evaluated exactly in decimal > > notation yields exactly 1.0D0 but a computer works in binary which can > > lead to roundoff error. Even one bit above 1.0D0 will give NaN for > > acos. You may want to consider restricting the argument with something > > like: > > double precision xc,xc2 > > xc = 1.0D0*(x1a)/x2) > > xc2 = max(1.0,min(1.0,xc)) > > x = acos(xc2) > > There is also a discussion not opened so far about why the use of dacos > > instead of acos. Does it have something to do with the Fortran you are > > using? What kind? Fortran77? Fortran 90? Fortran 95? > > Bob Stryk Hide quoted text  > > > Double Precision alone would not solve the problem of possible round > off errors in the limiting values of the inverse cosine argument. > The following works fine (F77 & g95); ensuring: > 0.0 <= arccos(x) <= PI [ 1.0 <= x <= 1.0] > Double Precision x, xh, x1, x2, a > xh = 0.20d0 > x1 = 0.50d0*(1.0d0 + xh) > x2 = 0.50d0*(1.0d0  xh) > a = 0.20d0 > x = 1.0d0*(x1a)/x2 > x = dmax1(1.0d0, dmin1(1.0d0, x)) > x = dacos(x) > The result is x = 0.0 > And, if a = 1.0d0, the result is x = PI ... Perfect! > Thank you very much, Mike Metcalf and Bob Stryk, for your help. > Monir
I agree with Bob here: use the generic names instead of the specific names. If the type changes, you would have to change the names of the functions as well. And now you need to know the precise specific names. Leave that to the compiler. Use: max instead of dmax1 etc. (FORTRAN 77 compilers already knew both the generic and the specific names). Regards, Arjen
Bob Stryk <rast @tcq.nospam.net> wrote: > There is also a discussion not opened so far about why the use of dacos > instead of acos. Does it have something to do with the Fortran you are > using? What kind? Fortran77? Fortran 90? Fortran 95? Generic intrinsics came into the language with f77. One would have to be using f66 or earlier (which doesn't seem likely) for that to be an issue. It is my general impression that most uses of specific intrinsics come from being trained by people who learned f66 and pass their f66 habits on.  Richard Maine  Good judgement comes from experience; email: last name at domain . net  experience comes from bad judgement. domain: summertriangle   Mark Twain
Bob Stryk wrote: > An easy way to get NaN is for the argument of acos to be greater than > one or less than minus one. Your values evaluated exactly in decimal > notation yields exactly 1.0D0 but a computer works in binary which can > lead to roundoff error. Even one bit above 1.0D0 will give NaN for > acos.
In days past this would flag a "domain error" message (with the name of the offending function), and now you're telling us we are fed with this evangelical nonsense, outrageous!
In article <464112B7.6E861@cyber.com>, subspace <n@cyber.com> writes: > Bob Stryk wrote: >> An easy way to get NaN is for the argument of acos to be greater than >> one or less than minus one. Your values evaluated exactly in decimal >> notation yields exactly 1.0D0 but a computer works in binary which can >> lead to roundoff error. Even one bit above 1.0D0 will give NaN for >> acos. > In days past this would flag a "domain error" message (with the name of > the offending function), and now you're telling us we are fed with this > evangelical nonsense, outrageous!
In the present time, people are expected to read the documentation that comes with one's compiler. mobile:kargl[203] cat > j.f90 real x x = 1.1 y = acos(x) print *, x, y end mobile:kargl[204] gfc4x o z j.f90 mobile:kargl[205] ./z 1.100000 NaN mobile:kargl[208] gfc4x o z ffpetrap=invalid g j.f90 mobile:kargl[209] ./z Floating exception (core dumped) mobile:kargl[210] gdb ./z z.core GNU gdb 6.1.1 [FreeBSD] Copyright 2004 Free Software Foundation, Inc. GDB is free software, covered by the GNU General Public License, and you are welcome to change it and/or distribute copies of it under certain conditions. Type "show copying" to see the conditions. There is absolutely no warranty for GDB. Type "show warranty" for details. This GDB was configured as "i386marcelfreebsd"... Core was generated by `z'. Program terminated with signal 8, Arithmetic exception. #0 0x080481f2 in MAIN__ () at j.f90:3  Steve http://troutmask.apl.washington.edu/~kargl/





