Home     |     .Net Programming    |     cSharp Home    |     Sql Server Home    |     Javascript / Client Side Development     |     Ajax Programming

Ruby on Rails Development     |     Perl Programming     |     C Programming Language     |     C++ Programming     |     IT Jobs

Python Programming Language     |     Laptop Suggestions?    |     TCL Scripting     |     Fortran Programming     |     Scheme Programming Language


 
 
Cervo Technologies
The Right Source to Outsource

MS Dynamics CRM 3.0

Fortran Programming Language

Epsilon, Precision or Tiny ???


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 ??

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*(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

On May 8, 11:31 am, Bob Stryk <rast@tcq.nospam.net> wrote:

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

On 8 mei, 20:41, monir <mon@mondenet.com> wrote:

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? Fortran-77? 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 -ffpe-trap=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 "i386-marcel-freebsd"...
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/

Add to del.icio.us | Digg this | Stumble it | Powered by Megasolutions Inc