COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

### Numeric enigma

• Follow

```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.  ;)

```
 0
Reply bvoh (245) 5/8/2004 12:58:52 AM

```"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

```
 0
Reply not_valid (1681) 5/8/2004 1:49:18 AM

```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...

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

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".
```
 0
Reply Ivan.Reid (496) 5/8/2004 10:04:12 AM

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

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

```
 0
Reply gah (12254) 5/8/2004 6:12:58 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
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.

```
 0
Reply bvoh (245) 5/11/2004 1:53:38 AM

4 Replies
30 Views

5/22/2013 11:25:38 AM