!!---------------------------------------------------------------
!!
!!           file:  demo.f90
!!
!!        version:  Elf90 1.00
!!
!!           date:  February 1995
!!
!!        purpose:  This Fortran source file tests a subset of the features
!!                  of the Lahey Elf90 compiler and its runtime support.
!!                  Successful execution of this program indicates that
!!                  installation has been successful.
!!
!!  compatibility:  This Fortran program can be compiled and run with all Lahey
!!                  Elf90 and ALL Fortran 90 Language Systems
!!
!!  to create/run:  Enter the following two commands in the order shown:
!!
!!                      [1]  elf90 demo<CR>
!!                      [2]  demo<CR>
!!
!!    portability:  This code is portable to any system running Fortran 90
!!
!!    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.
!!
!!----------------------------------------------------------------------------
!!
!! Unpublished - rights reserved under the Copyright Laws of the United States.
!! Use, duplication, or disclosure by the Government is subject to restrictions
!! as set forth in subparagraph (c)(1)(ii) of the Rights in Technical Data and
!! Computer Software clause at 252.227-7013.
!!
!!                       Lahey Computer Systems, Inc.
!!                                PO Box 6091
!!                         Incline Village, NV 89450
!!                            voice: 702-831-2500
!!                              fax: 702-831-8123
!!                        e-mail: support@lahey.com
!!
!!              Copyright(c) 1995-97, Lahey Computer Systems, Inc.
!!
!!----------------------------------------------------------------------------

        !! module constants_module

        !!  contains:

        !!      parameters for users to declare integer, real, and
        !!      complex types using KIND=

        !!      Useful constants
module constants_module
  implicit none
  integer, parameter                     ::  large_int_kind  = 4,       &
                                             large_real_kind = 8
  real(kind=large_real_kind), parameter  ::                             &
                        pi = 3.14159265358979323846_large_real_kind,    &
                        two_pi = 2*pi
end module constants_module


        !! module utilities_module
        !! for doing little common things
module utilities_module
  implicit none
  private
  public  ::  my_pause, factorial

  contains

    subroutine my_pause()
      write(*,*) "Press the Enter key to continue"
      read (*, *)
      return
    end subroutine my_pause

    recursive function factorial(n) result(fact_result)
      use constants_module
      integer(kind=large_int_kind)              ::  fact_result
      integer(kind=large_int_kind), intent(in)  ::  n

      if (n == 1) then
        fact_result = 1
      else
        fact_result = n * factorial( n-1 )
      end if
      return
    end function factorial

end module utilities_module


        !! module external_sub_procs
        !! contains the 8 demo subroutines and their internal procedures
