Fun with 40 year old code

1

Here’s a little puzzle. Consider the following code, which can be found in various forms in the two SLATEC quadrature routines dgaus8 and dqnc79:

K = I1MACH(14)
ANIB = D1MACH(5)*K/0.30102000D0
NBITS = ANIB

Note that:

  • I1MACH(10) = B, the base.
  • I1MACH(14) = T, the number of base-B digits.
  • D1MACH( 5) = LOG10(B)

For one thing, notice the 0.30102000D0 variable, is this supposed to be log10(2.0d0) (which is 0.30102999566398120)? Is this a typo? It looks like they ended it with 2000 when they meant to type 3000? In single precision, log10(2.0) is 0.30103001. Was it just rounded wrong and then somebody added a D0 to it at some point when they converted the routine to double precision? Or is there some deep reason to make this value very slightly less than log10(2.0)?

So, assuming that was a typo, then the expression for ANIB just reduces to:

ANIB = I1MACH(14)
NBITS = ANIB

So why didn’t they just do that? If this really was the intent, then in modern Fortran I believe we can just replace it with:

nbits = digits(1.0_wp)

(where wp is whatever real kind you want to use), which gives us 24 (real32), 53 (real64), and 113 (real128), the same as the original code gives. Is there any danger in making this replacement? Can anybody see any good reason why the original code was written like that?

References

1 Like
2

Lots of history here, which basically led to the model
intrinsics in Fortran. Assuming you have an old Fortran
version that has dozens of BOZ values for different machines
with architectures that do not even exist anymore, basically
on any likely machine you will encounter this C code should
give you reference values:

/* Standard C source for D1MACH */
#include <stdio.h>
#include <float.h>
#include <math.h>
double d1mach_(long *i)
{
        switch(*i){
          case 1: return DBL_MIN; 
          case 2: return DBL_MAX;
          case 3: return DBL_EPSILON/FLT_RADIX;
          case 4: return DBL_EPSILON;
          case 5: return log10((double)FLT_RADIX);
          }
        fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i);
        exit(1); return 0; /* some compilers demand return values */
}

So for an IEEE-compliant common platform

    case 3: /* = FLT_RADIX  ^ - DBL_MANT_DIG for IEEE:  = 2^-53 = 1.110223e-16 = .5*DBL_EPSILON */
    case 4: /* = FLT_RADIX  ^ (1- DBL_MANT_DIG) = for IEEE:  = 2^52 = 4503599627370496 = 1/DBL_EPSILON */
    case 5: return log10(2.0);/* = M_LOG10_2 in Rmath.h */

On a typical modern machine you will very likely see these values

    value = 1.112536929253601E-308;
    value = 4.494232837155789E+307;
    value = 1.110223024625157E-016;
    value = 2.220446049250313E-016;
    value = 0.301029995663981;

Reference:

Phyllis Fox, Andrew Hall, Norman Schryer,
Algorithm 528,
Framework for a Portable Library,
ACM Transactions on Mathematical Software,
Volume 4, Number 2, June 1978, page 176-188.

This can all be done now with intrinsics, so some fun things to try that
might make sense of it all if you play with it is

function d1mach ( i )
!*****************************************************************************80
!
!! D1MACH returns double precision real machine constants.
!
!  Discussion:
!
!    Assuming that the internal representation of a double precision real
!    number is in base B, with T the number of base-B digits in the mantissa,
!    and EMIN the smallest possible exponent and EMAX the largest possible
!    exponent, then
!
!      D1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
!      D1MACH(2) = B**EMAX*(1-B**(-T)), the largest magnitude.
!      D1MACH(3) = B**(-T), the smallest relative spacing.
!      D1MACH(4) = B**(1-T), the largest relative spacing.
!      D1MACH(5) = log10(B).
!
!  Modified:
!
!    24 April 2007
!
!  Author:
!
!    Phyllis Fox, Andrew Hall, Norman Schryer
!
!  Reference:
!
!    Phyllis Fox, Andrew Hall, Norman Schryer,
!    Algorithm 528:
!    Framework for a Portable Library,
!    ACM Transactions on Mathematical Software,
!    Volume 4, Number 2, June 1978, page 176-188.
!
!  Parameters:
!
!    Input, integer I, chooses the parameter to be returned.
!    1 <= I <= 5.
!
!    Output, real ( kind = 8 ) D1MACH, the value of the chosen parameter.
!
!====================================================================================================================================
! intrinsics
! exponent     -  Exponent function
! fraction     -  Fractional part of the model representation
! nearest      -  Nearest representable number
! rrspacing    -  Reciprocal of the relative spacing
! scale        -  Scale a real value
! set_exponent -  Set the exponent of the model
!
! spacing      -  Smallest distance between two numbers of a given type
! digits       -  Significant digits function
! epsilon      -  Epsilon function
! maxexponent  -  Maximum exponent of a real kind
! minexponent  -  Minimum exponent of a real kind
! precision    -  Decimal precision of a real kind
! radix        -  Base of a model number
! range        -  Decimal exponent range of a real kind
! tiny         -  Smallest positive number of a real kind
! huge         -  Largest number of a kind
!====================================================================================================================================
  implicit none

  real ( kind = 8 ) d1mach
  integer i

  select case(i)
  case(1); d1mach = 4.450147717014403D-308
  case(2); d1mach = 8.988465674311579D+307
  case(3); d1mach = 1.110223024625157D-016
  case(4); d1mach = 2.220446049250313D-016
  case(5); d1mach = 0.301029995663981D+000
  case default
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'D1MACH - Fatal error!'
    write ( *, '(a)' ) '  The input argument I is out of bounds.'
    write ( *, '(a)' ) '  Legal values satisfy 1 <= I <= 5.'
    write ( *, '(a,i12)' ) '  I = ', i
    d1mach = 0.0D+00
    stop
  end select

