Fun with 40 year old code

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