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

Tags :

Quality of prime seeking methods

Updated July 29, 2018 14:20 PM

Determining the next Twin Prime?

Updated July 30, 2017 03:20 AM