end function d1mach

program testit
use, intrinsic :: iso_fortran_env, only : real_kinds,real32,real64,real128
integer,parameter :: wp=kind(0.0d0)
!real(kind=real128) :: d1M(5)
real(kind=real64) :: d1M(5)
integer,parameter :: b=radix(0.0_wp)
integer,parameter :: t = digits (0.0_wp )
integer,       parameter :: emin = minexponent ( 0.0_wp )
integer,       parameter :: emax = maxexponent ( 0.0_wp )
real(kind=wp), parameter :: xmin = tiny ( 0.0_wp )
real(kind=wp), parameter :: xmax = huge ( 0.0_wp )
real(kind=wp), parameter :: eps = epsilon ( 0.0_wp )
real ( kind = 8 ) :: d1mach
external d1mach

     write(*,*)'WP=',WP ,'kind'
     write(*,*)'B=',B ,'radix'
     write(*,*)'T=',T ,'digits'
     write(*,*)'EMIN=',EMIN, 'minexponent'
     write(*,*)'EMAX=',EMAX, 'maxexponent'
     write(*,*)'XMIN=',XMIN, 'tiny'
     write(*,*)'XMAX=',XMAX, 'huge'
     write(*,*)'EPS =',EPS,  'epsilon'

     write(*,'(a)')repeat('=',80)
     write(*,*)'definition'
     D1M(1) = real(B,KIND=real64)**(EMIN-1)               ! the smallest positive magnitude.
     write(*,'(g0)')d1m(1)
     write(*,'(g0,T30,A)') tiny(0.0d0),         ' tiny        - smallest positive number of a real kind'

     D1M(2) = nearest(real(B,KIND=real64)**EMAX*(1-B**(-T)),-1.0_real64)       ! the largest magnitude.
     write(*,'(g0)')d1m(2)
     write(*,'(g0,T30,A)') huge(0.0d0),         ' huge        - largest number of a kind'

     D1M(3) = real(B,KIND=real64)**(-T)                   ! the smallest relative spacing.
     write(*,'(g0)')d1m(3)

     D1M(4) = real(B,KIND=real64)**(1-T)                  ! the largest relative spacing.
     write(*,'(g0)')d1m(4)
     write(*,'(g0,T30,A)') epsilon(0.0d0),      ' epsilon     - Epsilon function'

     D1M(5) = log10(real(B,kind=wp))
     write(*,'(g0)')d1m(5)

     write(*,'(a)')repeat('=',80)
     write(*,*)'d1mach'
     write(*,'(*(g0))')d1mach(1), ' D1MACH(1) = B**(EMIN-1), the smallest positive magnitude.'
     write(*,'(*(g0))')d1mach(2), ' D1MACH(2) = B**EMAX*(1-B**(-T)), the largest magnitude.'
     write(*,'(*(g0))')d1mach(3), ' D1MACH(3) = B**(-T), the smallest relative spacing.'
     write(*,'(*(g0))')d1mach(4), ' D1MACH(4) = B**(1-T), the largest relative spacing.'
     write(*,'(*(g0))')d1mach(5), ' D1MACH(5) = log10(B).'

     write(*,'(a)')repeat('=',80)
     write(*,*)'intrinsics'
     write(*,'(g0,T30,A)') spacing(0.0d0),      ' spacing     - Smallest distance between two numbers of a given type'
     write(*,'(g0,T30,A)') digits(0.0d0),       ' digits      - Significant digits function'
     write(*,'(g0,T30,A)') maxexponent(0.0d0),  ' maxexponent - Maximum exponent of a real kind'
     write(*,'(g0,T30,A)') minexponent(0.0d0),  ' minexponent - Minimum exponent of a real kind'
     write(*,'(g0,T30,A)') precision(0.0d0),    ' precision   - Decimal precision of a real kind'
     write(*,'(g0,T30,A)') radix(0.0d0),        ' radix       - Base of a model number'
     write(*,'(g0,T30,A)') range(0.0d0),        ' range       - Decimal exponent range of a real kind'
