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 `multiplicity`

are 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

- ServerfaultXchanger
- SuperuserXchanger
- UbuntuXchanger
- WebappsXchanger
- WebmastersXchanger
- ProgrammersXchanger
- DbaXchanger
- DrupalXchanger
- WordpressXchanger
- MagentoXchanger
- JoomlaXchanger
- AndroidXchanger
- AppleXchanger
- GameXchanger
- GamingXchanger
- BlenderXchanger
- UxXchanger
- CookingXchanger
- PhotoXchanger
- StatsXchanger
- MathXchanger
- DiyXchanger
- GisXchanger
- TexXchanger
- MetaXchanger
- ElectronicsXchanger
- StackoverflowXchanger
- BitcoinXchanger
- EthereumXcanger