euler.f90

! This file is part of EulerFTN.
! Copyright (C) 2013-2014 Adam Hirst <adam@aphirst.karoo.co.uk>
!
! EulerFTN is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! EulerFTN is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with EulerFTN. If not, see <http://www.gnu.org/licenses/>.

module Euler
  use Functions
  implicit none

contains

  subroutine Problem_01
    ! If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9.
    ! The sum of these multiples is 23.
    !
    ! Find the sum of all the multiples of 3 or 5 below 1000.
    integer(long)              :: n, upper, i, the_sum
    integer(long), allocatable :: factors(:)
    logical(1),    allocatable :: sieve(:)
    real                       :: time(2)

    print *, 'Problem #1: Find the sum of all the multiples of <FACTORS> below <INPUT>.'
    print *, ''
    print *, 'How many factors?'
    read *, n
    allocate(factors(n))
    print *, 'Please input the factors:'
    read *, factors
    print *, 'Please input the exclusive upper bound:'
    read *, upper
    call cpu_time(time(1))
    ! we use a simple sieving algorithm to mark the multiples, then sum all the .true. indices
    sieve = [( .false._1, i = 1,upper-1 )]
    do i = 1, size(factors)
      sieve(factors(i)::factors(i)) = .true._1
    end do
    the_sum = sum(pack( [( i, i = 1,upper )] , sieve ))
    call cpu_time(time(2))
    print '(a,i0)', 'The sum of the inputted factors is ', the_sum
    print '(a,f0.3,a)', 'Took ',time(2)-time(1),' seconds'
  end subroutine Problem_01

  subroutine Problem_02
    ! Each new term in the Fibonacci sequence is generated by adding the previous two terms.
    ! By starting with 1 and 1, the first 10 terms will be:
    ! 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, ...
    !
    ! By considering the terms in the Fibonacci sequence whose values do not exceed four million,
    ! find the sum of the even-valued terms.
    integer(long) :: fibmax, nmax, n
    real          :: time(2)

    print *, 'Problem #2: Find the sum of the even-valued terms of the Fibonacci sequence,'
    print *, 'using only the terms which do not exceed <INPUT>.'
    print *, ''
    print *, 'Please input the INCLUSIVE upper bound for the sum:'
    read *, fibmax
    call cpu_time(time(1))
    ! determine index 'nmax' of first F_n STRICTLY BELOW fibmax
    nmax = Inv_fib(fibmax)
    ! by inspection, only the (3n)th Fibonacci numbers are even
    print *, sum( [( Fib(n), n=3,nmax,3 )] )
    call cpu_time(time(2))
    print '(a,f0.3,a)', 'Took ',time(2)-time(1),' seconds'
  end subroutine Problem_02

  subroutine Problem_03
    ! The prime factors of 13195 are 5, 7, 13 and 29.
    ! What is the largest prime factor of the number 600851475143 ?
    integer(long)              :: max_prime, to_factorise, i ! up to 9223372036854775807 for me
    integer(long), allocatable :: primes(:), wheel(:), factors(:)
    real(dp)                   :: time(2)

    print *, 'Problem #3: What is the largest prime factor of the number <INPUT>?'
    print *, ''
    print *, 'Please input the number to factorise:'
    read *, to_factorise
    if (to_factorise <= 0) stop 'Invalid: please input a positive integer.'
    print *, 'What is the largest (prime) number to use for the factorisation wheel?'
    print *, 'Note: any higher than 13 generates too large a wheel for most machines.'
    read *, max_prime
    call cpu_time(time(1))
    ! use simple prime sieve to generate first few primes
    primes = Primes_sieve(max_prime)
    ! generate a 'wheel' using those primes
    wheel = Primes_wheel(primes)
    ! use that wheel to determine all the prime factors of 'to_factorise'
    factors = Prime_factors(to_factorise, primes, wheel)
    call cpu_time(time(2))
    print *, 'The complete list of prime factors is:'
    ! this is a bit awkward, but it does print it out nicely
    write(*,'(a)',advance='no') '[ '
    do i = 1, size(factors)-1
      write(*,'(i0)',advance='no') factors(i)
      write(*,'(a)',advance='no') ', '
    end do
    print '(i0,a)', factors(size(factors)),' ]'
    print '(a,i0)', 'and the largest prime factor is ',factors(size(factors)) ! generated in ascending order
    print '(a,f0.3,a)', 'Took: ',time(2)-time(1),' seconds'
  end subroutine Problem_03

  subroutine Problem_04
    ! A palindromic number reads the same both ways.
    ! The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 * 99.
    ! Find the largest palindrome made from the product of two 3-digit numbers.
    integer(long) :: n, half, div, min_div, max_div, current_pal
    real(dp)      :: time(2)

    print *, 'Problem #4: What is the largest palindromic number made from the product of two <INPUT>-digit numbers?'
    print *, ''
    print *, 'Please input the number of digits for both these factors:'
    read *, n
    ! make sure that the resulting palindrome will fit in the integer KIND being used
    if ( real(2*n) >= log10(real( huge(n) )) ) then
      stop 'Too many digits, product irrepresentable as current integer'
    end if
    ! start the computation proper
    call cpu_time(time(1))
    max_div = (10**n) - 1
    outer: do half = (10**n)-3, 10**(n-1), -1
      current_pal = Palindrome(half)
      ! we only have to check against divisors that could give n-digit co-divisors
      min_div = current_pal / max_div
      ! even-digit (2n) palindromes are all divisible by 11
      ! therefore we only have to check divisors also divisible by 11
      inner: do div = max_div - modulo(max_div,11_long), min_div, -11
        !print *, div
        if ( modulo(current_pal,div) == 0 ) then
          ! we're done
          exit outer ! so we break out of the outer loop as well
        end if
      end do inner
    end do outer
    call cpu_time(time(2))
    !print *, 'max:',max_div - modulo(max_div,11_long),' min:', min_div, div
    print *, 'The largest such palindrome is:'
    print '(i0,a,i0,a,i0)', current_pal,' = ',div,' * ',current_pal/div
    print '(a, f0.3, a)', 'Took: ',time(2)-time(1),' seconds'
  end subroutine Problem_04

  subroutine Problem_05
    ! 2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.
    ! What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?
    integer(long)              :: i, limit
    integer(long), allocatable :: divisors(:)
    real(dp)                   :: time(2)

    print *, 'Problem #5: What is the smallest positive number that is evently divisible by all the numbers from 1 to <INPUT>?'
    print *, ''
    print *, 'Upper bound (inclusive)?'
    read *, limit
    divisors = [( i, i=1,limit )]
    call cpu_time(time(1))
    print '(i0)', Lowest_common_multiple( divisors ) ! 232792560
    call cpu_time(time(2))
    print '(a,f0.3,a)', 'Took: ',time(2)-time(1),' seconds'
  end subroutine Problem_05

  subroutine Problem_06
    ! The difference between the sum of the squares and the square of the sums of the first ten natural numbers is 2640.
    ! Find the difference between the sum of the squares and the square of the sums for the first hundred natural numbers.
    integer(long) :: n

    print *, 'Problem #6: What is the difference between the sum of the squares and the square of the sums for the first <INPUT>&
             &natural numbers?'
    print *, ' '
    print *, 'Upper bound (inclusive)?'
    read *, n
    ! the order here ensures no division rounding, and minimises the chance of integer overflow
    print '(i0)', ( ( (n*(n+1))/2 ) * ( ((n-1)*(3*n + 2))/2 ) ) / 3 ! fails somewhere over N = 50000
  end subroutine Problem_06

  subroutine Problem_07
    ! By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
    ! What is the 10001st prime number?
    integer(long), allocatable :: primes(:)
    integer(long)              :: n
    real                       :: time(2)

    print *, 'Problem #7: What is the <INPUT>th prime number?'
    print *, ' '
    print *, 'Which prime number should we find?'
    read *, n
    call cpu_time(time(1))
    primes = Primes_sieve(Largest_pn(n))
    call cpu_time(time(2))
    print '(a,i0,a,i0)', 'The ',n,'th prime number is ',primes(n)
    print '(a,f0.3,a)', 'Took: ',time(2)-time(1),' seconds'
  end subroutine Problem_07

 subroutine Problem_08
    ! Find the greatest product of five consecutive digits in the 1000-digit number.
    integer, allocatable :: numbers(:)
    integer              :: lines = 20, linewidth = 50, i, max_product = 1
    real                 :: time(2)

    print *, 'Problem #8: What is the greatest product of 5 consecutive numbers in the 1000-digit number in file <INPUT>?'
    ! load the file into an array of integers
    print *, 'Default: number stored in "problem8.txt" with 50 numbers per line.'
    open(5, file='problem_08.txt', status='old')
    ! 50 integers per line, 20 lines, no spaces
    call cpu_time(time(1))
    allocate(numbers(lines*linewidth))
    do i = 1,lines
      read(5,'(*(i1))') numbers(linewidth*i-linewidth+1:linewidth*i)
    end do
    close(5)
    do i = 1, size(numbers)-4
      if (any( numbers(i:i+4) == 0 )) cycle
      max_product = max(max_product, product(numbers(i:i+4)))
    end do
    call cpu_time(time(2))
    print '(a,i0)', 'The greatest product of 5 consecutive numbers is ',max_product
    print '(a,f0.3,a)', 'Took: ',time(2)-time(1),' seconds'
  end subroutine Problem_08

  subroutine Problem_09
    ! There exists only one Pythagorean triplet for which `a+b+c=1000` Find the product `abc`.
    integer(long) :: total, a, b, c
    real          :: time(2)

    print *, 'Problem #9: What is the product `abc`, where `a^2 + b^2 = c^2` (Pythagorean triplet) and `a + b + c = <INPUT>?'
    print *, ''
    print *, 'Please input the desired sum for `a+b+c`:'
    read *, total
    call cpu_time(time(1))
    ! upper bound for `a` is `total/3` since `a < b < 3`, and it's the sum of *3* numbers...
    outer: do a = 1, total/3
      ! similarly for `b`, for the smallest `a` the sum `total` has to come from the remaining *2* numbers...
      do b = a, total/2
        ! clearly we know what `c` should be in each case
        c = total - a - b
        ! if the triplet is also Pythagorean, then we're done
        if (a**2 + b**2 == c**2) then
          exit outer
        end if
      end do
    end do outer
    call cpu_time(time(2))
    print *, a, b, c, a+b+c, a*b*c
    print '(a,i0)', 'The product of the triplet is ', a*b*c
    print '(a,f0.3,a)', 'Took: ',time(2)-time(1),' seconds'
  end subroutine Problem_09

  subroutine Problem_10
    ! The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.
    ! Find the sum of all the primes below two million.
    integer(long) :: n, the_sum
    real          :: time(2)

    print *, 'Problem #10: What is the sum of all the primes below <INPUT>?'
    print *, ''
    print *, 'Upper bound (exclusive)?'
    read *, n
    call cpu_time(time(1))
    the_sum = sum(Primes_sieve(n-1))
    call cpu_time(time(2))
    print '(a,i0)', 'The sum is ',the_sum
    print '(a,f0.3,a)', 'Time: ',time(2)-time(1),' seconds'
  end subroutine Problem_10

  subroutine Problem_11
    ! In the 20x20 grid below, four numbers along a diagonal line have been marked in red.
    ! The product of these numbers is 26 * 63 * 78 * 14 = 1788696.
    ! What is the greatest product of four adjacent numbers in the same direction (including diagonally) in the 20x20 grid?
    integer              :: width = 20, height = 20, i, j, k, temp(4), max_product = 0
    integer, allocatable :: grid(:,:)
    real                 :: time(2)

    print *, 'What is the greatest product of four adjacent numbers (including diagonally) in the number grid in file <INPUT>?'
    print *, 'Defaulting to 20x20 grid at "problem_11.txt".'
    allocate(grid(height, width))
    open(5, file='problem_11.txt', status='old')
    ! we'll be naughty and load each row from the file into a column of the array, because it's quicker
    ! in effect, we're storing the transpose of the grid
    do i = 1, width
      read(5,'(20(i3))') grid(:,i)
    end do
    close(5)
    call cpu_time(time(1))
    ! read horizontally
    do i = 1, height
      do j = 1, width-3
        if (any( grid(i,j:j+3) == 0 )) cycle
        max_product = max( max_product, product(grid(i,j:j+1)) )
      end do
    end do
    ! read vertically
    do i = 1, width
      do j = 1, height-3
        if (any( grid(j:j+3,i) == 0 )) cycle
        max_product = max( max_product, product(grid(j:j+3,i)) )
      end do
    end do
    ! read down-right diagonal
    do i = 1, width-3
      do j = 1, height-3
        temp = [( grid(j+k,i+k), k=0,3 )]
        if (any( temp == 0 )) cycle
        max_product = max( max_product, product(temp) )
      end do
    end do
    ! read up-right diagonal
    do i = 1, width-3
      do j = 4, height
        temp = [( grid(j-k,i+k), k=0,3 )]
        if (any( temp == 0 )) cycle
        max_product = max( max_product, product(temp) )
      end do
    end do
    call cpu_time(time(2))
    print *, max_product
    print '(a,f0.3,a)', 'Time: ',time(2)-time(1),' seconds'
  end subroutine Problem_11

  subroutine Problem_13
  end subroutine Problem_13

  subroutine Problem_14
    ! Which starting number, under one million, produces the longest Collatz chain?
    integer(long)              :: i, n, max_n, max_length
    integer(long), allocatable :: lengths(:)
    real(dp)                   :: time(2)

    print *, 'Problem #14: Find the starting number producing the longest Collatz chain, with the start below <INPUT>.'
    print *, ''
    print *, 'Upper bound (exclusive)?'
    read *, n
    lengths = [( 0_long, i=1,n-1 )]
    call cpu_time(time(1))
    do i = 1, n
      if (lengths(i) == 0) lengths(i) = Collatz(i,lengths)
    end do
    max_length = maxval(lengths)
    max_n = sum(maxloc(lengths))
    call cpu_time(time(2))
    print '(a,i0,a,i0,a)', 'Longest sequence starts at ',max_n,', and is ',max_length,' long.'
    print '(a,f0.3,a)', 'Took: ',time(2) - time(1),' seconds'
  end subroutine Problem_14

end module Euler