end program testit

urbanjs@venus:/tmp$

1 Like
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
4

Thank you! I think you are correct on all counts. I overlooked the fact that B was not always 2.

Very interesting topic… I’m going to expand this into a blog post…

5

So, if we compare two ways to do this:

! slatec way:
 anib_real32  = LOG10(real(B,real32))*K/0.30102000_real32;   nbits_real32  = anib_real32
 anib_real64  = LOG10(real(B,real64))*K/0.30102000_real64;   nbits_real64  = anib_real64
      
! naive way:      
 anib_real32_alt  = LOG10(real(B,real32))*K/log10(2.0_real32);   nbits_real32_alt  = anib_real32_alt
 anib_real64_alt  = LOG10(real(B,real64))*K/log10(2.0_real64);   nbits_real64_alt  = anib_real64_alt

! note that K = imach(11) for real32, and imach14 for real64

We see that, for the following architectures from the i1mach routine, the naive way will give the wrong answer:

  • BURROUGHS 5700 SYSTEM (real32: 38 instead of 39, and real64: 77 instead of 78)
  • BURROUGHS 6700/7700 SYSTEMS (real32: 38 instead of 39, and real64: 77 instead of 78)

And you are correct, if you use ANINT for the naive way:

 anib_real32_alt  = LOG10(real(B,real32))*K/log10(2.0_real32);   nbits_real32_alt  = anint(anib_real32_alt)
 anib_real64_alt  = LOG10(real(B,real64))*K/log10(2.0_real64);   nbits_real64_alt  = anint(anib_real64_alt)

Then they all match, and you don’t need to employ the “magic” number. I also wonder why they didn’t do that. Was ANINT not available on some compilers maybe?

6

NINT() and ANINT() were new f77 intrinsic functions. Many of these numerical libraries date back to before f77 compilers were available. So once the old pre-f77 codes were working, the “if it ain’t broke, don’t fix it” rule applied.

1 Like
7

The last comment in the prologue of d1mach.f at Netlib says:

C   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)

The Fortran 90 (or 8X, with X a hexadecimal digit, as some used to joke at the time) standard was released, I think, in 1991.

So, Fortran 77 and ANINT had been available for about a decade by the time of the last listed revision. However, as @RonShepard says, “It it ain’t broke…” was followed. They could have substituted an expression using IFIX(+0.5), considering that only positive arguments were needed in d1mach().

1 Like
8

Well, in the new Quadpack, I’ll be updating it to remove the magic 0.30102000 number and use the ANINT function that has been around for 45 years. I think it’s safe to do that now. :smile: It will even still work on the Burroughs, just in case those come back in style.

9

The second line after the one with “0.30102000” in dgaus8.f now becomes the main puzzle. Here are the three lines:

      ANIB = D1MACH(5)*K/0.30102000D0 ! will be replaced using ANINT
      NBITS = ANIB
      NLMX = MIN(60,(NBITS*5)/8)      ! Mystery here

What is the significance of the upper limit, 60, as well as the multiplier 5/8?

P.S. (1 day later). The statements

ANIB = D1MACH(5)*K/0.30102000D0
NIB = ANIB

can be replaced by the modern equivalent (at least on machines whose FP operations use base 2; please check and confirm!)

NIB = DIGITS(0d0)

I think that we could use these modernization occasions to replace the [d,r,i]1mach calls by their intrinsic function equivalents whenever possible.

1 Like
10

Yep, 60 is also the size of some of the arrays. Lot of magic numbers in there!

If you look at dqnc79 you will see something similar…but different:

NLMX = MIN(99,(NBITS*4)/5)

Which also has the comment that in 1993: “Increase array size 80->99, and KMX 2000->5000 for SUN -r8 wordlength.”

11

Here’s the writeup: Degenerate Conic | Quadrature Quandary

The Java version gave me a chuckle. That’s never going to work on the Burroughs 5700. :slight_smile:

2 Likes