Solve p15 in fortran

docs/src/fortran/lib/constants.rst

@@ -32,14 +32,10 @@ View source code :source:`fortran/include/constants.f90`
     :type: integer
 .. f:variable:: ERROR_ANSWER_TYPE_MISMATCH
     :type: integer
-.. f:variable:: ERROR_PRIME_ALLOCATE_FAILED
-    :type: integer
-.. f:variable:: ERROR_UTILS_ALLOCATE_FAILED
+.. f:variable:: ERROR_ALLOCATE_FAILED
     :type: integer
 .. f:variable:: ERROR_FILE_READ_FAILED
     :type: integer
-.. f:variable:: ERROR_PROB_ALLOCATE_FAILED
-    :type: integer
 
     Denotes the exit codes of different failure modes, counting up from 1
 

docs/src/fortran/lib/math.rst

@@ -11,6 +11,24 @@ View source code :source:`fortran/include/math.f90`
     :returns answer: n!
     :rtype answer: integer(i18t)
 
+.. f:function:: n_choose_r(n, r)
+
+    :p integer n:
+    :p integer r:
+    :returns answer: nCr
+    :rtype answer: integer(i18t)
+
+.. f:function:: n_choose_r_slow(n, r)
+
+    This is a helper function for numbers that would otherwise overflow using the naive calculation of
+    ``factorial(n) / factorial(r) / factorial(n - r)``. It's a slower path, but a much more reliable one
+    because it actually factorizes the numbers and cancels out as many as possible before multiplying.
+
+    :p integer n:
+    :p integer r:
+    :returns answer: nCr
+    :rtype answer: integer(i18t)
+
 .. literalinclude:: ../../../../fortran/src/include/math.f90
    :language: Fortran
    :linenos:

docs/src/fortran/p0010.rst

@@ -3,6 +3,14 @@ Fortran Implementation of Problem 10
 
 View source code :source:`fortran/src/p0010.f90`
 
+Includes
+--------
+
+- `primes.f90 <./lib/primes.html>`_
+
+Problem Solution
+----------------
+
 .. f:module:: Problem0010
 
 .. f:function:: integer Problem0010/p0010()

docs/src/fortran/p0015.rst

@@ -0,0 +1,22 @@
+Fortran Implementation of Problem 15
+====================================
+
+View source code :source:`fortran/src/p0015.f90`
+
+Includes
+--------
+
+- `math.f90 <./lib/math.html>`_
+
+Problem Solution
+----------------
+
+.. f:module:: Problem0015
+
+.. f:function:: integer Problem0015/p0015()
+
+.. literalinclude:: ../../../fortran/src/p0015.f90
+   :language: Fortran
+   :linenos:
+
+.. tags:: combinatorics

fortran/src/include/constants.f90

@@ -16,13 +16,12 @@ module constants
     integer, parameter :: ERROR_ANSWER_MISMATCH = 1
     integer, parameter :: ERROR_ANSWER_TIMEOUT = 2
     integer, parameter :: ERROR_ANSWER_TYPE_MISMATCH = 3
-    integer, parameter :: ERROR_PRIME_ALLOCATE_FAILED = 4
-    integer, parameter :: ERROR_UTILS_ALLOCATE_FAILED = 5
-    integer, parameter :: ERROR_FILE_READ_FAILED = 6
-    integer, parameter :: ERROR_PROB_ALLOCATE_FAILED = 7
+    integer, parameter :: ERROR_ALLOCATE_FAILED = 4
+    integer, parameter :: ERROR_FILE_READ_FAILED = 5
 
     ! file/string sizes
     integer, parameter :: DATA_MAX_NAME_SIZE = 32
     integer, parameter :: ANSWERT_STR_SIZE = 16
+    integer, parameter :: MAX_FACTORIAL_64 = 20
 end module constants
 

fortran/src/include/math.f90

@@ -10,5 +10,78 @@ pure integer(i18t) function factorial(n) result(answer)
             answer = answer * i
         end do
     end function
