! Fortran 95 free source form program kinds95.f90 by J F Harper, 
! Mathematics, Victoria University of Wellington, New Zealand.
! Comments and questions to john DOT harper AT vuw DOT ac DOT nz , please.
! It prints various properties of the available real and integer kinds,
! including whether each real kind is in IEEE754, x86 80-bit, or neither.
! It uses no intrinsic module because the F95 standard did not have them.
!
! The program worked with these compilers in Linux systems: gfortran 14.2.0,
! g95 0.94, ifort 2024.2, ifx 2025.0.0 and AMD flang 5.0.0, but both g95
! and flang raised ?? queries. Lfortran 0.41.0 can't compile it.
! 
! If there are at most 6 real kinds with different precisions or 6 integer 
! kinds with different ranges, all are tested (but see the Warnings below).
! If there may be more the program says so. IEEE754 offers 5 real kinds, and
! https://en.wikipedia.org/wiki/Extended_precision describes a 6th: x86 80-bit.
!
! Overflowing arithmetic operations and IEEE NaN values are not tested. Some
! compilers would give compile-time errors, others would crash at run time.
!
! Some things in this program comply with Fortran 95 where Fortran 2003 could
! be more concise, to make kinds95.f90 and kinds03.f90 minimally different.

! Warning 1. The lowest integer kind is selected_int_kind(1) because a compiler
!    disallowed selected_int_kind(0). Hence a 4-bit integer kind will not be
!    detected even if it is available.
! Warning 2: The Fortran standard requires exponent(inf)==huge(1) but with
!    inf read in from 'INF' that is tested only for the first three real kinds.
!    That allows this program to compile with all the compilers listed above.
! Warning 3: -0.0 should be written with a - sign in Fortran 2008 or later
! Recent changes:
! 06/11/24 exponent(inf)==huge(1) checked; see Warning 2.
! 29/04/25 "tiniest" mentioned but not calculated: some compilers can't
! 26/03/26 queries note wording improved
! 03/06/26 In realrealprops order epsilon,tiny,huge -> huge,epsilon,tiny
! 08/05/26 bugfix in checkreadinfnan
program kinds95
  implicit none
  integer,  parameter:: maxrk = 6, & ! max. number of real kinds 
       maxik = 6, & ! max. number of integer kinds
       nilist = 10  ! number of tested integer properties of reals
! Set up real kinds (there must be at least two. IEEE suggests 5, x86 one more)
  integer  ,parameter:: rk1 = selected_real_kind(1)
  real(rk1),parameter::  neg01 = -0.0_rk1, h1 = huge(neg01)
  integer  ,parameter::srk2 = selected_real_kind(precision(h1)+1), &
       rk2 = (srk2+rk1+sign(1,srk2)*(srk2-rk1))/2
!          = merge(srk2,rk1,srk2>0) but f95 needs the above workaround
  real(rk2),parameter::  neg02 = -0.0_rk2, h2 = huge(neg02)
  integer  ,parameter::srk3 = selected_real_kind(precision(h2)+1), &
       rk3 = (srk3+rk2+sign(1,srk3)*(srk3-rk2))/2
  real(rk3),parameter::  neg03 = -0.0_rk3, h3 = huge(neg03) 
! rk3 = srk3 if that's a valid real kind, rk2 if not
  integer  ,parameter::srk4 = selected_real_kind(precision(h3)+1), &
       rk4 = (srk4+rk3+sign(1,srk4)*(srk4-rk3))/2
  real(rk4),parameter::  neg04 = -0.0_rk4, h4 = huge(neg04)
! rk4 = srk4 if that's a valid real kind, rk3 if not
  integer  ,parameter::srk5 = selected_real_kind(precision(h4)+1), &
       rk5 = (srk5+rk4+sign(1,srk5)*(srk5-rk4))/2
  real(rk5),parameter::  neg05 = -0.0_rk5, h5 = huge(neg05)
! rk5 = srk5 if that's a valid real kind, rk4 if not
  integer  ,parameter::srk6 = selected_real_kind(precision(h5)+1), &
       rk6 = (srk6+rk5+sign(1,srk6)*(srk6-rk5))/2
  real(rk6),parameter::  neg06 = -0.0_rk6, h6 = huge(neg06)
