rskey's Calculator forensics in Fortran
06-30-2014, 12:44 PM (This post was last modified: 03-13-2016 02:11 PM by HP67.)
Post: #1
 HP67 Senior Member Posts: 592 Joined: Dec 2013
rskey's Calculator forensics in Fortran
This page http://www.rskey.org/~mwsebastian/miscprj/radians.htm made me wonder what the various versions of Fortran I have would produce for the calculator forensics test. Using the info there as a basic guide I prepared this small piece of Fortran code:

Code:
PROGRAM foren   IMPLICIT NONE   REAL(8) :: PI = 3.14159265358979323846264_8   REAL(8) :: DR, RD   DR = PI / 180.0_8   RD = 180.0_8 / PI   WRITE (*,*) RD * ASIN(RD * ACOS(RD * ATAN(TAN(DR * COS(DR * SIN(DR * 9.0_8)))))) END PROGRAM foren

For VAX/VMS I used this piece of F66-conforming code for the F77 compiler:

Code:
 C C RSKEY FORENSICS C       REAL*8 PI /3.14159265358979323846264D0/       REAL*8 DR, RD, RESULT C       DR = PI / 180.0D0       RD = 180.0D0 / PI C       RESULT = RD * DASIN(RD * DACOS(RD * DATAN(DTAN(DR *       X        DCOS(DR * DSIN(DR * 9.0D0)))))) C       WRITE (6,100) RESULT   100 FORMAT (' ',F17.15)       STOP       END

For IBM OS/360 FORTRAN G & H I used this piece of FORTRAN IV code which is very similar to the VAX code except I had to use the FORTRAN IV function subroutine names for two trig functions whereas on the VAX the F77 compiler was able to either find the right function or use the one from F66, however it is implemented.

Code:
 C C RSKEY FORENSICS C       REAL*8 PI /3.14159265358979323846264D0/       REAL*8 DR, RD, RESULT C       DR = PI / 180.0D0       RD = 180.0D0 / PI C       RESULT = RD * DARSIN(RD * DARCOS(RD * DATAN(DTAN(DR *      X         DCOS(DR * DSIN(DR * 9.0D0)))))) C       WRITE (6,100) RESULT   100 FORMAT (' ',F17.15)       STOP       END

Here are the results with the Fortran compilers I have: Sorry for the bad alignment. For whatever reason code tags don't help much and I can't set the font to a fixed with font so although it looks almost correct in preview/posting mode once it is posted it doesn't.

Code:
 g95              8.999999999967708 gfortran        8.9999999998325695 ifort             8.99999999983257 OpenUH        8.99999999983257 Open64         8.99999999983257 Solaris Studio 8.999999999832553 VMS/VAX       8.999999999984616 IBM G&H        8.999999999967723 MS compiler V3.2 8.99999999989670 [Thanks to Paul Berger]

This is the version info for the various compilers:

Code:
 g95 --version G95 (GCC 4.0.3 (g95 0.94!) Jan 17 2013) Copyright (C) 2002-2008 Free Software Foundation, Inc. gfortran --version GNU Fortran (GCC) 4.7.1 Copyright (C) 2012 Free Software Foundation, Inc. ifort --version ifort (IFORT) 13.1.1 20130313 Copyright (C) 1985-2013 Intel Corporation.  All rights reserved. openf95 --version OpenUH 3.0.29 (based on Open64 Compiler Suite: Version 5.0) Built on: 2013-06-14 15:03:42 -0500 Thread model: posix Configured with --prefix=/opt/openuh/3.0.29 --enable-coarrays --enable-dragon-support --disable-caf-runtime GNU gcc version 4.2.0 (Open64 5.0 driver) openf95 --version Open64 Compiler Suite: Version 5.0 Built on: 2011-11-09 11:16:36 +0800 Thread model: posix GNU gcc version 4.2.0 (Open64 5.0 driver) f95 -V f90: Sun Fortran 95 8.6 Linux_i386 2011/11/16 (DEC VAX/VMS) Compaq Fortran 77 V6.6-201 (IBM) FORTRAN IV G LEVEL  21 (IBM) LEVEL 21.8 ( JUN 74 ) OS/360  FORTRAN H

I don't know what I expected to find but I was surprised to see that Intel's Fortran produces the same results as the open64 suite. I have not benchmarked big programs but Intel's compiler is supposedly to beat anything else pretty consistently for runtimes on Intel hardware. If that is true then it means they really did some important work on the base open64 code.

