Bellard's algorithm for calculating Pi

For a university project I wanted to implement Bellard's algorithm for calculating the n-th digit of pi in Fortran. I stumbled across this question on math.stackexchange: https://math.stackexchange.com/questions/1776840/confusion-with-bellards-algorithm-for-pi

With the answer to that question I managed to implement the code, but I'm not getting a result and I can't figure out what I'm doing wrong:

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)
    v_max = INT(log(REAL(2*N))/log(REAL(primes(i))))
    m = primes(i)**v_max
    s = 0
    v = 0
    b = 1
    A = 1
    DO j = 1,(N)
        b = MOD((j/(primes(i)**multiplicity(digit, j)) * b), m)
        A = MOD((2*j-1)/(primes(i)**multiplicity(primes(i), (2*j-1)))*A, m)
        v = v - multiplicity(primes(i),j)+multiplicity(primes(i),(2*j-1))
        IF (v > 0) THEN
            s = MOD(s*j*b*(A**(-1))*a**(v_max-v), m)
        END IF
    END DO
    s = MOD(s * base**(digit-1), m)
    pi_sum = pi_sum + MOD((s/m), REAL(1))
END DO

PRINT*, pi_sum  

END PROGRAM

The custom functions find_number_of_primes, get_primes, log2 and multiplicityare working as intended. First one finds the number of prime numbers in the given interval, second one returns the prime numbers in an array, the third one calculates

log_2(x) = log(x)/log(2)

and the last one calculates the multiplicity (I guess that's what it is called) by checking how often the second number has to be divided by the first number until the rest of the division is no longer zero. Here are the source codes:

The logarithm:

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

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

Finding the number of prime numbers in interval:

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

Getting the actual prime numbers:

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

Pastebin for the whole file: https://pastebin.com/0F4zQYaH

Now in the question that I've linked, at the calculation of b in the inner loop there is a^{v(n, k)} in the denominator, but in the answer to the question the term in the denominator is a^{v(a, k)}. Also, the inner loop in the question runs from 1 to N while the loop in the answer goes from 1 to 2N.

The final result is NaN, here is some console output for digit = 2 and eps = 1:

************ i =           1 ***************************
 v_max =   2.00000000
 m =   9.00000000
 ------------- j =           1 -------------------------------
 b =   0.00000000
 A =   0.00000000
 v =   0.00000000
 ------------- j =           2 -------------------------------
 b =   0.00000000
 A =   0.00000000
 v =   0.00000000
 ------------- j =           3 -------------------------------
 b =   0.00000000
 A =   0.00000000
 v =   0.00000000
 ------------- j =           4 -------------------------------
 b =   0.00000000
 A =   0.00000000
 v =   0.00000000
 ------------- j =           5 -------------------------------
 b =   0.00000000
 A =   0.00000000
 v =   1.00000000
 v > 0 --> s =              NaN
 ------------- j =           6 -------------------------------
 b =   0.00000000
 A =   0.00000000
 v =   1.00000000
 v > 0 --> s =              NaN

and so one, complete output: https://pastebin.com/ucJNg6Vi

Can anybody help me figure out what I'm doing wrong?

Solution

There is a bug in your code:

v = v - multiplicity(primes(i),j)+multiplicity(primes(i),(2*j-1))

should be

v = v - multiplicity(primes(i),j)- multiplicity(primes(i),(2*j-1))

Also,

        s = MOD(s*j*b*(A**(-1))*a**(v_max-v), m)

should be

        s = MOD(s + j*b*(A**(-1))*a**(v_max-v), m)