Для университетского проекта я хотел реализовать алгоритм Белларда для вычисления n-й цифры числа пи в Фортране.Я наткнулся на этот вопрос на math.stackexchange: https://math.stackexchange.com/questions/1776840/confusion-with-bellards-algorithm-for-pi
С ответом на этот вопрос мне удалось реализовать код, но я не получил результат, и я не могу понять, что я 'Я делаю неправильно:
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
Пользовательские функции find_number_of_primes, get_primes, log2 и multiplicity работают как задумано.Первый находит число простых чисел в заданном интервале, второй возвращает простые числа в массиве, третий вычисляет
log_2(x) = log(x)/log(2)
, а последний вычисляет кратность (я думаю, что это и есть)вызывается), проверяя, как часто второе число должно делиться на первое число, пока остальная часть деления больше не станет нулевой.Вот исходные коды:
Логарифм:
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
Вставьте для всего файла: https://pastebin.com/0F4zQYaH
Теперь в вопросе, который я связал, при расчете b во внутреннем цикле есть a ^ {v (n, k)} в знаменателе, но в ответе на вопрос термин в знаменателе равен a ^ {v (a, k)} .Кроме того, внутренний цикл в вопросе проходит от 1 до N , тогда как цикл в ответе изменяется от 1 до 2N .
Окончательный результат - NaN, вот несколько выводов консоли для digit = 2 и 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
и так далее, полный вывод: https://pastebin.com/ucJNg6Vi
Может кто-нибудь помочь мне понять, что я делаю не так?