I have other compilers on other hardware platforms. I'll add the info to this thread when I get a chance.

The only thing I have left is for SOLARIS SPARC. Will possibly update that later.

It ain't OVER 'till it's 2 PICK
06-30-2014, 02:31 PM
Post: #2
 Paul Berger (Canada) Senior Member Posts: 431 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
The ancient IBM (Microsoft) ver 2 Fortran compiler gives a result of 8.99999999989670 that is with the math co-processor library linked running on a 486.
06-30-2014, 02:35 PM
Post: #3
 HP67 Senior Member Posts: 592 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
That's pretty good! Now I will have to test it on IBM's FORTRAN IV and see what we get.

It ain't OVER 'till it's 2 PICK
07-01-2014, 07:34 AM (This post was last modified: 07-01-2014 11:27 AM by Martin Paech.)
Post: #4
 Martin Paech Junior Member Posts: 19 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
Thank you for sharing this test!

(06-30-2014 02:35 PM)HP67 Wrote:  Now I will have to test it on IBM's FORTRAN IV and see what we get.

Unfortunately, I have currently no access to a FORTRAN compiler prior 77 standard, but I'm very interested in your result.

BTW: The F90 compiler (v3.1.3) on my PA-RISC workstation (HP J6750 running HP-UX 11i v1 with 64-bit a kernel) gives

8.99999999983257

Martin
07-01-2014, 07:53 AM
Post: #5
 HP67 Senior Member Posts: 592 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
Looks like a lot of stuff is based on Open64.

It ain't OVER 'till it's 2 PICK
03-08-2016, 09:09 AM (This post was last modified: 03-08-2016 09:12 AM by HP67.)
Post: #6
 HP67 Senior Member Posts: 592 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
Updated with new source and results of three additional compilers. Added DEC F77 for VMS/VAX and IBM G & H FORTRAN for OS/360. G and H results are identical since both compilers use the same library. If anybody has a copy of H Extended it will be illuminating to see what the results are with the HX library and REAL*16. There is an F77+ compiler for IBM mainframes also but I haven't had access to it in years.

I guess I should consolidate the results from the two other people who contributed. Guys, please let me know if it is ok and I'll copy your results into the first post also. And if you still have it, please post your source code. Thanks.

It ain't OVER 'till it's 2 PICK
Post: #7
 Paul Berger (Canada) Senior Member Posts: 431 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
(03-08-2016 09:09 AM)HP67 Wrote:  I guess I should consolidate the results from the two other people who contributed. Guys, please let me know if it is ok and I'll copy your results into the first post also. And if you still have it, please post your source code. Thanks.

Here you go you may consolidate my result if you wish. I had to look around a bit for it and found it on a disk that is not longer in a system.

Edit: I just compiled with a later version of the MS compiler V3.2 using the software floating point and I get the same result.

Code:
      PROGRAM foren C       IMPLICIT NONE        REAL*8  DR, RD, PI        PI = 3.14159265358979323846264        DR = PI / 180.0        RD = 180.0 / PI        WRITE (*,*) RD * DASIN(RD * DACOS(RD * DATAN(DTAN(DR * DCOS(DR *      CDSIN(DR * 9.0))))))       END
03-09-2016, 10:08 AM (This post was last modified: 03-09-2016 10:35 AM by Dieter.)
Post: #8
 Dieter Senior Member Posts: 2,398 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
(06-30-2014 12:44 PM)HP67 Wrote:  Here are the results with the Fortran compilers I have: Sorry for the bad alignment. For whatever reason code tags don't help much and I can't set the font to a fixed with font so although it looks almost correct in preview/posting mode once it is posted it doesn't.

Correct alignment is no problem. You just have to consider that the message editor always uses a proportional font, while everything within the two "code" tags is formatted with a fixed font. So simply use another (fixed font) editor for such tables and copy/paste them from there:

Code:
g95             8.999999999967708 gfortran        8.9999999998325695 ifort           8.99999999983257 OpenUH          8.99999999983257 Open64          8.99999999983257 Solaris Studio  8.999999999832553 VMS/VAX         8.999999999984616 IBM G&H         8.999999999967723

This leads to the question of the "perfect" forensic result for each compiler, depending on the implemented precision, i.e. the final value if every intermediate result is correctly rounded to working precision. No, it's not 9.000000.... ;-) You should expect to lose at least 5 digits in the final result.

For n-digit BCD calculators the perfect results have been discussed here earlier, but I think it's bit more tricky with n-bit binary arithmetics, once done in software and the other time by calling hardware FPU routines.

