Fortran - Numeric enigma

Author

Numeric enigma

2004-05-12, 9:08 pm

I have two *identical* programs, executed as a console and win app
respectively, which produce non identical results. This is probably a
question for the vendor (CVF) but you're welcome to pitch in with
whatever causes you may come up with for the discrepancy.

program ash
*------------------------------------------------
* inv hyperbolic sine: y = asinh(x)
*------------------------------------------------
data tol/1.e-3/, eps/1.e-6/, x/.001/
arsh(x) = log(x + sqrt(x**2+1))

y = x
do while (abs(y-yo) .gt. tol*abs(y)+eps)
yo = y
y = y - (sinh(y) - x)/cosh(y)
enddo
write(*,*) ' x =', x, ' asinh(x) =', y, ' arsh(x) =', arsh(x)
end

subroutine ash
*------------------------------------------------
* inv hyperbolic sine: y = asinh(x)
*------------------------------------------------
!dec$ attributes dllexport :: ash
include 'usr_sym.inc'
data tol/1.e-3/, eps/1.e-6/, x/.001/
arsh(x) = log(x + sqrt(x**2+1))

y = x
do while (abs(y-yo) .gt. tol*abs(y)+eps)
yo = y
y = y - (sinh(y) - x)/cosh(y)
enddo
write(7,*) ' x =', x,' asinh(x) =', y,' arsh(x) =', arsh(x)
end

x = 1.0000000E-03 asinh(x) = 9.9999993E-04 arsh(x) = 1.0000233E-03
x = 1.0000000E-03 asinh(x) = 9.9999993E-04 arsh(x) = 9.9999993E-04

Whatever the reason for the *difference* produced by the arsh() it's not
obvious, except for the possibility that CVF doesn't port itself. ;)

James Van Buskirk

2004-05-12, 9:08 pm

"bv" <bvoh@Xsdynamix.com> wrote in message
news:409C2F36.33D69B8C@Xsdynamix.com...

> do while (abs(y-yo) .gt. tol*abs(y)+eps)

> Whatever the reason for the *difference* produced by the arsh() it's not
> obvious, except for the possibility that CVF doesn't port itself. ;)

One obvious candidate is the nonconforming usage of the
uninitialized variable yo. What results do you get if
you rewrite the deprecated [:)] DO WHILE loop as an unrestricted
DO with an EXIT clause?

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

Dr Ivan D. Reid

2004-05-12, 9:09 pm

On Sat, 08 May 2004 01:49:18 GMT, James Van Buskirk <not_valid@comcast.net>
wrote in <y0Xmc.1557$iF6.219209@attbi_s02>:
> "bv" <bvoh@Xsdynamix.com> wrote in message
> news:409C2F36.33D69B8C@Xsdynamix.com...

[color=darkred]
[color=darkred]
> One obvious candidate is the nonconforming usage of the
> uninitialized variable yo. What results do you get if
> you rewrite the deprecated [:)] DO WHILE loop as an unrestricted
> DO with an EXIT clause?

On a practical level, I'd expect that to possibly cause an
exception on the first pass if yo were initialised to an invalid pattern,
or almost equally unlikely immediate termination if yo were accidentally
close enough to y.

I've played around a bit in g77, and small changes in y don't
produce the "large" difference seen in the two different statement
function results (using double precision makes a slight difference:
program ash
implicit none
double precision darsh
....
darsh(x) = dlog(x + dsqrt(x**2+1.0d0))
....
write(*,*) arsh(.001),darsh(.001d0)

0.000999976764 0.000999999881)

I'd suspect that the two different environments use
different math libraries -- that's a question for the compiler vendor, I'd
have thought.

--
Ivan Reid, Electronic & Computer Engineering, ___ CMS Collaboration,
Brunel University. Ivan.Reid@brunel.ac.uk Room 40-1-B12, CERN
KotPT -- "for stupidity above and beyond the call of duty".

glen herrmannsfeldt

2004-05-12, 9:09 pm

Dr Ivan D. Reid wrote:

(snip)

> I've played around a bit in g77, and small changes in y don't
> produce the "large" difference seen in the two different statement
> function results (using double precision makes a slight difference:
> program ash
> implicit none
> double precision darsh
> ...
> darsh(x) = dlog(x + dsqrt(x**2+1.0d0))
> ...
> write(*,*) arsh(.001),darsh(.001d0)
>
> 0.000999976764 0.000999999881)
>
> I'd suspect that the two different environments use
> different math libraries -- that's a question for the compiler vendor, I'd
> have thought.