! rk6 = srk6 if that's a valid real kind, rk5 if not
  real(rk1)::  ninf1(3)=0.0_rk1 ! to ensure initialization
  real(rk2)::  ninf2(3)=0.0_rk2
  real(rk3)::  ninf3(3)=0.0_rk3
  real(rk4)::  ninf4(3)=0.0_rk4
  real(rk5)::  ninf5(3)=0.0_rk5
  real(rk6)::  ninf6(3)=0.0_rk6
! set up integer kinds (compilers must offer at least one)
  integer    ,parameter:: &
       ik1 = selected_int_kind(1) ! Can't have (0): Silverfrost bug
  integer(ik1),parameter::  i1 = 1_ik1
  integer     ,parameter:: sik2 = selected_int_kind(range(i1)+1), &
       ik2 = (sik2+ik1+sign(1,sik2)*(sik2-ik1))/2
!          = merge(sik2,ik1,sik2>0) not an F95 constant, OK in F2003
  integer(ik2),parameter::  i2 = 1_ik2
! ik2 = sik2 if that's a valid integer kind, ik1 if not
  integer     ,parameter:: sik3 = selected_int_kind(range(i2)+1), &
       ik3 = (sik3+ik2+sign(1,sik3)*(sik3-ik2))/2
  integer(ik3),parameter::  i3  = 1_ik3
! ik3 = sik3 if that's a valid integer kind, ik2 if not
  integer     ,parameter:: sik4 = selected_int_kind(range(i3)+1), &
       ik4 = (sik4+ik3+sign(1,sik4)*(sik4-ik3))/2
  integer(ik4),parameter::  i4 = 1_ik4
! ik4 = sik4 if that's a valid integer kind, ik3 if not
  integer     ,parameter:: sik5 = selected_int_kind(range(i4)+1), &
       ik5 = (sik5+ik4+sign(1,sik5)*(sik5-ik4))/2
  integer(ik5),parameter::  i5 = 1_ik5
! ik5 = sik5 if that's a valid integer kind, ik4 if not
  integer     ,parameter:: sik6 = selected_int_kind(range(i5)+1), &
       ik6 = (sik6+ik5+sign(1,sik6)*(sik6-ik5))/2
  integer(ik6),parameter::  i6 = 1_ik6
! ik6 = sik6 if that's a valid integer kind, ik5 if not
  integer nrk, & ! number of real kinds, intent(out) in subroutine realkinds 
       realbits  ! number of bits in default real, intent(out) in realkinds 
  logical queries ! true if ?? written, evaluated in various subroutines
  character(132) line
  queries = .false.
  call realkinds(realbits,nrk)
  call intkinds
  call bits95check(realbits)
  if(queries) write(*,"( /,(1X,A))") 'NOTE: ?? indicates either'// &
         ' a real type neither IEEE754 nor x86 80-bit,', &
         ' or an option allowing nonstandard Fortran, or a compiler bug.'