BTW... wouldn't a really optimizing compiler replace x=asin(acos(atan(tan(cos(sin(9)))))) with a simple x=9 ? ;-)

Dieter
03-09-2016, 09:27 PM
Post: #9
 Paul Dale Senior Member Posts: 1,467 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
(03-09-2016 10:08 AM)Dieter Wrote:  This leads to the question of the "perfect" forensic result for each compiler, depending on the implemented precision, i.e. the final value if every intermediate result is correctly rounded to working precision. No, it's not 9.000000.... ;-) You should expect to lose at least 5 digits in the final result.

I'd expect most modern compilers would use IEEE 754 mathematics. While this standard doesn't guarantee correctly rounded trigonometric functions, a lot of libraries do, at least for 64 bit reals.

Quote:BTW... wouldn't a really optimizing compiler replace x=asin(acos(atan(tan(cos(sin(9)))))) with a simple x=9 ? ;-)

Absolutely not. Optimising compilers go to a lot of effort to preserve the floating point semantics of what the programmer wrote. In this case, it would be possible for the compiler to replace the expression by x = correctly rounded result near 9. There are usually compiler options to turn of this strict enforcement and then it might produce x=9.

Pauli
03-13-2016, 02:15 PM
Post: #10
 HP67 Senior Member Posts: 592 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
(03-08-2016 05:31 PM)Paul Berger (Canada) Wrote:
Code:
      PROGRAM foren C       IMPLICIT NONE        REAL*8  DR, RD, PI        PI = 3.14159265358979323846264        DR = PI / 180.0        RD = 180.0 / PI        WRITE (*,*) RD * DASIN(RD * DACOS(RD * DATAN(DTAN(DR * DCOS(DR *      CDSIN(DR * 9.0))))))       END

Hi Paul,

thanks I have updated the first post with your results. Do you see any difference if you specify double precision constants? The way the code is written it would normally cause standard reals to have to be promoted to REAL*8. For instance 180.0 in double precision in FORTRAN IV thru F90 would probably have to be specified as 180.0D0. I can't remember when KIND came into vogue but I think from F95 something like 180.0_8 would be the way to code it. I have seen accuracy fall off when constants are not declared with the correct precision but your result seems to contradict that. I wonder if there is some default in the compiler you're using that handles that automagically.

It ain't OVER 'till it's 2 PICK
03-13-2016, 02:24 PM
Post: #11
 HP67 Senior Member Posts: 592 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
(03-09-2016 10:08 AM)Dieter Wrote:  Correct alignment is no problem. You just have to consider that the message editor always uses a proportional font, while everything within the two "code" tags is formatted with a fixed font. So simply use another (fixed font) editor for such tables and copy/paste them from there:

Are you saying to align them outside of the website and then paste?

(03-09-2016 10:08 AM)Dieter Wrote:  BTW... wouldn't a really optimizing compiler replace x=asin(acos(atan(tan(cos(sin(9)))))) with a simple x=9 ? ;-)

That isn't as outlandish as it might appear. I looked over Solaris Studio's code generation one time and I was amazed at the level of program understanding. It certainly eliminates constant expressions to a very high degree. It should be possible to do exactly what you're suggesting.

It ain't OVER 'till it's 2 PICK
03-13-2016, 04:51 PM
Post: #12
 Paul Berger (Canada) Senior Member Posts: 431 Joined: Dec 2013
RE: rskey's Calculator forensics in Fortran
(03-13-2016 02:15 PM)HP67 Wrote:
(03-08-2016 05:31 PM)Paul Berger (Canada) Wrote:

Hi Paul,

thanks I have updated the first post with your results. Do you see any difference if you specify double precision constants? The way the code is written it would normally cause standard reals to have to be promoted to REAL*8. For instance 180.0 in double precision in FORTRAN IV thru F90 would probably have to be specified as 180.0D0. I can't remember when KIND came into vogue but I think from F95 something like 180.0_8 would be the way to code it. I have seen accuracy fall off when constants are not declared with the correct precision but your result seems to contradict that. I wonder if there is some default in the compiler you're using that handles that automagically.

The format for double precision constants in the compiler is 180.0D0, and it makes no difference. One thing that is unexpected the compiler has 3 math libraries 8087 hardware, 8087 emulation, and "REGMATH.LIB". The REGMATH library is supposed to be a faster library, with reduced precision but it still gives the same result.

Paul.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)