Fortran code snippets

22

Beatty’s theorem

Posted on 9 January 2022 by John D. Cook

Here’s a surprising theorem [1].

(Beatty’s theorem) Let a and b be any pair of irrational numbers greater than 1 with 1/ a + 1/ b = 1. Then the sequences { ⌊ na ⌋ } and { ⌊ nb ⌋ } contain every positive integer without repetition.

Fortran code (Cook used Python) 
module beatty_mod
private
public :: wp, floor_mult, check
integer , parameter :: wp = selected_real_kind(33, 4931) ! kind(1.0d0)
contains
pure function floor_mult(x,n) result(vec)
real(kind=wp), intent(in) :: x
integer      , intent(in) :: n
integer                   :: vec(n)
integer                   :: i
real(kind=wp)             :: xmult
xmult = 0.0_wp
do i=1,n
   xmult = xmult + x
   vec(i) = int(xmult)
end do
end function floor_mult
!
subroutine check(rvec,svec)
integer, intent(in) :: rvec(:),svec(:)
integer             :: i,j,n,r_old
logical, parameter  :: print_svec = .false.
n = size(rvec)
j = 0
do i=2,n
   r_old = rvec(i-1)
   if (rvec(i) /= r_old + 1) then
      j = j + 1
      if (print_svec) print*,svec(j),r_old+1
      if (svec(j) /= r_old + 1) then
         print*,i, j, r_old, rvec(i), svec(j)
         stop "did not work"
      end if
   end if 
end do
end subroutine check
end module beatty_mod
!
program test_beatty
use beatty_mod, only: wp, floor_mult,check
implicit none
real(kind=wp)    , parameter :: r = sqrt(2.0_wp), s = r/(r-1)
character (len=*), parameter :: fmt="(a,*(1x,i0))"
integer          , parameter :: n = 20 ! 10**7
logical          , parameter :: print_rs_vec = .true. ! .false.
integer                      :: rvec(n),svec(n)
rvec = floor_mult(r,n)
svec = floor_mult(s,n)
if (print_rs_vec) then
   print fmt,"{na}",rvec
   print fmt,"{nb}",svec
end if
print "('n, rvec(n) = ',i0,1x,i0)", rvec(n)
call check(rvec,svec)
end program test_beatty

Results:

{na} 1 2 4 5 7 8 9 11 12 14 15 16 18 19 21 22 24 25 26 28
{nb} 3 6 10 13 17 20 23 27 30 34 37 40 44 47 51 54 58 61 64 68
n, rvec(n) = 28

The numbers not in {na} are in {nb} and {na} and {nb} are disjoint.

Changing two lines in the main program to

integer          , parameter :: n = 10**8 ! n = 20
logical          , parameter :: print_rs_vec = .false.

the program, which uses quadruple precision, runs in a few seconds and verifies Beatty’s theorem up to that N.

1 Like