+
+    pure integer(i18t) function n_choose_r(n, r) result(answer)
+        integer, intent(in) :: n, r
+        if (n <= MAX_FACTORIAL_64) then
+            answer = factorial(n) / factorial(r) / factorial(n-r)  ! fast path if small enough
+        else
+            answer = n_choose_r_slow(n, r)  ! slow path for larger numbers
+        end if
+    end function
+
+    pure integer(i18t) function n_choose_r_slow(n, r) result(answer)
+        integer(i2t), allocatable :: factors(:)
+        integer(i18t) :: tmp
+        integer :: i, j
+        allocate(factors(n))
+        if (.not. allocated(factors)) then
+            print *, "n_choose_r allocation failed. Exiting."
+            stop ERROR_ALLOCATE_FAILED
+        end if
+        factors = 0
+        ! collect factors of final number
+        do i = 2, n
+            factors(i) = 1
+        end do
+        ! negative factor values indicate need to divide
+        do i = 2, r
+            factors(i) = factors(i) - 1_i2t
+        end do
+        do i = 2, n - r
+            factors(i) = factors(i) - 1_i2t
+        end do
+        do i = n, 2, -1  ! this loop reduces to prime factors only
+            do j = 2, i - 1
+                if (mod(i, j) == 0) then
+                    factors(j) += factors(i)
+                    factors(i / j) += factors(i)
+                    factors(i) = 0;
+                    exit
+                end if
+            end do
+        end do
+        i = 2
+        j = 2
+        answer = 1
+        do while (i <= n)
+            do while (factors(i) > 0)
+                tmp = answer
+                answer = answer * i
+                do while (answer < tmp .and. j <= n)
+                    do while (factors(j) < 0)
+                        tmp = tmp / j
+                        factors(j) = factors(j) + 1_i2t
+                    end do
+                    j = j + 1_i2t
+                    answer = tmp * i
+                end do
+                if (answer < tmp) then
+                    answer = -1  ! this indicates an overflow
+                    return
+                end if
+                factors(i) = factors(i) - 1_i2t
+            end do
+            i = i + 1_i2t
+        end do
+        do while (j <= n)
+            do while (factors(j) < 0)
+                answer = answer / j
+                factors(j) = factors(j) + 1_i2t
+            end do
+            j = j + 1_i2t
+        end do
+        deallocate(factors)
+    end function
 end module math
 

fortran/src/include/primes.f90

@@ -56,7 +56,7 @@ subroutine expand_sieve(potential_n)
         allocate(prime_sieve(new_size))
         if (.not. allocated(prime_sieve)) then
             print *, "Memory allocation failed for prime sieve. Exiting."
-            stop ERROR_PRIME_ALLOCATE_FAILED
+            stop ERROR_ALLOCATE_FAILED
         end if
         prime_sieve = -1
         call clear_prime_bit(0_i18t)

fortran/src/p0014.f90

@@ -29,7 +29,7 @@ integer(i18t) function p0014() result(answer)
         allocate(collatz_len_cache(collatz_cache_size))
         if (.not. allocated(collatz_len_cache)) then
             print *, "Allocation of collatz length cache failed. Stopping run."
-            stop ERROR_PROB_ALLOCATE_FAILED
+            stop ERROR_ALLOCATE_FAILED
         end if
         collatz_len_cache = 0
         collatz_len_cache(1) = 1

fortran/src/p0015.f90

@@ -0,0 +1,24 @@
+! Project Euler Problem 15
+!
+! Turns out this is easy, if you think sideways a bit
+!
+! You can only go down or right. If we say right=1, then you can only have 20 1s, since otherwise you go off the grid.
+! You also can't have fewer than 20 1s, since then you go off the grid the other way. This means you can look at it as a
+! bit string, and the number of 40-bit strings with 20 1s is 40c20.
+!
+! Problem:
+!
+! Starting in the top left corner of a 2×2 grid, and only being able to move to the right and down, there are exactly 6
+! routes to the bottom right corner.
+!
+! How many such routes are there through a 20×20 grid?
+
+module Problem0015
+    use constants
+    use math
+    implicit none
+contains
+    integer(i18t) function p0015() result(answer)
+        answer = n_choose_r(40, 20)
+    end function p0015
+end module Problem0015