!         Author: Jordan Schettler
! Date Submitted: 4/26/05
!        Purpose: Final Project
!
MODULE primality_tests

  TYPE :: rep
    INTEGER(kind=8) :: s, t
  END TYPE rep

  CONTAINS

! Recursive powering algorithm with modular reduction;
! we compute a^n (mod m) for positive integers a, n, m, under
! the assumption that 1 < a < m.

  RECURSIVE FUNCTION pow(a,n,m) RESULT(r)
    IMPLICIT none
    INTEGER(kind=8), INTENT(in) :: a, n, m
    INTEGER(kind=8) :: r

    IF (n == 1) THEN
      r = a
    ELSE IF (MOD(n, 2) == 0) THEN
      r = MOD(pow(a, n/2, m)**2, m)
    ELSE
      r = MOD(a*(pow(a, (n - 1)/2, m)**2), m)
    END IF
  END FUNCTION pow

! Decomposition function for odd n > 1;
! we find the unique representation n = 1 + (2^r%s)r%t
! where r%t is odd.

  TYPE(rep) FUNCTION decomp(n) RESULT(r)
    IMPLICIT none
    INTEGER(kind=8), INTENT(in) :: n

    r%s = 1
    r%t = (n - 1)/2
    DO WHILE(MOD(r%t, 2) == 0)
      r%t = r%t/2
      r%s = r%s + 1
    END DO
  END FUNCTION decomp

! Strong probable prime test; given a positive integer n and an integer a, we determine
! if n is a strong probable prime base a.

  LOGICAL FUNCTION is_strong(n, a) RESULT(r)
    IMPLICIT none
    INTEGER(kind=8), INTENT(in) :: n, a
    TYPE(rep) :: D
    INTEGER(kind=8) :: i, b

    D = decomp(n)
    b = pow(a, D%t, n)
    r = .FALSE.
    
    IF (b == 1 .OR. b == n - 1) THEN
      r = .TRUE.
    ELSE
      even_part: DO i = 1, D%s - 1
        b = MOD(b**2, n)
        IF (b == n - 1) THEN
          r = .TRUE.
          EXIT even_part
        END IF
      END DO even_part
    END IF
  END FUNCTION is_strong

! Miller's primality test. Given an odd integer n > 3, we check if n is a strong probable prime
! to each base within a certain range. If not, then n is definitely composite; if so, then n
! is either prime or the extended Riemann hypothesis is false. Here we assume the ERH is true.

  LOGICAL FUNCTION is_prime(n) RESULT(r)
    IMPLICIT none
    INTEGER(kind=8), INTENT(in) :: n
    INTEGER(kind=8) :: i, W

    W = MIN(INT(2*(LOG(REAL(n))**2)), n - 1)
    r = .TRUE.

    comp_test: DO i = 2, W
      IF (.NOT. is_strong(n, i)) THEN
        r = .FALSE.
        EXIT comp_test
      END IF
    END DO comp_test
  END FUNCTION is_prime

! Full reptend prime test. Given a prime p, we determine if 10 is a primitive root mod p.
! That is, 10^k is congruent to 1 for k = p - 1 = phi(p) mod p, while this does not hold
! for any smaller choice of k.

  LOGICAL FUNCTION is_full_reptend(p) RESULT(r)
    IMPLICIT none
    INTEGER(kind=8), INTENT(in) :: p
    INTEGER(kind=8) :: i

    r = .FALSE.
    IF (pow(10, p - 1, p) == 1) THEN
      r = .TRUE.
      check: DO i = 1, p - 2
        IF (pow(10, i, p) == 1) THEN
          r = .FALSE.
          EXIT check
        END IF
      END DO check
    END IF
  END FUNCTION is_full_reptend
END MODULE primality_tests

PROGRAM main
  USE primality_tests
  IMPLICIT none
  INTEGER(kind=8) :: i, j, k, status
  INTEGER, PARAMETER :: unit = 10
  REAL :: x
  CHARACTER(len=*), PARAMETER :: filename = "full_reptend_primes"

  OPEN (UNIT = unit, FILE = filename, STATUS = 'REPLACE', ACTION = 'WRITE', IOSTAT = status)
  IF (status /= 0) THEN
    WRITE(*,*) "Oops!"
    STOP
  END IF

  WRITE(*, '("Enter a real number x to create or replace a file named full_reptend_primes&
  & containing a list (in increasing order) of full reptend primes which are less than or equal&
  & to x: ")', ADVANCE = 'NO')
  READ(*,*) x

! We need only check for odd integers > 3 since that is the lower bound of Miller's test,
! and clearly 1/2 = .5 and 1/3 = .333... do not have maximal periods in their decimal exapansions.

  k = 0
  j = 0
  DO i = 5, FLOOR(x), 2
    IF (is_prime(i)) THEN
      j = j + 1
      IF (is_full_reptend(i)) THEN
        WRITE(unit, '(I)') i
        k = k + 1
      END IF
    END IF
  END DO

!  WRITE(*,*) "An approximation to Artin's constant is ", k/REAL(j + 2)

  CLOSE(UNIT = unit)

END PROGRAM main
