Bellard's algorithm for calculating Pi

by brettljausn   Last Updated October 16, 2018 13:20 PM

For a university project I wanted to implement Bellard's algorithm for calculating the $n$-th digit $\pi$ in FORTRAN. I stumbled across this question on math.stackexchange: Confusion with Bellard's 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) = \frac{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$.

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

EDIT: the final result is NaN, heres some console output for digit = 2 and eps = 1: https://pastebin.com/ucJNg6Vi



Related Questions


Quality of prime seeking methods

Updated July 29, 2018 14:20 PM

Determining the next Twin Prime?

Updated July 30, 2017 03:20 AM

Lower bounds for integer factorization

Updated July 11, 2017 12:20 PM