I thought it was the same executable run in two different
environments, though that could be two different dynamic
link libraries.

It might be different rounding mode for the math processor.
When y=0.001, y**2=.000001, 1+y**2=1.000001,
sqrt(1+y**2)=1.0000005, to eight digits, or 1. in single precision.
log(1.001) is pretty much 0.001, but is the 1.001 rounded up or
down from the nominal value?

Note that this function exaggerates any rounding errors, and it gets
much worse as y gets even smaller.

I am not sure how this is usually evaluated. Most of the common
math functions have well understood algorithms, including doing
it different ways for different parts of the domain.

-- glen

robin

2004-05-12, 9:09 pm

Subject: Numeric enigma
From: bv <bvoh@Xsdynamix.com>,Eclipse Software
Date: Sat, 08 May 2004 00:58:52 GMT
..
| I have two *identical* programs, executed as a console and win app
| respectively, which produce non identical results.
..
The code contains at least two bugs.
..
| This is probably a
| question for the vendor (CVF) but you're welcome to pitch in with
| whatever causes you may come up with for the discrepancy.
|
| program ash
| *------------------------------------------------
| * inv hyperbolic sine: y = asinh(x)
| *------------------------------------------------
| data tol/1.e-3/, eps/1.e-6/, x/.001/
| arsh(x) = log(x + sqrt(x**2+1))
|
| y = x
| do while (abs(y-yo) .gt. tol*abs(y)+eps)
..
This variable 'yo' is uninitialized.
..
| yo = y
| y = y - (sinh(y) - x)/cosh(y)
| enddo
| write(*,*) ' x =', x, ' asinh(x) =', y, ' arsh(x) =', arsh(x)
| end
..
This code appears incomplete -- there is no call for subroutine 'ash'.
..
| subroutine ash
| *------------------------------------------------
| * inv hyperbolic sine: y = asinh(x)
| *------------------------------------------------
| !dec$ attributes dllexport :: ash
| include 'usr_sym.inc'
..
The code is missing 'usr_sym.inc'.
..
| data tol/1.e-3/, eps/1.e-6/, x/.001/
| arsh(x) = log(x + sqrt(x**2+1))
|
| y = x
| do while (abs(y-yo) .gt. tol*abs(y)+eps)
..
Maybe 'yo' is uninitialised here too?
..
| yo = y
| y = y - (sinh(y) - x)/cosh(y)
| enddo
| write(7,*) ' x =', x,' asinh(x) =', y,' arsh(x) =', arsh(x)
..
Are we meant to assume that unit 7 produces the second
output line below?
..
| end
|
| x = 1.0000000E-03 asinh(x) = 9.9999993E-04 arsh(x) = 1.0000233E-03
| x = 1.0000000E-03 asinh(x) = 9.9999993E-04 arsh(x) = 9.9999993E-04
..
| Whatever the reason for the *difference* produced by the arsh() it's not
| obvious, except for the possibility that CVF doesn't port itself. ;)
..
Wouldn't it depend on what the definitions in use_sym.inc?

2004-05-12, 9:09 pm

"Dr Ivan D. Reid" wrote:
>
> I'd suspect that the two different environments use
> different math libraries -- that's a question for the compiler vendor, I'd
> have thought.

I checked both (exe & dll),

Microsoft (R) COFF Binary File Dumper Version 6.00.8447
Copyright (C) Microsoft Corp 1992-1998. All rights reserved.
Dump of file ash.exe
File Type: EXECUTABLE IMAGE
Image has the following dependencies:
DFORRT.dll
MSVCRT.dll

and then verified that "system" dir has only one copy of each; this
seems to eliminate, what could have been, the most plausible
explanation.

2004-05-12, 9:09 pm

robin wrote:
>
> This code appears incomplete -- there is no call for subroutine 'ash'.

That call is buried somewhere in WinMain, and (hopefully) has no effect
on numbers within the app.

> Are we meant to assume that unit 7 produces the second
> output line below?
>
> | x = 1.0000000E-03 asinh(x) = 9.9999993E-04 arsh(x) = 1.0000233E-03
> | x = 1.0000000E-03 asinh(x) = 9.9999993E-04 arsh(x) = 9.9999993E-04

Yes, win app may not write to unit "*".

> Wouldn't it depend on what the definitions in use_sym.inc?

No, it contains symbol export defs only; it has no bearing on the value
of "x" or evaluation of arsh(x). Note, it's the console app that failed.

Sponsored Links