contains
  
  subroutine realkinds(   realbits,nrk)
    integer,intent(out):: realbits,nrk
    integer:: i,stdpos,stdrk,ilist(maxrk,nilist),iolen=0,iolen1e0=0
    integer,parameter:: rkarray(maxik) = (/ rk1,rk2,rk3,rk4,rk5,rk6 /)
    logical,dimension(maxrk):: stdOK
    character:: neg0(maxrk)*4,crlist(3)*11,cilist(nilist)*22,fmt*32, &
         stdnames(maxrk)*10, ck(maxrk)*10,ck2(maxrk)*10
    real(rk6) rlist(maxrk,3)
    integer,parameter:: & ! half,single,double,extended,quad,octuple  
       stdbits(maxrk)   = (/  16,    32,    64,      80, 128,   256 /), &
       stddigits(maxrk) = (/  11,    24,    53,      64, 113,   237 /), &
       stdmaxexp(maxrk) = (/   4,     7,    10,      14,  14,    18 /)
    data stdnames/'half','single','double','extended','quad','octuple'/
    data cilist/'kind','minexponent','maxexponent  . . . . .','range', &
         'digits','precision  . . . . . .','radix','bits (from above)',&
         'bits (from IEEE+x86) .','bits (from iolength)'/
    data crlist/'huge','epsilon','tiny'/

    nrk = maxrk-count(rkarray(1:maxrk-1)==rkarray(2:maxrk))
    write(*,"( /A/ )") ' Properties of real kinds:'
    rlist(1,:) = real((/ h1,epsilon(h1),tiny(h1) /),rk6)
    rlist(2,:) = real((/ h2,epsilon(h2),tiny(h2) /),rk6)
    rlist(3,:) = real((/ h3,epsilon(h3),tiny(h3) /),rk6)
    rlist(4,:) = real((/ h4,epsilon(h4),tiny(h4) /),rk6)
    rlist(5,:) = real((/ h5,epsilon(h5),tiny(h5) /),rk6)
    rlist(6,:) = real((/ h6,epsilon(h6),tiny(h6) /),rk6)
    ilist(1,1:7) = (/ rk1,minexponent(h1),maxexponent(h1),range(h1), &
         digits(h1),precision(h1),radix(h1) /)
    ilist(2,1:7) = (/ rk2,minexponent(h2),maxexponent(h2),range(h2), &
         digits(h2),precision(h2),radix(h2) /)
    ilist(3,1:7) = (/ rk3,minexponent(h3),maxexponent(h3),range(h3), &
         digits(h3),precision(h3),radix(h3) /)
    ilist(4,1:7) = (/ rk4,minexponent(h4),maxexponent(h4),range(h4), &
         digits(h4),precision(h4),radix(h4) /)
    ilist(5,1:7) = (/ rk5,minexponent(h5),maxexponent(h5),range(h5), &
         digits(h5),precision(h5),radix(h5) /)
    ilist(6,1:7) = (/ rk6,minexponent(h6),maxexponent(h6),range(h6), &
         digits(h6),precision(h6),radix(h6) /)
    realbits = digits(1.0)+ilogbase(maxexponent(1.0),2) + 1
    inquire(iolength=iolen1e0) 1e0 ! e but not . allowed in name; 1e0=1.0
    if(iolen1e0<1.or.iolen1e0>4*realbits) iolen1e0 = huge(1)
    iolen = -1 ! in case inquire(iolength=iolen) fails
    ck2 = 'non-std.'
    do i = 1,maxrk
       if (i==1) inquire(iolength=iolen) h1
       if (i==1) ilist(i,8) = digits(h1)+ilogbase(maxexponent(h1),2) + 1
       if (i==2) inquire(iolength=iolen) h2
       if (i==2) ilist(i,8) = digits(h2)+ilogbase(maxexponent(h2),2) + 1
       if (i==3) inquire(iolength=iolen) h3
       if (i==3) ilist(i,8) = digits(h3)+ilogbase(maxexponent(h3),2) + 1
       if (i==4) inquire(iolength=iolen) h4
       if (i==4) ilist(i,8) = digits(h4)+ilogbase(maxexponent(h4),2) + 1
       if (i==5) inquire(iolength=iolen) h5
       if (i==5) ilist(i,8) = digits(h5)+ilogbase(maxexponent(h5),2) + 1
       if (i==6) inquire(iolength=iolen) h6
       if (i==6) ilist(i,8) = digits(h5)+ilogbase(maxexponent(h6),2) + 1
       ilist(i,10) = iolen*realbits/iolen1e0 ! probably 0 if iolen1e0 bad
       stdrk = minloc(abs(ilist(i,5)-stddigits),1)
       ilist(i, 9) = stdbits(stdrk)
       ck(i) = stdnames(stdrk)
       if(ilist(i,5)==64) then
          ck(i) = 'extended'
       else if (ilist(i,8)==bit_size(1)) then
          ck(i) = 'single'
       else if (ilist(i,8)==2*bit_size(1)) then
          ck(i) = 'double'
       else if (ilist(i,8)==4*bit_size(1)) then
          ck(i) = 'quad'
       else if (ilist(i,8)<bit_size(1))then
          ck(i) = merge('half ','lower',stddigits(stdrk)==ilist(i,5))
       else if (ilist(i,8)>4*bit_size(1)) then
          ck(i) = 'higher'
       else
          ck(i) = ' ??'
          queries = .true.
       end if
       if(ilist(i,9)==ilist(i,10))then
          ck2(i) = 'IEEE 754'
       else if(abs(ilist(i,8)-80)<=1) then
          ck2(i) = 'x86 80-bit'
       end if
    end do ! i = 1,maxrk
    do i = 1,nrk
       stdpos = minloc(abs(stddigits-ilist(i,5)),1)
       stdOK(i) = stddigits(stdpos)==ilist(i,5) .and. &
            stdmaxexp(stdpos)==ilogbase(ilist(i,3),radix(1.0))
    end do
    write(line,"(99(A,:,1X))")' precision called',(adjustr(ck(i)),i=1,nrk)
    call writeandquery(.not. all(stdOK(1:nrk)),line,nrk)
    if(any(len_trim(ck2)>0)) write(*,"(99(A,:,1X))") &
         ' standard obeyed ',(adjustr(ck2(i)),i=1,nrk)
