#Author Matthew Watson. This file is in the public domain.
#!/bin/bash
#Self compiling fortran program in a bash script. Compiles and runs happyab
#Call it sh happy.sh 3 1 for max of 10^3 and output to file on.
#sh happy.sh 8 0 for max of 10^8 and no output
if [ -e "./.happy$1$2" ]
then ./.happy$1$2>gzip -f --fast>happy.txt.gz && exit
else //usr/bin/tail -n +10 $0 | gfortran -O3 -cpp -DPTEN=$1 -DOUT=$2 -o .happy$1$2 -x f95 - && ./.happy$1$2>gzip -f --fast>happy.txt.gz && exit
fi
#if PTEN < 3 || !(defined PTEN)
#define POWTEN 3
#else
#define POWTEN PTEN
#endif
#if OUT==1 || !(defined OUT)
#define DO_OUTPUT .TRUE.
#else
#define DO_OUTPUT .FALSE.
#endif
!Calculate all of the happy numbers up to 10^POWTEN and write the
!result to happy.txt if DO_OUTPUT is set
program happy
implicit none
logical,parameter :: writeoutput = DO_OUTPUT
integer,parameter :: maxnum = 10**POWTEN
!Happy numbers are not very dense, for any value we can possibly
!calculate, there are going to be less than 0.15*n happy numbers.
integer,parameter :: numhap = 15*10**(POWTEN-2)
!We only need to cache a number if we're going to look it up.
!So if our largest number is 99..n..9 or smaller, we only need n*9*9
!cache values.
integer,parameter :: numcache = (POWTEN) * 90
integer :: idx = 1, w = 2, num = 2
integer :: happyI(numhap)
!Cache unhappy numbers rather than happy ones. They're ~10 times as
!common and short cirtuiting the if block speeds things up marginally
!and implies that our number is happy iff the first check is false and
!w<num.
logical :: unHappyArr(numcache) = .FALSE.
unHappyArr(4) = .TRUE. !We need one unhappy number to get started
happyI(1) = 1
do !In this block we calculate the sum of squares
if (num > numcache) then !of digits, check if it's in our unhappy cache.
exit !if so, cache it. If not, check if we've
end if !checked that number already.
!If we have, it must be happy, so num must be
w = sum_sq(w) !happy. Add it to our list.
if(unHappyArr(w)) then
if(num < numcache) then
unHappyArr(num) = .TRUE.
end if
num = num + 1
w = num
else if(w < num) then !Take advantage of the fact that if a number
!is less than num, we must have already
idx = idx + 1 !checked it.
happyI(idx) = num
num = num + 1
w = num
end if
end do
do !Same as the above loop, but we're done caching
if (num > maxnum) then !numbers.
exit !Also, once we are past 162, the possible
end if !values for w are all less than num
!so we no longer need to check w<num
w = sum_sq(w)
if(unHappyArr(w)) then
num = num + 1
w = num
else
idx = idx + 1
happyI(idx) = num
num = num + 1
w = num
end if
end do
if(writeoutput) then
print *,happyI(1:idx)
elseif(idx < 1000) then !Don't print if we have lots, too slow.
print *, happyI(1:idx)
end if
contains
!Calculate the sum of the squares of the digits
!of a number, then return that sum.
function sum_sq (number) result (result)
implicit none
integer, intent (in) :: number
integer :: result,digit,rest,work
result = 0
work = number
do
if (work == 0) then !exit if the remaining number is 0
exit !then we have nothing to do. return
end if
rest = work / 10 !Integer divide by 10 to shift right
digit = work - 10 * rest !Take mod 10 (one digit)
result = result + digit * digit !add square
work = rest !Set number to shifted and repeat
end do
end function sum_sq
end program happy