Want more features on Pastebin? Sign Up, it's FREE!
Guest

Self compiling fortran happy number program in a bash script

By: a guest on May 10th, 2012  |  syntax: None  |  size: 4.01 KB  |  views: 26  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #Author Matthew Watson. This file is in the public domain.
  2. #!/bin/bash
  3. #Self compiling fortran program in a bash script. Compiles and runs happyab
  4. #Call it sh happy.sh 3 1 for max of 10^3 and output to file on.
  5. #sh happy.sh 8 0 for max of 10^8 and no output
  6. if [ -e "./.happy$1$2" ]
  7.   then ./.happy$1$2>gzip -f --fast>happy.txt.gz && exit
  8.   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
  9. fi
  10. #if PTEN < 3 || !(defined PTEN)
  11. #define POWTEN 3
  12. #else
  13. #define POWTEN PTEN
  14. #endif
  15. #if OUT==1 || !(defined OUT)
  16. #define DO_OUTPUT .TRUE.
  17. #else
  18. #define DO_OUTPUT .FALSE.
  19. #endif
  20. !Calculate all of the happy numbers up to 10^POWTEN and write the
  21. !result to happy.txt if DO_OUTPUT is set
  22. program happy
  23.   implicit none
  24.   logical,parameter :: writeoutput = DO_OUTPUT
  25.   integer,parameter :: maxnum = 10**POWTEN
  26.  
  27.   !Happy numbers are not very dense, for any value we can possibly
  28.   !calculate, there are going to be less than 0.15*n happy numbers.
  29.   integer,parameter :: numhap = 15*10**(POWTEN-2)
  30.  
  31.   !We only need to cache a number if we're going to look it up.
  32.   !So if our largest number is 99..n..9 or smaller, we only need n*9*9
  33.   !cache values.
  34.   integer,parameter :: numcache = (POWTEN) * 90
  35.  
  36.   integer :: idx = 1, w = 2, num = 2
  37.   integer :: happyI(numhap)
  38.  
  39.   !Cache unhappy numbers rather than happy ones. They're ~10 times as
  40.   !common and short cirtuiting the if block speeds things up marginally
  41.   !and implies that our number is happy iff the first check is false and
  42.   !w<num.
  43.  
  44.   logical :: unHappyArr(numcache) = .FALSE.
  45.  
  46.   unHappyArr(4) = .TRUE.        !We need one unhappy number to get started
  47.   happyI(1) = 1
  48.  
  49.   do                            !In this block we calculate the sum of squares
  50.     if (num > numcache) then    !of digits, check if it's in our unhappy cache.
  51.       exit                      !if so, cache it. If not, check if we've
  52.     end if                      !checked that number already.
  53.                                 !If we have, it must be happy, so num must be
  54.     w = sum_sq(w)               !happy. Add it to our list.
  55.  
  56.     if(unHappyArr(w)) then
  57.  
  58.       if(num < numcache) then
  59.         unHappyArr(num) = .TRUE.
  60.       end if
  61.       num = num + 1
  62.       w = num
  63.  
  64.     else if(w < num) then       !Take advantage of the fact that if a number
  65.                                 !is less than num, we must have already
  66.       idx = idx + 1             !checked it.
  67.       happyI(idx) = num
  68.       num = num + 1
  69.       w = num
  70.     end if
  71.  
  72.   end do
  73.   do                            !Same as the above loop, but we're done caching
  74.     if (num > maxnum) then      !numbers.
  75.       exit                      !Also, once we are past 162, the possible
  76.     end if                      !values for w are all less than num
  77.                                 !so we no longer need to check w<num
  78.     w = sum_sq(w)
  79.  
  80.     if(unHappyArr(w)) then
  81.       num = num + 1
  82.       w = num
  83.     else
  84.       idx = idx + 1
  85.       happyI(idx) = num
  86.       num = num + 1
  87.       w = num
  88.     end if
  89.  
  90.   end do
  91.   if(writeoutput) then
  92.     print *,happyI(1:idx)
  93.   elseif(idx < 1000) then         !Don't print if we have lots, too slow.
  94.     print *, happyI(1:idx)
  95.   end if
  96.  
  97. contains
  98.  
  99.   !Calculate the sum of the squares of the digits
  100.   !of a number, then return that sum.
  101.  
  102.  
  103.   function sum_sq (number) result (result)
  104.  
  105.     implicit none
  106.     integer, intent (in) :: number
  107.     integer :: result,digit,rest,work
  108.     result = 0
  109.     work = number
  110.     do
  111.       if (work == 0) then             !exit if the remaining number is 0
  112.         exit                          !then we have nothing to do. return
  113.       end if
  114.       rest = work / 10                !Integer divide by 10 to shift right
  115.       digit = work - 10 * rest        !Take mod 10 (one digit)
  116.       result = result + digit * digit !add square
  117.       work = rest                     !Set number to shifted and repeat
  118.     end do
  119.  
  120.   end function sum_sq
  121.  
  122. end program happy
clone this paste RAW Paste Data