! Give integer properties of real kinds (kind,minexponent,maxexponent,
!     range, digits, precision, radix, bits needed (3 ways)
    call realintprops(ilist,cilist)
! Give real properties of real kinds (epsilon,tiny,huge) ~ real(radix)**radpwr
    call realrealprops(ilist,rlist,crlist)
! Test whether -0.0 is printed with its sign
    write(neg0,'(F4.1,1X)') neg01,neg02,neg03,neg04,neg05,neg06
    write(*,"(1X,A,5X,99(A4,:,7X))")'-0.0 is written as',neg0(1:nrk)
    fmt = "(1X,A,T26,99(A3,:,8X))"
    call checkreadinfnan(fmt,ck)
    if(nrk>3) write(*,"(1X,A)") 'This program omits Inf tests that would' &
         //' require features missing from','some Fortran 95 compilers.'
    write(*,"( /,1X,A)") 'There '//merge('are no','may be',nrk<maxrk)// &
         ' real kinds with higher precision.'
  end subroutine realkinds

  elemental integer function ilogbase(number,base) ! e.g. ilogbase(8,2)==3
    integer,intent(in)::              number,base
    ilogbase = nint(log(real(number,kind(1d0)))/log(real(base,kind(1d0))))
  end function ilogbase

  subroutine writeandquery(bad, string,n,stringout) 
    logical,intent(in)::   bad
    character(*),intent(inout)::string
    integer,intent(in)::               n
    integer,intent(in),optional::        stringout
    integer blanksbeforequery
    blanksbeforequery = 11*(nrk-n)+1
! If(bad) queries = .true. and string has '??' added after some blanks. Then    
! it is written, internally if stringout is present, but to * (stdout) if not.
    string = trim(string)//repeat(' ',blanksbeforequery)//merge('??','  ',bad)
    if(present(stringout)) then
       write(string,"(A)")trim(string)
    else
       write(  *   ,"(A)")trim(string)
    end if
    if(bad) queries = .true.
  end subroutine writeandquery
  
  subroutine realintprops(ilist,cilist)
    integer,intent(in)::  ilist(:,:)
    character(*),intent(in)::  cilist(:)
    character(24) fmt
    integer i,trouble
    write(fmt,"(A,I0,A)") "(1X,A,I5,",nrk-1,"(3X,I8,:),A)"
    do i = 1,size(cilist)
       trouble = 0
       if(index(cilist(i),'radix')>0)  &
            trouble = trouble+merge(1,0,any(ilist(1:nrk,i)/=ilist(1,i)))
       if(index(cilist(i),'iolength')>0) &
            trouble = trouble+merge(1,0,any(ilist(1:nrk,i)<1))
       write(line,fmt) cilist(i),ilist(1:nrk,i)
       call writeandquery(trouble>0,line,nrk)
    end do
  end subroutine realintprops

  subroutine realrealprops(ilist,rlist,crlist)
    integer,intent(in)::   ilist(:,:)
    real(rk6),intent(in)::       rlist(:,:)
    character(*),intent(in)::          crlist(:)    
    integer:: i,j,bpower(nrk,3), b(nrk)
    logical rlistbad, differentb
    character cvalue(nrk,3)*11, cb*(10)
    b = ilist(1:nrk,7) ! b = radix = base of floating-point numbers
    differentb = any(b/=ilist(1,7))
    write(cb,"(F0.1)") real(b(1))
    do i = 1,nrk
       do j = 1,3
          write(cvalue(i,j),"(99(ES11.2E4,:))") rlist(i,j)
       end do
       bpower(i,1) = ilist(i,3)   ! huge = (1-epsilon/b)*b**maxexponent
       bpower(i,2) = 1-ilist(i,5) ! epsilon = b**(1-digits)
       bpower(i,3) = ilist(i,2)-1 ! tiny = b**(minexponent -1)
    end do
    do j = 1,3
       write(line,"(1X,A,5X,99(A,:))") crlist(j),cvalue(:,j)
       rlistbad = any(index(cvalue(:,j),'E')==0) .or. &
            any(index(cvalue(:,j),'0.00E')>0)
       call writeandquery(rlistbad,line,nrk)
       if(differentb) cycle ! b** line written only if all b are equal
       write(*,"(7X,A,T23,I6,99(I11,:))") &
            merge('near','   =',j==3)//' '//trim(cb)//' **',bpower(:,j)      
    end do
    write(*,"(A)") '(Least subnormal number is epsilon*tiny with '//&
    "compilers that allow it.)"
  end subroutine realrealprops

  subroutine checkreadinfnan(fmt,ck)
    character(*),intent(in)::fmt,ck(:)
    logical,dimension(maxrk):: nanOK,posOK,negOK,expOK(3)
! Read INF,-INF,NAN from cninf into ninf1, ninf2 etc. They can't be one array
    integer:: i, ios(maxrk)
    character(14):: cninf = ' INF -INF NAN '
    ios = 0
    write(*,"(1X,A)")'On reading real x from:'
    do i = 1,maxrk
       if(i==1)read(cninf,*,iostat=ios(i)) ninf1
       if(i==2)read(cninf,*,iostat=ios(i)) ninf2! kind(ninf2)/=kind(ninf1) etc
       if(i==3)read(cninf,*,iostat=ios(i)) ninf3
       if(i==4)read(cninf,*,iostat=ios(i)) ninf4
       if(i==5)read(cninf,*,iostat=ios(i)) ninf5
       if(i==6)read(cninf,*,iostat=ios(i)) ninf6
    end do
    if(all(ios(1:nrk)==0)) then
       do i = 1,3
          call readcheck(i,cninf)
       end do
       posOK = (/ ninf1(1)>h1,ninf2(1)>h2,ninf3(1)>h3,&
            ninf4(1)>h4,ninf5(1)>h5,ninf6(1)>h6 /)
       negOK = (/ ninf1(2)<-h1,ninf2(2)<-h2,ninf3(2)<-h3,&
            ninf4(2)<-h4,ninf5(2)<-h5,ninf6(2)<h6 /)
       nanOK = (/ ninf1(3)/=ninf1(3),ninf2(3)/=ninf2(3),ninf3(3)/=ninf3(3),&
            ninf4(3)/=ninf4(3),ninf5(3)/=ninf5(3),ninf6(3)/=ninf6(3) /)
       expOK = huge(1) == (/ exponent(ninf1(1)),exponent(ninf2(1)), &
            exponent(ninf3(1)) /) ! see Warning 2
       call askcheck(' Inf >  huge?',fmt,posOK,ios)
       call askcheck('-Inf < -huge?',fmt,negOK,ios)
       call askcheck(' NaN /= NaN?' ,fmt,nanOK,ios)
       call askcheck('exponent(Inf)==huge(1)?',fmt,expOK,ios)
    else
       do i = 1,nrk
          if(ios(i)/=0) write(*,*)'"'//trim(cninf)//'" cannot be read '//&
               'as three '//trim(ck(i))//'-precision "numbers"'
       end do
    end if
    if(.not.(all((/ nanOK,posOK,negOK /)))) write(*,*)'NaN or Inf errors '// &
         'are not compiler bugs if the -Ofast option was used.'
  end subroutine checkreadinfnan

  subroutine readcheck(  n,  cninf)
    integer,intent(in):: n
    character(*),intent(in)::cninf
    integer j
    character ninfkinds(maxrk)*11
    logical       infOK(nrk)
    write(ninfkinds,"(7X,F4.0)") & 
         ninf1(n),ninf2(n),ninf3(n),ninf4(n),ninf5(n),ninf6(n)
    infOK(1:nrk) = ninfkinds(1)==ninfkinds(1:nrk)
    write(line,"(1X,A,T20,99(A9,:,2X))") &
         "'"//adjustr(cninf(5*n-4:5*n-1))//"' x becomes",&
         (trim(adjustl(ninfkinds(j))),j=1,nrk) 
    call writeandquery(.not.all(infOK(1:nrk)),line,nrk)
  end subroutine readcheck
  
  subroutine askcheck(       question,format,isitOK,ios)
    character(*),intent(in)::question,format
    logical,dimension(:),intent(in)::        isitOK
    integer,dimension(:),intent(in)::               ios
    integer i,n
    n = min(nrk,size(isitOK))
    write(line, format) question,(tf(isitOK(i),ios(i)),i=1,n)
    call writeandquery(.not.all(isitOK(1:n)),line,n)
  end subroutine askcheck
  
  elemental character(3) function tf(OK,ios) ! returns 'Yes', ' No', or '   '
    logical, intent(in)::            OK
    integer, intent(in)::               ios
    tf =  merge(merge('Yes',' No',OK),'   ',ios==0)
  end function tf

  subroutine intkinds
    integer,parameter::ikarray(maxik) = (/ ik1,ik2,ik3,ik4,ik5,ik6 /)
    integer:: i,iolen=-1,k,nik ! nik = min(6,number of integer kinds)
    character:: approx*7,ck(maxik)*22 = ' '
    integer(ik6) ilist(maxik,6)
    write(*,"( /A/ )") ' Properties of integer kinds:'
    ilist(1,1:3) = (/ digits(i1),radix(i1),range(i1) /) ! RHS default integer
    ilist(2,1:3) = (/ digits(i2),radix(i2),range(i2) /) ! 
    ilist(3,1:3) = (/ digits(i3),radix(i3),range(i3) /)
    ilist(4,1:3) = (/ digits(i4),radix(i4),range(i4) /)
    ilist(5,1:3) = (/ digits(i5),radix(i5),range(i5) /)
    ilist(6,1:3) = (/ digits(i6),radix(i6),range(i6) /)
    do k = 1,maxik
       if (k==1) inquire(iolength=iolen) i1
       if (k==2) inquire(iolength=iolen) i2
       if (k==3) inquire(iolength=iolen) i3
       if (k==4) inquire(iolength=iolen) i4
       if (k==5) inquire(iolength=iolen) i5
       if (k==6) inquire(iolength=iolen) i6 
       ilist(k,4) = iolen
       ck(k) = trim(merge('(default integer)','                 ', &
            ikarray(k)==kind(1)))
       call writeandquery(bit_size(k)/=32.or.iolen<1,ck(k),nrk,666)
    end do
    ilist(1,5:6) = (/ bit_size(i1),huge(i1) /) ! RHS kind ik1, LHS kind(1)
    ilist(2,5:6) = (/ bit_size(i2),huge(i2) /) ! RHS kind ik2
    ilist(3,5:6) = (/ bit_size(i3),huge(i3) /) ! RHS kind ik3
    ilist(4,5:6) = (/ bit_size(i4),huge(i4) /) ! RHS kind ik4
    ilist(5,5:6) = (/ bit_size(i5),huge(i5) /) ! RHS kind ik5
    ilist(6,5:6) = (/ bit_size(i6),huge(i6) /) ! RHS kind ik6
    nik = maxik - count(ikarray(1:maxik-1)==ikarray(2:maxik))
    write(*,"(A)")' kind digits radix range iolen bit_  huge '
    write(*,"(A)")'                               size '
    do k = 1,nik
       write(*, "(I4,4I6,I5,4X,I0,1X,A)") &
            ikarray(k),(ilist(k,i),i=1,6),trim(ck(k))
       approx = merge(' approx','       ',ilist(k,6)>999)
       write(line,"(T42,A,T54,A,ES10.3)" ) &
            trim(hugesize(k,ilist)),' =',real(ilist(k,6),kind(1d0))
       write(*, "(A)")trim(line)//approx
    end do
    write(*,"(A)")' There '//merge('are no','may be',nik<maxik)//&
         ' integer kinds with wider range.'
  end subroutine intkinds

  character(16) function hugesize(k,ilist) ! len=16 in case radix is large
    integer     , intent(in)::    k
    integer(ik6), intent(in)::      ilist(:,:)
! k = kind; writes huge as radix**digits - 1 = ilist(k,2)**ilist(k,1) - 1
    write(hugesize,"(2(A,I0),A)") '= ',ilist(k,2),'**',ilist(k,1),' - 1'
  end function hugesize

  subroutine bits95check(realbits) ! can't use numeric_storage_size in f95
    integer,intent(in):: realbits
! Do default integers and reals use the same number of bits? They should.
    integer bits(2)
    bits(:) = (/ bit_size(1),realbits /)
    if(bits(1)/=bits(2)) &
         write(*,"( /,3(A,I0),A)") ' Bit sizes '// &
         'of 1,1.0 are ',bits(1),',',bits(2),': not equal ??'
  end subroutine bits95check

end program kinds95