100 lines
3.1 KiB
Fortran
100 lines
3.1 KiB
Fortran
! Farey Sequence
|
|
!
|
|
! tested with Intel ifx (IFX) 2025.2.1 20250806 on Kubuntu 25.10
|
|
! GNU gfortran (Ubuntu 15.2.0-4ubuntu4) 15.2.0 on Kubuntu 25.10
|
|
! VSI Fortran x86-64 V8.7-001 on OpenVMS V9.2-3
|
|
! U.B., March 2026
|
|
!
|
|
program Farey
|
|
|
|
implicit none
|
|
|
|
! These variables are used in "contained" subtoutines and functions
|
|
integer :: n,a,b,c,d ! Order, and current 2 terms of a Farey Sequence
|
|
|
|
! Local variables here in Main program:
|
|
integer :: ii, count ! Loop index, counter for terms in sequence
|
|
|
|
! compute and show the Farey sequence for orders 1 through 11 (inclusive).
|
|
|
|
do ii=1, 11
|
|
call startNew (ii)
|
|
write (*, '(I2, ": ")', advance='no') n ! Advance 'no' to print all in 1 line
|
|
count = 0
|
|
do
|
|
count = count + 1 ! Count to limit line length
|
|
write (*,'(i0, "/",i0)', advance='no') a, b ! Advance 'no' to print them in 1 line
|
|
if (.not. expectNext()) Exit ! Finish when sequence complete
|
|
write (*, '(", ")', advance='no') ! more is coming - print separator
|
|
if (count .eq. 10) then ! By default, OpenVMS does not like
|
|
count = 0 ! long output lines and produces a run-time
|
|
write (*,'(/,4x)', advance='no') ! error when line length exceeds 132
|
|
endif
|
|
call next() ! Calculate next term
|
|
enddo
|
|
write (*, *) ! Terminate line
|
|
enddo
|
|
|
|
! compute and display the number of fractions in the Farey sequence for order 100 through 1,000 inclusive) by hundreds.
|
|
write (*,*) ! empty line
|
|
do ii=100, 1000, 100
|
|
count = 0
|
|
call startNew (ii)
|
|
write (*, '(I4, ": ")', advance='no') n
|
|
do
|
|
count = count + 1 ! Dont print the terms, just count them
|
|
if (.not. expectNext()) Exit
|
|
call next()
|
|
enddo
|
|
write (*, '(I0)') count
|
|
end do
|
|
contains
|
|
|
|
|
|
! =========================================
|
|
! Start a new Farey Sequence or order newN
|
|
! =========================================
|
|
subroutine startNew(newN)
|
|
integer, intent(in) :: newN
|
|
|
|
n = newN
|
|
a = 0 ! 1st term is a/b=0/1
|
|
b = 1
|
|
c = 1 ! 2nd term is c/d=(1/n)
|
|
d = n
|
|
|
|
end subroutine startNew
|
|
|
|
! =================================================
|
|
! Calculate next term in the current Farey Sequence
|
|
! =================================================
|
|
subroutine next()
|
|
!
|
|
! The algorithm to calculate next term a/b of the sequence is
|
|
! described in Wikipedia, https://en.wikipedia.org/wiki/Farey_sequence
|
|
|
|
integer :: fact, nextC, nextD
|
|
|
|
fact = (n+b) / d
|
|
nextC = fact*c - a
|
|
nextD = fact*d - b
|
|
|
|
a = c
|
|
b = d
|
|
c = nextC
|
|
d = nextD
|
|
|
|
end subroutine next
|
|
|
|
! ==========================================================
|
|
! Say Yes if more terms in current sequence can be expected
|
|
! ==========================================================
|
|
function expectNext() result (YN)
|
|
|
|
! The end is reached when the last term (1/1) is encountered.
|
|
logical :: YN
|
|
|
|
YN = a .ne. 1 .or. b .ne. 1
|
|
end function expectNext
|
|
|
|
end program Farey
|