iraf-v216 · Code · Issues (50) · Pull requests (81)

iraf.net pull request #106

Replace d1mach.f and r1mach.f by C sources

olebole merged 1 commit to iraf-community/iraf

---

olebole commented on 2017-10-16

The Fortran sources for d1mach.f and r1mach.f don’t work on 64 bit due to a non-standard length of INTEGER in IRAF. Take the following test source code (file machtest.x):

task machtest = t_machtest  
procedure t_machtest ()  
real r1mach()  
double d1mach()  
int i  
int failures  
begin  
  
# Check that all values are positive numbers  
  
    failures = 0  
    do i=1, 5 {  
        if (! (r1mach(i) == r1mach(i))) {  
            call printf("R1MACH(%d) is not a number\n")  
	         call pargi(i)  
            failures = failures + 1  
        } else if (r1mach(i) <= 0) {  
            call printf("R1MACH(%d) is not positive\n")  
	         call pargi(i)  
            failures = failures + 1  
	}  
    }  
    call printf("\n")  
    do i=1, 5 {  
        if (! (d1mach(i) == d1mach(i))) {  
            call printf("D1MACH(%d) is not a number\n")  
	         call pargi(i)  
            failures = failures + 1  
        } else if (d1mach(i) <= 0) {  
            call printf("D1MACH(%d) is not positive\n")  
	         call pargi(i)  
            failures = failures + 1  
	}  
    }  
    if (failures == 0) {  
        call printf("Simple consistency check passed.\n")  
    } else {  
        call printf("Simple consistency check has %d failures\n")  
	     call pargi(failures)  
    }  
  
# Print a summary of all values  
  
    call printf("\n%9c | %15s | DESCRIPTION\n")  
         call pargstr("R1MACH")  
         call pargstr("D1MACH")  
    call printf("%9e | %15e | THE SMALLEST POSITIVE MAGNITUDE\n")  
         call pargr(r1mach(1))  
         call pargd(d1mach(1))  
    call printf("%9e | %15e | THE LARGEST MAGNITUDE\n")  
         call pargr(r1mach(2))  
         call pargd(d1mach(2))  
    call printf("%9e | %15e | THE SMALLEST RELATIVE SPACING\n")  
         call pargr(r1mach(3))  
         call pargd(d1mach(3))  
    call printf("%9e | %15e | THE LARGEST RELATIVE SPACING\n")  
         call pargr(r1mach(4))  
         call pargd(d1mach(4))  
    call printf("%9e | %15e | LOG10(B)\n")  
         call pargr(r1mach(5))  
         call pargd(d1mach(5))  
end  

Checking against NaN is done here by simply comparing the value with itself (which is true for all numbers, and false otherwise).

Compiling and running with the unchanged code results in:

cl> xc machtest.x  
cl> task $machtest = machtest.e  
cl> machtest  
R1MACH(2) is not positive  
R1MACH(4) is not positive  
  
D1MACH(1) is not positive  
D1MACH(3) is not a number  
D1MACH(5) is not positive  
Simple consistency check has 5 failures  
  
   R1MACH |          D1MACH | DESCRIPTION  
1.175E-38 | 0.00000000000E0 | THE SMALLEST POSITIVE MAGNITUDE  
0.00000E0 | 0.00000000052E0 | THE LARGEST MAGNITUDE  
3.4028E38 | 0.0000000000E-1 | THE SMALLEST RELATIVE SPACING  
0.00000E0 | 0.00000106048E0 | THE LARGEST RELATIVE SPACING  
5.9605E-8 | 0.00000000000E0 | LOG10(B)  

The reason is that the code of R1MACH and D1MACH uses implicite assumptions about the lengths of INTEGER, REAL and DOUBLE PRECISION in the EQUIVALENCE statements (specifically that the length of REAL and DOUBLE PRECISION is twice the length of INTEGER), which are obviously not true. There is also no reason to put this into a machine dependent FORTRAN source file, since a C function can return the required values without having machine dependent source code.

When running with this PR, one gets:

cl> machtest  
  
Simple consistency check passed.  
  
   R1MACH |          D1MACH | DESCRIPTION  
1.175E-38 | 2.22507385851E0 | THE SMALLEST POSITIVE MAGNITUDE  
3.4028E38 | 1.797693135E308 | THE LARGEST MAGNITUDE  
5.9605E-8 | 1.110223025E-16 | THE SMALLEST RELATIVE SPACING  
1.1921E-7 | 2.220446049E-16 | THE LARGEST RELATIVE SPACING  
3.0103E-1 | 3.0102999566E-1 | LOG10(B)  

The wrongly printed value for D1MACH(1) comes from a bug in the printf subroutine; when f.e. multiplying the value by 10 the result is as expected.

A simple test (just check for positive numbers) is added to #36.

---

Commits

• Replace d1mach.f and r1mach.f by C sources [92b6dc61]

---

Last updated on 2017-10-16