PROGRAM Main
    USE OMP_LIB
    IMPLICIT NONE
    INTEGER, PARAMETER :: big = SELECTED_INT_KIND(18), nprimes = 100
    INTEGER :: threads, count, space, width, k, n
    CHARACTER(LEN=100) :: buffer
    INTEGER(KIND=big) :: limit
    INTEGER, DIMENSION(nprimes) :: primes
    INTEGER(KIND=big), ALLOCATABLE, DIMENSION(:) :: numbers
    REAL(KIND=KIND(0.0D0)) :: temp

    CALL Erastothenes(primes)
    IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN
        WRITE (*,'("Wrong number of arguments")')
        STOP 1
    END IF
    CALL GET_COMMAND_ARGUMENT(1,buffer)
    READ (buffer,"(I100)",IOSTAT=k) limit
    IF (k /= 0 .OR. limit < 10) THEN
        WRITE (*,'("Invalid limiting value: ",A)') buffer
        STOP 1
    END IF
    ! Allocate a big enough array - see Wikipedia
    temp = LOG(REAL(limit,KIND=KIND(0.0D0)))
    space = INT(limit*EXP(-1.8D0*temp*LOG(LOG(temp))/LOG(temp)))
    ALLOCATE (numbers(space))
    count = 0
! The OpenMP clauses are utterly obscene, but seem to be what is
! required, though leaving them out would also work.
!$OMP PARALLEL DEFAULT(NONE), SHARED(threads,limit,k), PRIVATE(n)
!$OMP SINGLE
    threads = OMP_GET_NUM_THREADS()
    IF (threads == 1) THEN
        CALL Search(3_big,limit,2_big)
    ELSE
        threads = 9*threads/10
        DO n = 1,threads
        k = n    ! SHARED loop variables are often rejected
!$OMP TASK DEFAULT(NONE), FIRSTPRIVATE(threads,limit,k)
            CALL Search(INT(2*k+1,KIND=big),limit,INT(2*threads,KIND=big))
!$OMP END TASK
        END DO
!$OMP TASKWAIT
    END IF
!$OMP END SINGLE
!$OMP END PARALLEL
    CALL Sort(numbers(:count))
    PRINT *, "Carmichael numbers:"
    width = CEILING(LOG10(REAL(limit+1,KIND=KIND(0.0D0))))+1;
    WRITE (buffer, "('(',I0,'I',I0,')')") 80/width,width
    PRINT buffer, numbers(:count)
    PRINT *

CONTAINS

    SUBROUTINE Sort (vector)
        ! Shellsort is quite good enough
        INTEGER(KIND=big), INTENT(INOUT), DIMENSION(:) :: vector
        INTEGER(KIND=big) :: temp
        INTEGER :: limit, i, j, k

        limit = UBOUND(vector,1)
        k = (limit+1)/3
        DO WHILE (k > 0)
            DO i = k+1,limit
                IF (vector(i-k) > vector(i)) THEN
                    temp = vector(i)
                    j = i-k
                    DO WHILE (.TRUE.)
                        vector(j+k) = vector(j)
                        j = j-k
                        IF (j <= 0) EXIT
                        IF (vector(j) <= temp) EXIT
                    END DO
                    vector(j+k) = temp
                END IF
            END DO
            IF (k > 1) k = k+6/k
            k = (k-1)/3
        END DO
    END SUBROUTINE Sort

    RECURSIVE SUBROUTINE Search (from, to, step)
        INTEGER(KIND=big), INTENT(IN) :: from, to, step
        INTEGER(KIND=big) :: n
        LOGICAL, TARGET :: weak, strong

        DO n = from, to, step
            IF (Fermat(2,n) == 1) THEN
!$OMP TASK DEFAULT(NONE), SHARED(weak,n)
                weak = Weak_test(n)
!$OMP END TASK
!$OMP TASK DEFAULT(NONE), SHARED(strong,n)
                strong = Strong_test(n)
!$OMP END TASK
!$OMP TASKWAIT
                IF (weak .NEQV. strong) THEN
!$OMP CRITICAL
                    IF (count >= space) THEN
                        WRITE (*,'("Too many Carmichael numbers!")')
                        STOP 1
                    END IF
                    count = count+1
                    numbers(count) = n
!$OMP END CRITICAL
                END IF
            END IF
         END DO
    END SUBROUTINE Search

    PURE RECURSIVE SUBROUTINE Erastothenes (table)
        INTEGER, DIMENSION(:), INTENT(OUT) :: table
        INTEGER :: i, j, n

        table(1) = 2
        i = 1
        n = 3
        DO WHILE (i < UBOUND(table,1))
            DO j = 1,i
                IF (MOD(n,table(j)) == 0) GOTO 10
            END DO
            table(i+1) = n
            i = i+1
10          n = n+2
        END DO
    END SUBROUTINE Erastothenes

    PURE RECURSIVE FUNCTION Weak_test (trial)
        LOGICAL :: Weak_test
        INTEGER(KIND=big), INTENT(IN) :: trial
        INTEGER :: n

        Weak_test = .TRUE.
        DO n = 1,nprimes
            IF (primes(n) >= trial-1) EXIT
            IF (MOD(trial,INT(primes(n),KIND=big)) == 0) CYCLE
            Weak_test = (Fermat(primes(n),trial) == 1)
            IF (.NOT. Weak_test) EXIT
        END DO
    END FUNCTION Weak_test

    PURE RECURSIVE FUNCTION Strong_test (trial)
        LOGICAL :: Strong_test
        INTEGER(KIND=big), INTENT(IN) :: trial
        INTEGER :: n

        Strong_test = .TRUE.
        DO n = 1,nprimes
            IF (primes(n) >= trial-1) EXIT
            IF (MOD(trial,INT(primes(n),KIND=big)) == 0) CYCLE
            Strong_test = Euler(primes(n),trial)
            IF (.NOT. Strong_test) EXIT
        END DO
    END FUNCTION Strong_test

    PURE RECURSIVE FUNCTION Fermat (base, trial, init)
        INTEGER(KIND=big) :: Fermat
        INTEGER, INTENT(IN) :: base
        INTEGER(KIND=big), INTENT(IN) :: trial
        INTEGER(KIND=big) :: value, mult, exp
        INTEGER(KIND=big), INTENT(IN), OPTIONAL :: init

        IF (PRESENT(init)) THEN
            exp = init
        ELSE
            exp = trial-1
        END IF
        value = 1
        mult = base
        DO WHILE (exp > 0)
            IF (IBITS(exp,0,1) == 1) value = MOD(value*mult,trial)
            exp = ISHFT(exp,-1)
            mult = MOD(mult*mult,trial)
         END DO
         Fermat = value
    END FUNCTION Fermat

    PURE RECURSIVE FUNCTION Euler (base, trial)
        LOGICAL :: Euler
        INTEGER, INTENT(IN) :: base
        INTEGER(KIND=big), INTENT(IN) :: trial
        INTEGER(KIND=big) :: value, exp
        INTEGER :: shift, i

        exp = trial-1
        shift = 0
        DO WHILE (IBITS(exp,0,1) == 0)
            exp = ISHFT(exp,-1)
            shift = shift+1
        END DO
        value = Fermat(base,trial,exp)
        Euler = .TRUE.
        IF (value == 1) RETURN
        DO i = 1,shift
            IF (value == trial-1) RETURN
            value = MOD(value*value,trial)
        END DO
        Euler = .FALSE.
    END FUNCTION Euler

END PROGRAM Main