module external_sub_procs
  use constants_module
  use utilities_module

  implicit none
  private
  public  ::  factorial_demo,  conversion_demo,   carmichael_demo
  public  ::  ramanujan_demo,  stirling_demo,     chi_square_demo
  public  ::  pythagoras_demo, runtime_funcs_demo

  contains

        !!  subroutine:  FACTORIAL_DEMO
        !!      sets up the factorial demonstration and governs the
        !!      invocation of 'factorial'
    subroutine factorial_demo()
      integer(kind=large_int_kind) :: i

      write (*,*) "    demo #1: factorials"
      write (*,*) "    ___________________"
      write (*,*)
      write (*,*)
      write (*,*) "    Calculate and print out the factorials"
      write (*,*) "    of the integers 1 through 12."
      write (*,*)
      call my_pause()
      write(*,"(/////////////////////////)")     !  Clear the screen
      write (*,*)
      do i = 1, 12
        write(*,"('     n = ',I2.1,                              &
     &            '      factorial(n) = ',I9.1)") i, factorial(i)
      end do

      write (*,*)
      call my_pause()
      return
    end subroutine factorial_demo


        !!  subroutine:  CONVERSION_DEMO

        !!      sets up the Fahrenheit to Celsius conversion demonstration
    subroutine conversion_demo()
      integer(kind=large_int_kind), parameter  ::  l=5, h=100, i=5
      integer(kind=large_int_kind)             ::  fahr
      real(kind=large_real_kind)               ::  celsius

      write (*,*) "    demo #2: Fahrenheit to Celsius conversion"
      write (*,*) "    _________________________________________"
      write (*,*)
      write (*,*)
      write (*,*) "    This routine calculates and prints out the"
      write (*,*) "    Celsius equivalent of every fifth integral"
      write (*,*) "    Fahrenheit value from 5 to 100 degrees."
      write (*,*)
      call my_pause()
      write (*,"(/////////////////////////)")     !      Clear the screen
      write (*,*) '     deg F      deg C'
      write (*,*) '    -------    -------'
      do fahr = l, h, i
        celsius = (5.0/9.0) * (fahr - 32.0)
        write (*,"('     ',F7.2,'    ',F7.2)") real(fahr), celsius
      end do

      call my_pause()
      return
    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
        !!      Computer Programs" provided the inspiration for this example
    subroutine carmichael_demo()
      write (*,*) "    demo #3: Carmichael numbers"
      write (*,*) "    ___________________________"
      write (*,*)
      write (*,*)
      write (*,*) "    A Carmichael number is a non-prime that fools"
      write (*,*) "    the Fermat test for primality.  These numbers"
      write (*,*) "    are rare.  This program examines the integers"
      write (*,*) "    from 2 to 3,000 and lists out the Carmichael"
      write (*,*) "    numbers discovered."
      write (*,*)
      call my_pause()
      write (*,"(/////////////////////////)")   !      Clear the screen
      write (*,*)
      call carmichaels_in_range(2, 3000)
      write (*,*)
      call my_pause()
      return
    end subroutine carmichael_demo


    subroutine carmichaels_in_range(alpha, omega)
      integer(kind=large_int_kind), intent(in)  ::  alpha, omega
      integer(kind=large_int_kind)              ::  ialpha

      do  ialpha = alpha, omega
        if ( fast_prime(ialpha, 9)  .and.  .not. naive_test(ialpha) )      &
          write (*,*) ialpha, " is a Carmichael number!"
      end do
      return

      contains
        recursive function fast_prime(n, times) result(fast_result)
          integer(kind=large_int_kind), intent(in)  ::  n, times
          logical                                   ::  fast_result

          fast_result = ( times == 0 )
          if ( .not. fast_result )  then
            if ( fermat(n)  .eqv.  .true. )            &
                        fast_result = fast_prime(n, times-1)
            end if
          return
        end function fast_prime


        function fermat(n)
          logical                                   ::  fermat
          integer(kind=large_int_kind), intent(in)  ::  n
          integer(kind=large_int_kind)              ::  alpha
          real                                      ::  rnd

          call random_number(rnd)
          alpha = 2 + ( (n - 2) * rnd )
          fermat = ( alpha == expmod(alpha, n, n) )
          return
        end function fermat


        recursive function expmod(b, e, m) result(exp_result)
          integer(kind=large_int_kind)              ::  exp_result
          integer(kind=large_int_kind), intent(in)  ::  b, e, m
          integer(kind=large_int_kind)              ::  tmp

          if ( e == 0 ) then
            exp_result = 1
          else if ( 0 == mod(e, 2) ) then!             If e is even
            tmp = expmod(b, e/2, m)
            exp_result = mod(tmp*tmp, m)
          else!                                        If e is odd
            tmp = expmod(b, e-1, m)
            exp_result = mod(b*tmp, m)
          end if
          return
        end function expmod


        function naive_test(n)
          logical                                    ::   naive_test
          integer(kind=large_int_kind), intent(in)   ::   n

          naive_test = ( n == find_divisor(n, 2) )
          return
        end function naive_test


        recursive function find_divisor(n, t) result(divisor)
          integer(kind=large_int_kind)              ::  divisor
          integer(kind=large_int_kind), intent(in)  ::  n, t

          if ( t**2 > n) then
            divisor = n
          else if ( 0 == mod(n, t) ) then
            divisor = t
          else
            divisor = find_divisor(n, t+1)
          end if
          return
        end function find_divisor

      end subroutine carmichaels_in_range


        !!  subroutine:  RAMANUJAN_DEMO

        !!      sets up the Ramanujan series demonstration and
        !!      governs the invocation of 'search_ramanujans'

        !!      Abelson & Sussman's "Structure and Interpretation of
        !!      Computer Programs" provided the inspiration for this example
      subroutine ramanujan_demo()

        write (*,*) "    demo #4: Ramanujan's series"
        write (*,*) "    ___________________________"
        write (*,*)
        write (*,*)
        write (*,*) "    G. H. Hardy, in his obituary for Srinavasa"
        write (*,*) "    Ramanujan, told of the time that he visited"
        write (*,*) "    Ramanujan when the latter was ill.  Hardy "
        write (*,*) "    took taxi number '1729' to where Ramanujan"
        write (*,*) "    was staying, and in passing told Ramanujan"
        write (*,*) "    that it seemed a very uninteresting number."
        write (*,*) "    Ramanujan instantly took issue with Hardy,"
        write (*,*) "    pointing out that it is on the contrary"
        write (*,*) "    very interesting,  because it is the smallest"
        write (*,*) "    positive integer that can be expressed as the"
        write (*,*) "    sum of 2 cubes in exactly 2 different ways."
        write (*,*) "    The number '1729' has since been known as"
        write (*,*) "    Ramanujan's number and the series it begins"
        write (*,*) "    as Ramanujan's series."
        write (*,*)
        write (*,*) "    This program prints out the Ramanujan numbers"
        write (*,*) "    discovered between 1 and 15,000."
        write (*,*)
        call my_pause()
        write (*,"(/////////////////////////)")   !       Clear the screen
        write (*,*)
        call search_ramanujans(1, 15000)
        write (*,*)
        call my_pause()
        return
      end subroutine ramanujan_demo


      subroutine search_ramanujans(alpha, omega)
        integer(kind=large_int_kind), intent(in)  ::  alpha, omega
        integer(kind=large_int_kind)              ::  n

        do  n = alpha, omega
          if ( 0 == mod(n, 1000) )  write (*,*) "at:", n
          if ( is_ramanujan(n) )    write (*,*) n, " is a Ramanujan number!"
        end do
        return
      end subroutine search_ramanujans


      function is_ramanujan(n)
        logical                                   ::  is_ramanujan
        integer(kind=large_int_kind), intent(in)  ::  n
        integer(kind=large_int_kind)              ::  lim

        lim = 1 + int( (n/2)**(1.0/3.0) )
        is_ramanujan = search(n, 1, lim, 0)
        return
      end function is_ramanujan


      recursive function search(n, a, o, cubes) result(searchres)
        integer(kind=large_int_kind), intent(in)  ::  n, a, o, cubes
        integer(kind=large_int_kind)              ::  ca, re, rc
        logical                                   ::  searchres

        if (a >= o) then
          searchres = ( 2 == cubes )
        else
          ca = a**3
          re = n - ca
          rc = int( re**(1.0/3.0) )
          if ( rc**3 == re ) then
            searchres = search(n, a+1, o, cubes+1)
          else
            searchres = search(n, a+1, o, cubes)
          end if
        end if
        return
      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()
        write (*,*) "    demo #5: Stirling numbers of the 2nd kind"
        write (*,*) "    _________________________________________"
        write (*,*)
        write (*,*)
        write (*,*) "    Stirling numbers of the 2nd kind, S(n,k),"
        write (*,*) "    are the numbers of possible partitions of"
        write (*,*) "    n items into k groups.  They are useful in"
        write (*,*) "    calculating the expected probability of"
        write (*,*) "    some kinds of combinations of events."
        write (*,*)
        write (*,*) "    This program prints out Stirling numbers"
        write (*,*) "    of the 2nd kind from S(12,0) to S(12,12)."
        write (*,*)
        call my_pause()
        write (*,"(/////////////////////////)")   !        Clear the screen
        write (*,*)
        call stirling_numbers(12, 12, 0, 12)
        write (*,*)
        call my_pause()
        return
      end subroutine stirling_demo


      subroutine stirling_numbers(ia1, io1, ia2, io2)
        integer(kind=large_int_kind), intent(in)  ::  ia1, io1, ia2, io2
        integer(kind=large_int_kind)              ::  a1, ia
   
        do  a1 = ia1, io1
          do  ia = ia2, io2
            write (*,"('    n: ',I2.1,'    k: ',I2.1,'    S(n,k):',I7.1)") &
                     a1, ia, s2(a1, ia)
          end do
        end do
        return
      end subroutine stirling_numbers


      recursive function s2(n, k) result(s2_result)
        integer(kind=large_int_kind), intent(in)  ::  n, k
        integer(kind=large_int_kind)              ::  s2_result

        if ( n == 0  .and.  k == 0 ) then
          s2_result = 1
        else if ( n == 0  .or.  k == 0 ) then
          s2_result = 0
        else
          s2_result = s2(n-1, k-1) + (k * s2(n-1, k))
        end if
        return
      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()
        write (*,*) "    demo #6: chi-square quantiles"
        write (*,*) "    _____________________________"
        write (*,*)
        write (*,*)
        write (*,*) "    This program calculates the critical values"
        write (*,*) "    of the chi-square distribution for degrees of"
        write (*,*) "    freedom 10 through 200, with an increment of"
        write (*,*) "    10, at significance level 0.90.  The code is"
        write (*,*) "    based on ACM algorithm #451."
        write (*,*)
        call my_pause()
        write (*,"(/////////////////////////)")   !      Clear the screen
        write (*,*)
        call get_chi_quantiles()
        write (*,*)
        call my_pause()
        return
      end subroutine chi_square_demO


      subroutine get_chi_quantiles()
        real(kind=large_real_kind),   parameter  ::  pa=0.1
        integer(kind=large_int_kind), parameter  ::  df=10, omega_df=200
        integer(kind=large_int_kind)             ::  n

        do  n = df, omega_df, df
          write (*,*) n, chisqd(pa, n)
        end do
        return
      end subroutine get_chi_quantiles


      function chisqd(p, n)
        real(kind=large_real_kind)                ::  chisqd
        real(kind=large_real_kind), intent(in)    ::  p
        real(kind=large_real_kind)                ::  f, f1, f2, t
        integer(kind=large_int_kind), intent(in)  ::  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
        else if (n-2 == 0) then
          chisqd = -2.0 * log(p)
        else
          f = real(n)
          f1 = 1.0 / f
          t = gaussd(p)
          f2 = t * sqrt(f1)
          if ( n >= (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 = f * chisqd**3
          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 = f * chisqd**3
          end if
        end if
        return
      end function chisqd


      function gaussd(prob)
        real(kind=large_real_kind)              ::  gaussd
        real(kind=large_real_kind), intent(in)  ::  prob
        real(kind=large_real_kind), parameter   ::  factor = 1000.0
     
        gaussd = acc(prob*factor, factor)
        return
      end function gaussd


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


        !!  subroutine:  PYTHAGORAS_DEMO

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


      recursive subroutine gen_triplets(alpha, omega)
        integer(kind=large_int_kind), intent(in)  ::  alpha, omega

        if (alpha <= omega) then
          call cycle_over_v(v=1, u=alpha)
          call gen_triplets(alpha+1, omega)
        end if
        return
      end subroutine gen_triplets


      recursive subroutine cycle_over_v(v, u)
        integer(kind=large_int_kind), intent(in)  ::  u, v

        if (v < u) then
          if ((gcd(u, v) == 1)  .and.  (odd_even(u, v)))                     &
            call print_triplet(v, u)
          call cycle_over_v(v+1, u)
        end if
        return
      end subroutine cycle_over_v


      function odd_even(u, v)!              Both even or both odd then false
        logical                                   ::  odd_even
        integer(kind=large_int_kind), intent(in)  ::  u, v

        odd_even = ( mod(u+v, 2) == 1 )!   Sum even if both odd or both even
        return
      end function odd_even


      recursive function gcd(x, y) result(gcd_result)
        integer(kind=large_int_kind)              ::  gcd_result
        integer(kind=large_int_kind), intent(in)  ::  x, y
   
        if ( y == 0 ) then
          gcd_result = x
        else
          gcd_result = gcd(y, mod(x,y))
        end if
        return
      end function gcd


      subroutine print_triplet(v, u)
        integer(kind=large_int_kind), intent(in)  ::  u, v

        write (*,"('     u: ',I2.1,'    v: ',I2.1,'    triplet: {',               &
     &           I5.1,',',I5.1,',',I5.1,'}')")                                 &
                 u, v, 2*u*v, (u*u) - (v*v), (u*u) + (v*v)
        return

      end subroutine print_triplet


        !!  subroutine:  RUNTIME_FUNCS_DEMO

        !!      runs a miscellaneous collection of runtime calls
      subroutine runtime_funcs_demo()
        character (len=10)            ::  time_result
        character (len=8)             ::  date_result
        integer(kind=large_int_kind)  ::  at, ot, f, i

        write (*,*) "    demo #8: date, time, and system clock calls"
        write (*,*) "    ___________________________________________"
        write (*,*)
        write (*,*)
        write (*,*) "   This program exercises the following calls:"
        write (*,*)
        write (*,*) "               call date_and_time,"
        write (*,*) "               call system_clock"
        write (*,*)
        call my_pause()
        write (*,"(/////////////////////////)")   !        Clear the screen
        call date_and_time (date=date_result,time=time_result)
        write (*,*) "  The current time is:", time_result(1:2), ":",              &
                                           time_result(3:4), ":",              &
                                           time_result(5:9)
        write (*,*) "  The current date is: ", date_result(5:6), "/",             &
                                            date_result(7:8), "/",             &
                                            date_result(3:4)
        write (*,*)
        write (*,*)
   
        call system_clock(at)
        do  i = 100000, 0, -1
          f = factorial(12)
        end do
        call system_clock(ot)

        write (*,*) "  The number of ticks it takes to recursively"
        write (*,*) "  calculate and store 'factorial(12)' 100000"
        write (*,*) "  times is:", ot-at
    
        call system_clock(at)
        do  i = 100000, 0, -1
          f = fact_iter(12)
        end do
        call system_clock(ot)
  
        write (*,*)
        write (*,*) "  The number of ticks it takes to iteratively"
        write (*,*) "  calculate and store 'factorial(12)' 100000"
        write (*,*) "  times is:", ot-at
        write (*,*)
        write (*,*) "  The value of factorial(12) is:", f
        write (*,*)
        write (*,*)
        write (*,*)
        write (*,*)
        write (*,*) '            High thoughts must have high language.'
        write (*,*) '                       -- Aristophanes, The Frogs'
        write (*,*)
        call my_pause()
        return
      end subroutine runtime_funcs_demo


      function fact_iter(n)
        integer(kind=large_int_kind)              ::  fact_iter
        integer(kind=large_int_kind), intent(in)  ::  n
        integer(kind=large_int_kind)              ::  i

        fact_iter = 1

        do  i = n, 2, -1
          fact_iter = fact_iter * i
        end do
        return
      end function fact_iter

end module external_sub_procs


        !!   main program unit:  LAHEY_DEMO

        !!      sets up and governs the selection between the
        !!      demonstration subprograms
program lahey_demo
  use constants_module
  use utilities_module
  use external_sub_procs

  implicit none
  integer :: selection

  write (*,"(/////////////////////////)")   !          Clear the screen
  write (*,*) "             Lahey Elf90 Compiler"
  write (*,*) "             --------------------"
  write (*,*)
  write (*,*) "  installation test and demonstration program"
  write (*,*)
  write (*,*) "              Copyright(c) 1995-96"
  write (*,*) "         Lahey Computer Systems, Inc."
  do
    write (*,*)
    write (*,*) "    -----------------"
    write (*,*) "    Test/Action List:"
    write (*,*) "    -----------------"
    write (*,*) "     1 - factorials"
    write (*,*) "     2 - Fahrenheit to Celsius conversion"
    write (*,*) "     3 - Carmichael numbers"
    write (*,*) "     4 - Ramanujan's series"
    write (*,*) "     5 - Stirling numbers of the 2nd kind"
    write (*,*) "     6 - chi-square quantiles"
    write (*,*) "     7 - Pythagorean triplets"
    write (*,*) "     8 - date_and_time, and system_clock"
    write (*,*) "     0 - <stop this program>"
    write (*,*)
    write (*,*) "    Please select an option by entering the"
    write (*,*) "    associated number followed by <enter>.    "
    read (*,*) selection
    write (*,"(/////////////////////////)")   !       Clear the screen
    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()
      case default
        write (*,*) "Digits must be in the closed interval 0:8"
      call my_pause()
    end select
    write (*,"(/////////////////)")   !              Clear the screen
  end do
  stop
end program lahey_demo
