! fpio-test.f90
!
! K. Myneni, 16 Sep 2026
!
! Test decimal string to IEEE double precision floating
! point (binary64).
! 
! This is a port of the Forth floating point input tests.
!
! Build under FTN95:
!   ftn95 fpio-test.f90
!   slink fpio-test.obj -file:fpio-test.exe
!
! Build under gfortran:
!   gfortran -o fpio-test fpio-test.f90
!
program fpiotest
    implicit none
    integer idp(2), ref(2), test_num, e_s1, e_s2, e_tot
    integer, parameter :: dpkind = selected_real_kind (15, 307)
    real (kind=dpkind) dp
    equivalence(dp,idp)
    character(2) section

    e_tot = 0

    section = ' I'
    test_num = 1
    e_s1 = 0
    print *, 'Section ', section, ' Testing Conversion of Exactly Representable Numbers'
    ! I.1
    dp = 0.0D0
    ref = [0,0]
    call cmp_idp_ref()
    ! I.2
    dp = 9.99999935045640392457461415399766451285519391957298315801212d-39
    ref = [int(z'80000000'),int(z'380b38fb')]
    call cmp_idp_ref()
    ! I.3
    dp = -1.00000001335143196001808973960578441619873046875d-10
    ref = [int(z'e0000000'),int(z'bddb7cdf')]
    call cmp_idp_ref()
    ! I.4
    dp = 9.99999974737875163555145263671875d-05
    ref = [int(z'e0000000'),int(z'3f1a36e2')]
    call cmp_idp_ref()
    ! I.5
    dp = 0.100000001490116119384765625d0
    ref = [int(z'a0000000'),int(z'3fb99999')]
    call cmp_idp_ref()
    ! I.6
    dp = 1.0d0
    ref = [int(z'00000000'),int(z'3ff00000')]
    call cmp_idp_ref()
    ! I.7
    dp = -1.0d0
    ref = [int(z'00000000'),int(z'bff00000')]
    call cmp_idp_ref()
    ! I.8
    dp = 3.926990926265716552734375d-1
    ref = [int(z'60000000'),int(z'3fd921fb')]
    call cmp_idp_ref()
    ! I.9
    dp = 5.235987901687622070312500d-1
    ref = [int(z'40000000'),int(z'3fe0c152')]
    call cmp_idp_ref()
    ! I.10
    dp = 7.853981852531433105468750d-1
    ref = [int(z'60000000'),int(z'3fe921fb')]
    call cmp_idp_ref()
    ! I.11
    dp = 1.047197580337524414062500d0
    ref = [int(z'40000000'),int(z'3ff0c152')]
    call cmp_idp_ref()
    ! I.12
    dp = 1.178097248077392578125000d0
    ref = [int(z'80000000'),int(z'3ff2d97c')]
    call cmp_idp_ref()
    ! I.13
    dp = 1.570796370506286621093750d0
    ref = [int(z'60000000'),int(z'3ff921fb')]
    call cmp_idp_ref()
    ! I.14
    dp = 1.963495373725891113281250d0
    ref = [int(z'20000000'),int(z'3fff6a7a')]
    call cmp_idp_ref()
    ! I.15
    dp = 2.094395160675048828125000d0
    ref = [int(z'40000000'),int(z'4000c152')]
    call cmp_idp_ref()
    ! I.16
    dp = 2.356194496154785156250000d0
    ref = [int(z'80000000'),int(z'4002d97c')]
    call cmp_idp_ref()
    ! I.17
    dp = 2.617993831634521484375000d0
    ref = [int(z'c0000000'),int(z'4004f1a6')]
    call cmp_idp_ref()
    ! I.18
    dp = 2.748893499374389648437500d0
    ref = [int(z'e0000000'),int(z'4005fdbb')]
    call cmp_idp_ref()
    ! I.19
    dp = 3.141592741012573242187500d0
    ref = [int(z'60000000'),int(z'400921fb')]
    call cmp_idp_ref()
    ! I.20
    dp = 10.0d0
    ref = [int(z'00000000'),int(z'40240000')]
    call cmp_idp_ref()
    ! I.21
    dp = 1.0d1
    call cmp_idp_ref()
    ! I.22
    dp = 0.10d2
    call cmp_idp_ref()
    ! I.23
    dp = 0.010d3
    call cmp_idp_ref()
    ! I.24
    dp = 0.0000010d7
    call cmp_idp_ref()
    ! I.25
    dp = 0.000000000000010d15
    call cmp_idp_ref()
    ! I.26
    dp = 0.0000000000000000000000000000000000010d37
    call cmp_idp_ref()
    ! I.27
    dp = 1.0d10
    ref = [int(z'20000000'),int(z'4202a05f')]
    call cmp_idp_ref()
    ! I.28
    dp = 9999999933815812510711506376257961984d0
    ref = [int(z'40000000'),int(z'479e17b8')]
    call cmp_idp_ref()
    e_s1 = e_tot

    section = 'II'
    test_num = 1  ! reset test number to 1
    e_s2 = 0
    print *, 'Section ',section,' TESTING Rounding of Numbers'
    ! II.1
    dp = 1.0d-10
    ref = [int(z'd9d7bdbb'),int(z'3ddb7cdf')]
    call cmp_idp_ref()
    ! II.2
    dp = 2.71828182845904523536d0
    ref = [int(z'8b145769'),int(z'4005bf0a')]
    call cmp_idp_ref()
    ! II.3
    dp = 3.14159265358979323846d0
    ref = [int(z'54442d18'),int(z'400921fb')]
    call cmp_idp_ref()
    ! II.4
    dp = 3.518437208883201171875d+013
    ref = [int(z'00000002'),int(z'42c00000')]
    call cmp_idp_ref()
    ! II.5
    dp = 1.00000005960464477550d0
    ref = [int(z'10000000'),int(z'3ff00000')]
    call cmp_idp_ref()
    ! II.6
    dp = 5.00000000000000166533453693773481063544750213623046875d-1
    ref = [int(z'00000002'),int(z'3fe00000')]
    call cmp_idp_ref()
    ! II.7
    dp = 62.5364939768271845828d0
    ref = [int(z'd5aa7ca4'),int(z'404f44ab')]
    call cmp_idp_ref()
    ! II.8
    dp = 8.10109172351d-10
    ref = [int(z'aef0fd0c'),int(z'3e0bd5cb')]
    call cmp_idp_ref()
    ! II.9
    dp = 9214843084008499d0
    ref = [int(z'ec57761a'),int(z'43405e6c')]
    call cmp_idp_ref()
    ! II.10
    dp = 1.50000000000000011102230246251565404236316680908203125d0
    ref = [0, int(z'3ff80000')]
    call cmp_idp_ref()
    ! II.11
    dp = 9007199254740991.4999999999999999999999999999999995d0
    ref = [int(z'ffffffff'),int(z'433fffff')]
    call cmp_idp_ref()
    ! II.12
    dp = 1.000000000000000111022302462515654042363166809082031250d+00
    ref = [0, int(z'3ff00000')]
    call cmp_idp_ref()
    ! II.13
    dp = 1.000000000000000111022302462515654042363166809082031251d+00
    ref = [1, int(z'3ff00000')]
    call cmp_idp_ref
    e_s2 = e_tot - e_s1

    ! Show test results
    print *
    write(*, '(a,1X,I3)') 'Section  I Error Count:', e_s1
    write(*, '(a,1X,I3)') 'Section II Error Count:', e_s2
    write(*, '(a,1X,I3)') 'Total      Error Count:', e_tot

    CONTAINS
      subroutine cmp_idp_ref()
        character(15), parameter :: fmt2='(a,2X,z8,2X,z8)'
        if (.NOT.(idp(1)==ref(1)).AND.(idp(2)==ref(2))) then
          print *
          write(*,'(a,1X,a2,a,I3)') 'ERROR in test',section,'.',test_num
          ! output high 32, low 32 hex values
          write(*,fmt2) 'Obtained ', idp(2), idp(1)
          write(*,fmt2) 'Should be', ref(2), ref(1)
          e_tot = e_tot + 1
        endif
        test_num = test_num + 1 
      end subroutine
end program