> > 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*(x1-a)/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*(x1-a)/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*(x1-a)/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? Fortran-77? 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*(x1-a)/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