How can we implement continued fraction statements in fortran?
I’m not sure, about your question.
- Do you want, the value of “f” from the two lists of integers a_i and b_i?
=> it is relatively easy.
- or do you want the decomposition of “f” into the two lists of integers?
=> You need an algorithm for that (see wikipedia)
By the way, for continued fraction, normally, the a_i are equal to one.
As stated by @gardhor, it is relatively easy if you use the Euler-Wallis recurrence relations. Then you just need to write a recursive function:
module continued_fractions implicit none integer, parameter :: pi_seq(9) = [3, 7, 15, 1, 292, 1, 1, 1, 2] !! Simple continued fraction expansion of Pi taken from https://oeis.org/A001203 contains !> Returns the *simple* continued fraction of the sequence b recursive function cf(b) result(pq) integer, intent(in) :: b(:) real :: pq(2) real :: pq_n1(2), pq_n2(2) integer :: n n = size(b) if (n > 2) then pq_n2 = cf(b(1:n-2)) pq_n1 = cf(b(1:n-1)) pq = b(n)*pq_n1 + pq_n2 else if (n == 2) then pq_n1 = cf(b(1:1)) pq(1) = b(2)*pq_n1(1) + 1.0 pq(2) = b(2) else if (n == 1) then pq(1) = b(1) pq(2) = 1.0 end if end function end module program main use continued_fractions, only: cf, pi_seq implicit none real :: f, pq(2) real, parameter :: pi = 4.0*atan(1.0) pq = cf(pi_seq) f = pq(1)/pq(2) print *, "continued fraction:", f print *, "true value:", pi end program
I did not implement the scaling part to prevent overflow. The reason you need it however, is the values of p_n and q_n quickly become large integers as you add more terms. In the example above with just 9 terms of the continued fractions, the values in
To reach the full 20 fractions, it seems like you will need the rescaling.
As a cool side note, the fourth convergent of \pi, [3;7,15,1] = 355/113, is also known as Milü after a Chinese astronomer born 429 AD. It is within 0.000009% of the value of \pi.