pi

real function log2(x)
  implicit none
  real, intent(in) :: x

  log2 = log(x) / log(2.)
end function


integer function find_number_of_primes(limit) result(number_primes)
    implicit none
    INTEGER, INTENT(IN) :: limit
    INTEGER :: i, j
    LOGICAL :: is_prime

    number_primes = 0

    DO i = 3,(limit-1)
        is_prime = .TRUE.
        DO j = 2, (i-1)
            IF (MODULO(i, j) == 0) THEN
                is_prime = .FALSE.
                EXIT
            END IF
        END DO
        IF (is_prime .EQV. .TRUE.) THEN
        number_primes = number_primes + 1
        END IF
    END DO
end function find_number_of_primes

SUBROUTINE get_primes(number_primes, primes)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: number_primes
    INTEGER :: i, found_primes, j
    LOGICAL:: is_prime
    INTEGER, INTENT(OUT) :: primes(number_primes)
    i = 3
    found_primes = 0
    DO
    is_prime = .TRUE.
    DO j = 2, (i-1)
        IF (MODULO(i,j) == 0) THEN
            is_prime = .FALSE.
        END IF
    END DO
    IF (is_prime .EQV. .TRUE.) THEN
        found_primes = found_primes + 1
        primes(found_primes) = i
        IF (found_primes == number_primes) EXIT
    END IF
    i = i + 1
    END DO


end SUBROUTINE get_primes

integer function multiplicity(a, b)
    implicit none
    INTEGER, INTENT(IN) :: a, b

    multiplicity = 0
    DO
        multiplicity = multiplicity + 1
        IF (MOD(b, (a**(multiplicity+1))) /= 0) THEN
            EXIT
        END IF
    END DO
end function multiplicity



PROGRAM pi
    IMPLICIT NONE
    INTEGER :: find_number_of_primes, number_primes, last_digit, digit, N, &
    eps, i, j, multiplicity, test
    REAL :: log2, base, v_max, m, s, v, b, A, pi_sum
    INTEGER, ALLOCATABLE :: primes(:)


    digit = 3
    eps = 2
    base = 10
    N = INT((digit+eps)*log2(base))
    pi_sum = 0
    number_primes = find_number_of_primes(2*N)
    ALLOCATE(primes(number_primes))
    CALL get_primes(number_primes, primes)

    DO i=1,SIZE(primes)
    PRINT*,"************ ", "i =",i,"***************************"
        v_max = INT(log(REAL(2*N))/log(REAL(primes(i))))
        m = primes(i)**v_max
        s = 0
        v = 0
        b = 1
        A = 1
        PRINT*, "v_max =", v_max
        PRINT*, "m =", m
        DO j = 1,(N)
            PRINT*,"------------- ", "j =",j,"-------------------------------"
            b = MOD((j/(primes(i)**multiplicity(digit, j)) * b), m)
            PRINT*, "b =", b
            A = MOD((2*j-1)/(primes(i)**multiplicity(primes(i), (2*j-1)))*A, m)
            PRINT*, "A =", A
            v = v - multiplicity(primes(i),j)+multiplicity(primes(i),(2*j-1))
            PRINT*, "v =", v
            IF (v > 0) THEN
                s = MOD(s*j*b*(A**(-1))*a**(v_max-v), m)
                PRINT*, "v > 0 --> s =", s
            END IF
        END DO
        PRINT*, "------------------------------------------"
        s = MOD(s * base**(digit-1), m)
        PRINT*, "s nach schleife =", s
        pi_sum = pi_sum + MOD((s/m), REAL(1))
        PRINT*, "summe =", pi_sum
    END DO

    PRINT*, pi_sum

    PRINT*, primes

    PRINT*, MOD(6, 4)

END PROGRAM