Hi all,
How can we implement continued fraction statements in fortran?
For example:
Hi,
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 pq
are:
833719.000 265381.000
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.