Fun with 40 year old code

3

That would not work on an old IBM mainframe, which uses base 16, not base 2. You can check in the file i1mach.f, and you will find a small number of mostly obsolete machines for which I1mach(14) = 14, rather than the value 53 for IEEE double precision. For the menagerie of machines of the 70s and 80s, they had to convert between bits, nybbles, bytes, digits, etc.

Next, about the rounding of log10(2d0). We may speculate that some fudging was needed to get the answer to come out equal to a known result for the common case where I1mach(14) = 53. My speculation: the programmer tried at first

ANIB = D1MACH(5)*K/0.30103000D0

but that gave ANIB = 52.99999923659105…, and that made NIB = 52 (truncation to next lower integer), which was not “right”. Using 0.30102000 gives ANIB = 53.00175991, which leads to NIB = 53, which is what was wanted. I don’t understand why the ANINT intrinsic was not used.

Those were the days!

Additional comments:
In this context, we may take note of a hidden thorn that awaits us in the Netlib file d1mach.f. Once upon a time, in order to avoid the messy equivalence statements in that file, I replaced the executable statements by

      SELECT CASE (I)
      CASE (1)
         D1MACH = TINY(1d0)
      CASE (2)
         D1MACH = HUGE(1d0)
      CASE (3)
         D1MACH = EPSILON(1d0)/2
      CASE (4)
         D1MACH = EPSILON(1d0)
      CASE (5)
         D1MACH = LOG10(2d0)
      END SELECT

and that gave mysterious errors in my programs.

The problem is that LOG10 is declared as an integer array, and many compilers allowed 2d0 as an index into that array (allowing a double precision index as an extension to the language), so LOG10(2) is equivalenced to the upper 32 bits of DMACH(5). One fix is to replace the variable LOG10 by ILOG10, so that it does not clash with the intrinsic function LOG10. Alternatively, one could remove all the now-superfluous local variables from the routine.

3 Likes