!!----------------------------------------------------------------------------
!!
!!           file : demo.f90
!!
!!       function : This is a FORTRAN source file intended to test a subset
!!                  of the features of the Lahey/Fujitsu LF95 compiler and its
!!                  runtime support. Successful execution of this program
!!                  indicates that installation has been successful.
!!
!!  compatibility : This FORTRAN file can be compiled and run with all Lahey
!!                  LF95 language systems.  Note that the code in this
!!                  use of extensions to the FORTRAN 77 standard.
!!
!!  to create/run : Enter the following two commands in the order shown:
!!
!!                      [1]  lf95 demo.f90 <CR>
!!                      [2]  demo <CR>
!!
!!    portability : This code is portable to any system running Lahey LF95
!!                  or LF90
!!
!!    limitations :
!!
!!        warnings:
!!
!!        remarks : Many of the functions below are written in recursive form.
!!                  Those concerned with performance may wish to recast them
!!                  in the iterative mode.
!!
!!                  Three of the examples in this demonstration program were
!!                  inspired by Abelson and Sussman, The Structure and
!!                  Interpretation of Computer Programs, Cambridge, MA: The
!!                  MIT Press, 1985.
!!
!!
!!     DEMO.F90  Copyright(c) 1994 - 1998, Lahey Computer Systems, Inc.
!!     Copying for sale requires permission from Lahey Computers Systems.
!!     Otherwise, distribution of all or part of this file is permitted 
!!     if these four lines are included.
!!
!!----------------------------------------------------------------------------
        !! module constants_module
        !! contains parameters for declaring KIND= and the like
        module constants_module
            integer, parameter :: large_int_kind = 4
            integer, parameter :: large_real_kind = 8
        end module constants_module

        !! module utilities_module
        !! for doing little common things
        module utilities_module
        contains
            subroutine mypause()
            print *, ""
            call system("PAUSE")
            end subroutine mypause
        end module utilities_module

        !!
        !!   main program unit : LAHEY_DEMO
        !!
        !!      sets up and governs the selection between a defined set of
        !!      demonstration subprograms
        !!
        program lahey_demo
            integer ::selection
            logical :: toggle
            toggle = .false.
            call system("cls")
            print *, " "
            print *, "             Lahey/Fujitsu LF95 Compiler"
            print *, "             ---------------------------"
            print *, " "
            print *, "  installation test and demonstration program"
            print *, " "
            print *, "              Copyright(c) 1998"
            print *, "         Lahey Computer Systems, Inc."
  131       if (toggle) then
                call system("cls")
            end if
            toggle = .true.
            print *, " "
            print *, "    -----------------"
            print *, "    Test/Action List:"
            print *, "    -----------------"
            print *, "     1 - factorials"
            print *, "     2 - Fahrenheit to Celsius conversion"
            print *, "     3 - Carmichael numbers"
            print *, "     4 - Ramanujan's series"
            print *, "     5 - Stirling numbers of the 2nd kind"
            print *, "     6 - chi-square quantiles"
            print *, "     7 - Pythagorean triplets"
            print *, "     8 - date_and_time, and other system calls"
            print *, "     0 - <stop this program>"
            print *, " "
            print *, "    Please select an option by entering the"
            print *, "    associated number followed by <return>.    "
            read *, selection
            call system("cls")
            select case (selection)
                case (0)
                    stop
                case (1)
                    call factorial_demo
                case (2)
                    call conversion_demo
                case (3)
                    call carmichael_demo
                case (4)
                    call ramanujan_demo
                case (5)
                    call stirling_demo
                case (6)
                    call chi_square_demo
                case (7)
                    call pythagoras_demo
                case (8)
                    call runtime_funcs_demo
            end select
            go to 131
        end program lahey_demo



        !!
        !!  subroutine : FACTORIAL_DEMO
        !!
        !!      sets up the factorial demonstration and governs the
        !!      invocation of 'factorizal'
        !!
        subroutine factorial_demo
            use constants_module
            use utilities_module
            integer(kind=large_int_kind) :: i,f,factorial
            print *, "    demo #1: factorials"
            print *, "    ___________________"
            print *, " "
            print *, " "
            print *, "    This routine calculates and prints out the"
            print *, "    factorials of the integers from 1 through 12."
            print *, " "
            call mypause
            call system("cls")
            print *, " "
            i = 1
            do 
                if (i > 12) exit  
                f = factorial(i)
                print "('     n = ',I2.1,                               &
     &                      '      factorial(n) = ',I9.1)", i,f
                i = i + 1
            end do
            print *, " "
            call mypause
        end subroutine factorial_demo


        recursive function factorial(n) result(factresult)
            use constants_module
            integer(kind=large_int_kind) :: factresult
            integer(kind=large_int_kind) ::  n
            if (n .eq. 1) then
                factresult = 1
            else
                factresult = n * factorial(n - 1)
            end if
        end function factorial



        !!
        !!  subroutine : CONVERSION_DEMO
        !!
        !!      sets up the Fahrenheit to Celsius conversion demonstration
        !!
        subroutine conversion_demo
            use constants_module
            use utilities_module
            real(kind=large_real_kind) ::  h, l, i, celsius, fahr
            print *, "    demo #2: Fahrenheit to Celsius conversion"
            print *, "    _________________________________________"
            print *, " "
            print *, " "
            print *, "    This routine calculates and prints out the"
            print *, "    Celsius equivalent of every fifth integral"
            print *, "    Fahrenheit value from 5 to 100 degrees."
            print *, " "
            call mypause
            call system("cls")
            l = 5
            h = 100
            i = 5
            print *, '     deg F      deg C'
            print *, '    -------    -------'
            do fahr = l, h, i
                celsius = (5.0/9.0) * (fahr - 32.0)
                print "('     ',F7.2,'    ',F7.2)", fahr, celsius
            end do
            call mypause
        end subroutine conversion_demo



        !!
        !!  subroutine : CARMICHAEL_DEMO
        !!
        !!      sets up the Carmichael number demonstration and governs
        !!      the invocation of 'carmichaels_in_range'
        !!
        !!      Abelson & Sussman's "Structure and Interpretation of
        !!      of Computer Programs" provided the inspiration for
        !!      this example
        !!
        subroutine carmichael_demo
            use utilities_module
            print *, "    demo #3: Carmichael numbers"
            print *, "    ___________________________"
            print *, " "
            print *, " "
            print *, "    A Carmichael number is a non-prime that fools"
            print *, "    the Fermat test for primality.  These numbers"
            print *, "    are rare.  This program examines the integers"
            print *, "    from 2 to 3,000 and lists out the Carmichael"
            print *, "    numbers discovered."
            print *, " "
            call mypause
            call system("cls")
            print *, " "
            call carmichaels_in_range(2,3000)
            print *, " "
            call mypause
        end subroutine carmichael_demo


        subroutine carmichaels_in_range(ialpha,omega)
            use constants_module
            integer(kind=large_int_kind) :: alpha,omega,ialpha
            alpha=ialpha
            do 
                if (alpha > omega) exit
                call test_carmichael(alpha)
                alpha = alpha+1
            end do
        end subroutine carmichaels_in_range


        subroutine test_carmichael(n)
            use constants_module
            integer(kind=large_int_kind) :: n
            logical :: fermat_test,naive_test,f
            f = fermat_test(n)
            if (f .and. .not. naive_test(n)) then
                print*,n," is a Carmichael number!"
            end if
        end subroutine test_carmichael


        logical function fermat_test(int)
            use constants_module
            integer(kind=large_int_kind) :: int
            logical :: fast_prime
            fermat_test = fast_prime(int,9)
        end function fermat_test


        recursive function fast_prime(i,times) result(fast_result)
            use constants_module
            integer(kind=large_int_kind) :: i, times
            logical :: fermat,fast_result
            if (times .eq. 0) then
                fast_result = .true.
            else
                if (fermat(i) .eqv. .true.) then
                    fast_result = fast_prime(i,times-1)
                else
                    fast_result = .false.
                end if
            end if
        end function fast_prime


        logical function fermat(int)
            use constants_module
            integer(kind=large_int_kind) :: int, alpha,expmod
            alpha = 2 + ((int - 2) * rnd())
            fermat = (alpha .eq. expmod(alpha,int,int))
        end function fermat


        recursive function expmod(b,e,m) result(exp_result)
            use constants_module
            integer(kind=large_int_kind) :: b,e,m,tmp,exp_result
            if (0 .eq. e) then
                exp_result = 1
            else if (0 .eq. mod(e,2)) then
                tmp = expmod(b,e/2,m)
                exp_result = mod(tmp*tmp, m)
            else
                tmp = expmod(b,e-1,m)
                exp_result = mod(b*tmp, m)
            end if
        end function expmod


        logical function naive_test(n)
            use constants_module
            integer(kind=large_int_kind) :: n,smallest_divisor
            naive_test = (n .eq. smallest_divisor(n))
        end function naive_test


        function smallest_divisor(n)
            use constants_module
            integer(kind=large_int_kind) :: smallest_divisor
            integer(kind=large_int_kind) :: n,find_divisor
            smallest_divisor = find_divisor(n,2)
        end function smallest_divisor


        recursive function find_divisor(n, t) result(divisor)
            use constants_module
            integer(kind=large_int_kind) :: n,t,divisor
            if ((t * t) > n) then
                divisor = n
            else if (0 .eq. mod(n,t)) then
                divisor = t
            else
                divisor = find_divisor(n,t+1)
            end if
        end function find_divisor





        !!
        !!  subroutine : RAMANUJAN_DEMO
        !!
        !!      sets up the Ramanujan series demonstration and governs
        !!      the invocation of 'search_ramanujans'
        !!
        !!      Abelson & Sussman's "Structure and Interpretation of
        !!      of Computer Programs" provided the inspiration for
        !!      this example
        !!
        subroutine ramanujan_demo
            use utilities_module
            print *, "    demo #4: Ramanujan's series"
            print *, "    ___________________________"
            print *, " "
            print *, " "
            print *, "    G. H. Hardy, in his obituary for Srinavasa"
            print *, "    Ramanujan, told of the time that he visited"
            print *, "    Ramanujan when the latter was ill.  Hardy "
            print *, "    took taxi number '1729' to where Ramanujan"
            print *, "    was staying, and in passing told Ramanujan"
            print *, "    that it seemed a very uninteresting number."
            print *, "    Ramanujan instantly took issue with Hardy,"
            print *, "    pointing out that it is on the contrary"
            print *, "    very interesting,  because it is the smallest"
            print *, "    positive integer that can be expressed as the"
            print *, "    sum of 2 cubes in exactly 2 different ways."
            print *, "    The number '1729' has since been known as"
            print *, "    Ramanujan's number and the series it begins"
            print *, "    as Ramanujan's series."
            print *, " "
            print *, "    This program prints out the Ramanujan numbers"
            print *, "    discovered between 1 and 15,000."
            print *, " "
            call mypause
            call system("cls")
            print *, " "
            call search_ramanujans(1,15000)
            print *, " "
            call mypause
        end subroutine ramanujan_demo


        subroutine search_ramanujans(ialpha,omega)
            use constants_module
            integer(kind=large_int_kind) :: alpha,omega,ialpha
            logical :: is_ramanujan
            alpha=ialpha
            do 
                if (alpha > omega) exit
                if (0 .eq. mod(alpha,1000)) then
                    print*,"at:",alpha
                end if
                if (is_ramanujan(alpha)) then
                    print*,alpha," is a Ramanujan number!"
                end if
                alpha = alpha+1
            end do
        end subroutine search_ramanujans


        logical function is_ramanujan(n)
            use constants_module
            integer(kind=large_int_kind) :: n,lim
            logical :: search
            lim = 1 + int((n/2)**(1.0/3.0))
            is_ramanujan = search(n,1,lim,0)
        end function is_ramanujan


         recursive function search(n,a,o,cubes) result(searchres)
            use constants_module
            integer(kind=large_int_kind) :: n,a,o,cubes,ca,re,rc
            logical :: searchres
            if (a >= o) then
                if (2 .eq. cubes) then
                    searchres = .true.
                else
                    searchres = .false.
                end if
            else
                ca = a*a*a
                re = n - ca
                rc = int(re**(1.0/3.0))
                if (rc*rc*rc .eq. re) then
                    searchres = search(n,a+1,o,cubes+1)
                else
                    searchres = search(n,a+1,o,cubes)
                end if
            end if
        end function search



        !!
        !!  subroutine : STIRLING_DEMO
        !!
        !!      sets up the demonstration that calculates Stirling
        !!      numbers of the second kind and governs the invocation
        !!      of 'stirling_numbers'
        !!
        subroutine stirling_demo
            use utilities_module
            print *, "    demo #5: Stirling numbers of the 2nd kind"
            print *, "    _________________________________________"
            print *, " "
            print *, " "
            print *, "    Stirling numbers of the 2nd kind, S(n,k),"
            print *, "    are the numbers of possible partitions of"
            print *, "    n items into k groups.  They are useful in"
            print *, "    calculating the expected probability of"
            print *, "    some kinds of combinations of events."
            print *, " "
            print *, "    This program prints out Stirling numbers"
            print *, "    of the 2nd kind from S(12,0) to S(12,12)."
            print *, " "
            call mypause
            call system("cls")
            print *, " "
            call stirling_numbers(12,12,0,12)
            print *, " "
            call mypause
        end subroutine stirling_demo


        subroutine stirling_numbers(ia1,io1,ia2,io2)
            use constants_module
            integer(kind=large_int_kind) :: ia1,io1,ia2,io2
            integer(kind=large_int_kind) ::  a1, o1, a2, o2, ia
            integer(kind=large_int_kind) :: s,s2
            a1=ia1
            o1=io1
            a2=ia2
            o2=io2
            do 
                if (a1 > o1) exit
                ia = a2
                do 
                    if (ia > o2) exit
                    s = s2(a1,ia)
            print "('    n: ',I2.1,'    k: ',I2.1,'    S(n,k): ',I7.1)",&
     &              a1,ia,s
                    ia=ia+1
                end do
                a1=a1+1
            end do
        end subroutine stirling_numbers


        recursive function s2(n,k) result(s2_result)
            use constants_module
            integer(kind=large_int_kind) :: n,k,s2_result
            if (n .eq. 0 .and. k .eq. 0) then
                s2_result = 1
            else if (n .eq. 0 .or. k .eq. 0) then
                s2_result = 0
            else
                s2_result = s2(n-1,k-1) + (k * s2(n-1,k))
            end if
        end function s2



        !!
        !!  subroutine : CHI_SQUARE_DEMO
        !!
        !!      sets up the chi-square quantiles demonstration and governs
        !!      the invocation of 'get_chi_quantiles'
        !!
        subroutine chi_square_demo
            use utilities_module
            print *, "    demo #6: chi-square quantiles"
            print *, "    _____________________________"
            print *, " "
            print *, " "
            print *, "    This program calculates the critical values"
            print *, "    of the chi-square distribution for degrees of"
            print *, "    freedom 10 through 200, with an increment of"
            print *, "    10, at significance level 0.90.  The code is"
            print *, "    based on ACM algorithm #451."
            print *, " "
            call mypause
            call system("cls")
            print *, " "
            call get_chi_quantiles
            print *, " "
            call mypause
        end subroutine chi_square_demO


        subroutine get_chi_quantiles
            use constants_module
            real(kind=large_real_kind) :: pa,csqa,chisqd
            integer(kind=large_int_kind) :: df,omega_df
            pa = 0.1
            df = 10
            omega_df = 200
            do 
                if (df > omega_df) exit
                csqa = chisqd(pa,df)
                print*,df,csqa
                df=df+10
            end do
        end subroutine get_chi_quantiles


        function chisqd(p,n)
            use constants_module
            real(kind=large_real_kind) :: chisqd
            real(kind=large_real_kind) :: p,gaussd,f,f1,f2,t
            integer(kind=large_int_kind) :: n
            real(kind=large_real_kind) :: c(21),a(19)
            c = (/                                                      &
     &           1.565326e-3,                                           &  ! c(01) 
     &           1.060438e-3,                                           &  ! c(02) 
     &          -6.950356e-3,                                           &  ! c(03) 
     &          -1.323293e-2,                                           &  ! c(04) 
     &           2.277679e-2,                                           &  ! c(05) 
     &          -8.986007e-3,                                           &  ! c(06) 
     &          -1.513904e-2,                                           &  ! c(07) 
     &           2.530010e-3,                                           &  ! c(08) 
     &          -1.450117e-3,                                           &  ! c(09) 
     &           5.169654e-3,                                           &  ! c(10) 
     &          -1.153761e-2,                                           &  ! c(11) 
     &           1.128186e-2,                                           &  ! c(12) 
     &           2.607083e-2,                                           &  ! c(13) 
     &          -0.2237368  ,                                           &  ! c(14) 
     &           9.780499e-5,                                           &  ! c(15) 
     &          -8.426812e-4,                                           &  ! c(16) 
     &           3.125580e-3,                                           &  ! c(17) 
     &          -8.553069e-3,                                           &  ! c(18) 
     &           1.348028e-4,                                           &  ! c(19) 
     &           0.4713941  ,                                           &  ! c(20) 
     &           1.0000886  /)                                             ! c(21) 

            a =  (/                                                     &
     &            1.264616e-2,                                          &  ! a(01)
     &           -1.425296e-2,                                          &  ! a(02)
     &            1.400483e-2,                                          &  ! a(03)
     &           -5.886090e-3,                                          &  ! a(04)
     &           -1.091214e-2,                                          &  ! a(05)
     &           -2.304527e-2,                                          &  ! a(06)
     &            3.135411e-3,                                          &  ! a(07)
     &           -2.728494e-4,                                          &  ! a(08)
     &           -9.699681e-3,                                          &  ! a(09)
     &            1.316872e-2,                                          &  ! a(10)
     &            2.618914e-2,                                          &  ! a(11)
     &           -0.2222222  ,                                          &  ! a(12)
     &            5.406674e-5,                                          &  ! a(13)
     &            3.483789e-5,                                          &  ! a(14)
     &           -7.274761e-4,                                          &  ! a(15)
     &            3.292181e-3,                                          &  ! a(16)
     &           -8.729713e-3,                                          &  ! a(17)
     &            0.4714045  ,                                          &  ! a(18)
     &            1.0000000  /)                                            ! a(19)

            if (n-2 < 0) then
                chisqd = gaussd(0.5 * p)
                chisqd = chisqd * chisqd
                return
            else if (n-2 .eq. 0) then
                chisqd = -2.0 * log(p)
                return
            else
                f = real(n)
                f1 = 1.0 / f
                t = gaussd(p)
                f2 = sqrt(f1) * t
                if (n .ge. (2 + int(4.0 * abs(t)))) then
                    chisqd = (((a(1) + a(2) * f2) * f1 +(((a(3) +       &
     &                      a(4) * f2) * f2                             &
     &                      + a(5)) * f2 + a(6))) * f1 + (((((          &
     &                      a(7) + a(8) * f2) * f2 + a(9)) * f2         &
     &                      + a(10)) * f2 + a(11)) * f2 + a(12)))       &
     &                      * f1 + (((((a(13) * f2                      &
     &                      + a(14)) * f2 + a(15)) * f2 + a(16)) * f2   &
     &                      + a(17)) * f2 * f2                          &
     &                      + a(18)) * f2 + a(19)
                    chisqd = chisqd * chisqd * chisqd * f
                    return
                else
                    chisqd = (((((((c(1) * f2 + c(2)) * f2 + c(3))      &
     &                      * f2 + c(4)) * f2                           &
     &                      + c(5)) * f2 + c(6)) * f2 + c(7)) * f1      &
     &                      + ((((((c(8) + c(9) * f2) * f2              &
     &                      + c(10)) * f2 + c(11)) * f2 + c(12)) *      &
     &                      f2 + c(13)) * f2 + c(14))) * f1 +           &
     &                      (((((c(15) * f2 + c(16)) * f2 + c(17))      &
     &                      * f2 + c(18)) * f2                          &
     &                      + c(19)) * f2 + c(20)) * f2 + c(21)
                    chisqd = chisqd * chisqd * chisqd * f
                    return
                end if
            end if
        end function chisqd


        function gaussd(prob)
            use constants_module
            real(kind=large_real_kind) :: gaussd
            real(kind=large_real_kind) :: prob,factor,acc
            factor = 1000.0
            gaussd = acc(prob*factor,factor)
        end function gaussd


        function acc(tar,fct)
            use constants_module
            real(kind=large_real_kind) :: acc
            real(kind=large_real_kind) :: tar,act,fct,f,e,val
            act = 0.0
            val = -5.0 - (1.0/fct)
            f = 1 / sqrt(2 * 3.14159265358979323846)
            do 
                if (act >= tar) exit
                val = val + (1.0/fct)
                e = - ((val*val) / 2)
                act = act + (f * exp(e))
            end do
            acc = val
        end function acc



        !!
        !!  subroutine : PYTHAGORAS_DEMO
        !!
        !!      sets up the Pythagorean triplets demonstration and governs
        !!      the invocation of 'gen_triplets'
        !!
        subroutine pythagoras_demo
            use utilities_module
            print *, "    demo #7: Pythagorean triplets"
            print *, "    _____________________________"
            print *, " "
            print *, " "
            print *, "    A Pythagorean triplet is a set of 3 integers,"
            print *, "    {x,y,z}, that satisfies the formula x*x + y*y"
            print *, "    = z*z.   All primitive Pythagorean triplets"
            print *, "    can be derived from the formulas"
            print *, " "
            print *, "                 x = 2uv"
            print *, "                 y = (u**2) - (v**2)"
            print *, "                 z = (u**2) + (v**2)"
            print *, " "
            print *, "    where u > v, u and v have no common factor,"
            print *, "    and one of them is odd and the other even."
            print *, "    These formulas are presented geometrically in"
            print *, "    the 10th book of Euclid's Elements."
            print *, " "
            print *, "    This program prints out all of the primitive"
            print *, "    Pythagorean triplets for values of u and v"
            print *, "    from 1 through 9."
            print *, " "
            call mypause
            call system("cls")
            print *, " "
            call gen_triplets(1,9)
            print *, " "
            call mypause
        end subroutine pythagoras_demo


        recursive subroutine gen_triplets(ialpha,omega)
            use constants_module
            integer(kind=large_int_kind) :: alpha,omega,ialpha
            alpha=ialpha
            if (alpha <= omega) then
                call cycle_over_v(1,alpha)
                call gen_triplets(alpha+1,omega)
            end if
        end subroutine gen_triplets


        recursive subroutine cycle_over_v(v,u)
            use constants_module
            integer(kind=large_int_kind) :: u,v,gcd
            logical :: o_e
            if (v < u) then
                if ((gcd(u,v) .eq. 1) .and. (o_e(u,v))) then
                    call print_triplet(v,u)
                end if
                call cycle_over_v(v+1,u)
            end if
        end subroutine cycle_over_v


        logical function o_e(u,v)
            use constants_module
            integer(kind=large_int_kind) :: u,v,mu,mv
            mu = mod(u,2)
            mv = mod(v,2)
            if ((mu .eq. 0) .and. (mv .eq. 0)) then
                o_e = .false.
            else if ((mu .ne. 0) .and. (mv .ne. 0)) then
                o_e = .false.
            else
                o_e = .true.
            end if
        end function o_e


        recursive function gcd(x,y) result(gcd_result)
            use constants_module
            integer(kind=large_int_kind) :: x,y,gcd_result
            if (0 .eq. y) then
                gcd_result = x
            else
                gcd_result = gcd(y,mod(x,y))
            end if
        end function gcd


        subroutine print_triplet(v,u)
            use constants_module
            integer(kind=large_int_kind) :: u,v,x,y,z
            x = 2*u*v
            y = (u*u) - (v*v)
            z = (u*u) + (v*v)
            print "('     u: ',I2.1,'    v: ',I2.1,'    triplet: {',    &
     &              I5.1,',',I5.1,',',I5.1,'}')",                       &
     &              u,v,x,y,z
        end subroutine print_triplet



        !!
        !!  subroutine : RUNTIME_FUNCS_DEMO
        !!
        !!      runs a miscellaneous collection of runtime calls
        !!
        subroutine runtime_funcs_demo
          use constants_module
          use utilities_module
          character (len=10) :: time_result
          character (len=8) :: date_result
          integer(kind=large_int_kind) :: at, ot, f, i
          integer(kind=large_int_kind) :: factorial, fact_iter
            print *, "    demo #8: date, time, and system clock calls"
            print *, "    ___________________________________________"
            print *, " "
            print *, " "
            print *, "   This program exercises the following calls:"
            print *, " "
            print *, "               call date_and_time,"
            print *, "               call system_clock, and"
            print *, "               call system."
            print *, " "
            call mypause
            call system("cls")
            call date_and_time (date=date_result,time=time_result)
            print *, "  The current time is: ",time_result(1:2),":",    &
     &time_result(3:4),":",time_result(5:9)
            print *, "  The current date is: ",date_result(5:6),"/",    &
     &date_result(7:8),"/",date_result(3:4)
            print *, " "
            print *, " "

            i = 100000
            call system_clock(at)
            do 
                if (i <= 0) exit
                f = factorial(12)
                i = i - 1
            end do
            call system_clock(ot)
            print *, "  The number of ticks it takes to recursively"
            print *, "  calculate and store 'factorial(12)' 100000"
            print *, "  times is:", ot-at
            i = 100000
            call system_clock(at)
            do 
                if (i <= 0) exit
                f = fact_iter(12)
                i = i - 1
            end do
            call system_clock(ot)
            print *, " "
            print *, "  The number of ticks it takes to iteratively"
            print *, "  calculate and store 'factorial(12)' 100000"
            print *, "  times is:", ot-at
            print *, " "
            print *, "  The value of factorial(12) is:",f
            print *, " "
            print *, " "
            print *, "  The following quotation is printed with system"
            print *, "  calls:"
            print *, " "
            print *, " "
            call system('echo       High thoughts must have high language.')
            call system('echo                  -- Aristophanes, The Frogs')
            print *, " "
            call mypause
        end subroutine runtime_funcs_demo


        function fact_iter(n)
          use constants_module
          integer(kind=large_int_kind) :: fact_iter
          integer(kind=large_int_kind) :: n, i, r
            r = 1
            i = n
            do 
                if (i <= 0) exit
                r = r * i
                i = i - 1
            end do
            fact_iter = r
        end function fact_iter