|
| 1 | + |
| 2 | +!Copyright (C) 2021 Richard Weed. |
| 3 | +!All rights reserved. |
| 4 | + |
| 5 | +!Redistribution and use in source and binary forms, with or without |
| 6 | +!modification, are permitted provided that the following conditions are met: |
| 7 | + |
| 8 | +!1. Redistributions of source code, in whole or in part, must retain the |
| 9 | +!above copyright notice, this list of conditions and the following |
| 10 | +!disclaimer. |
| 11 | + |
| 12 | +!2. Redistributions in binary form, in whole or in part, must reproduce the |
| 13 | +!above copyright notice, this list of conditions and the following disclaimer |
| 14 | +!in the documentation and/or other materials provided with the distribution. |
| 15 | + |
| 16 | +!3. The names of the contributors may not be used to endorse or promote from |
| 17 | +!products derived from this software without specific prior written |
| 18 | +!permission. |
| 19 | + |
| 20 | +!4. Redistributions of this software, in whole or in part, in any form, |
| 21 | +!must be freely available and licensed under this original License. The |
| 22 | +!U.S. Government may add additional restrictions to their modified and |
| 23 | +!redistributed software as required by Law. However, these restrictions |
| 24 | +!do not apply to the original software distribution. |
| 25 | + |
| 26 | +!5. Redistribution of this source code, including any modifications, may |
| 27 | +!not be intentionally obfuscated. |
| 28 | + |
| 29 | +!6. Other code may make use of this software, in whole or in part, without |
| 30 | +!restriction, provided that it does not apply any restriction to this |
| 31 | +!software other than outlined above. |
| 32 | + |
| 33 | +!THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS |
| 34 | +!IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, |
| 35 | +!THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR |
| 36 | +!PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND |
| 37 | +!CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| 38 | +!EXEMPLARARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| 39 | +!PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; |
| 40 | +!OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, |
| 41 | +!WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR |
| 42 | +!OTHERWISE), ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF |
| 43 | +!ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| 44 | + |
| 45 | +Module binomial_coefficients |
| 46 | + |
| 47 | +! Utilities to compute/return binomial coefficients. Utilities will |
| 48 | +! compute a table of n binomial coefficients, an array containing the |
| 49 | +! n+1 coefficients for a degree n polynomial, and single binomial |
| 50 | +! coefficients. Only double precision supported for now. |
| 51 | + |
| 52 | +USE ISO_FORTRAN_ENV, ONLY: DP=>REAL64 |
| 53 | + |
| 54 | + Integer, PRIVATE :: mm |
| 55 | + Real(DP), PRIVATE, SAVE :: binomialCoefsTo20(0:20, 0:20) = TRANSPOSE( & |
| 56 | + RESHAPE(REAL([ & |
| 57 | + !n=0 |
| 58 | + 1, (0, mm=1,20), & |
| 59 | + !n=1 |
| 60 | + 1, 1, (0, mm=1,19), & |
| 61 | + !n=2 |
| 62 | + 1, 2, 1, (0, mm=1,18), & |
| 63 | + !n=3 |
| 64 | + 1, 3, 3, 1, (0, mm=1,17), & |
| 65 | +! n=4 |
| 66 | + 1, 4, 6, 4, 1, (0, mm=1,16), & |
| 67 | +! n=5 |
| 68 | + 1, 5, 10, 10, 5, 1, (0, mm=1,15), & |
| 69 | +! n=6 |
| 70 | + 1, 6, 15, 20, 15, 6, 1, (0, mm=1,14), & |
| 71 | +! n=7 |
| 72 | + 1, 7, 21, 35, 35, 21, 7, 1, (0, mm=1,13), & |
| 73 | +! n=8 |
| 74 | + 1, 8, 28, 56, 70, 56, 28, 8, 1, (0, mm=1,12), & |
| 75 | +! n=9 |
| 76 | + 1, 9, 36, 84, 126, 126, 84, 36, 9, 1, (0, mm=1,11), & |
| 77 | +! n=10 |
| 78 | + 1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1, (0, mm=1,10), & |
| 79 | +! n=11 |
| 80 | + 1, 11, 55, 165, 330, 462, 462, 330, 165, 55, 11, 1, (0, mm=1,9), & |
| 81 | +! n=12 |
| 82 | + 1, 12, 66, 220, 495, 792, 924, 792, 495, 220, 66, 12, 1, (0, mm=1,8), & |
| 83 | +! n=13 |
| 84 | + 1, 13, 78, 286, 715, 1287, 1716, 1716, 1287, 715, 286, 78, 13, 1, (0, mm=1,7), & |
| 85 | +! n=14 |
| 86 | + 1, 14, 91, 364, 1001, 2002, 3003, 3432, 3003, 2002, 1001, 364, 91, 14, 1, (0, mm=1,6), & |
| 87 | +! n=15 |
| 88 | + 1, 15, 105, 455, 1365, 3003, 5005, 6435, 6435, 5005, 3003, 1365, 455, 105, 15, 1, (0, mm=1,5), & |
| 89 | +! n=16 |
| 90 | + 1, 16, 120, 560, 1820, 4368, 8008, 11440, 12870, 11440, 8008, 4368, 1820, 560, 120, 16, 1, (0, mm=1,4), & |
| 91 | +! n=17 |
| 92 | + 1, 17, 136, 680, 2380, 6188, 12376, 19448, 24310, 24310, 19448, 12376, 6188, 2380, 680, 136, 17, 1, (0, mm=1,3), & |
| 93 | +! n=18 |
| 94 | + 1, 18, 153, 816, 3060, 8568, 18564, 31824, 43758, 48620, 43758, 31824, 18564, 8568, 3060, 816, 153, 18, 1, 0, 0, & |
| 95 | +! n=19 |
| 96 | + 1, 19, 171, 969, 3876, 11628, 27132, 50388, 75582, 92378, 92378, 75582, 50388, 27132, 11628, 3876, 969, 171, 19, 1, 0, & |
| 97 | +! n=20 |
| 98 | + 1, 20, 190, 1140, 4845, 15504, 38760, 77520, 125970, 167960, 184756, 167960, 125970, 77520, 38760, 15504, 4845, 1140, 190, 20, & |
| 99 | + 1],DP),[21,21])) |
| 100 | + |
| 101 | + |
| 102 | + PRIVATE :: DP |
| 103 | + |
| 104 | +Contains |
| 105 | + |
| 106 | +! Define some support routines for generating binomial coefficients etc. |
| 107 | + |
| 108 | + Pure Subroutine binomialCoefsArray(n, bincoefs) |
| 109 | + |
| 110 | +! Get array of all binomial coefficient for a given n . |
| 111 | +! bincoefs must be dimensioned (0:n), For n <= 20 the |
| 112 | +! binonialCoefsTo20 is used to create the array. For n>20 |
| 113 | +! the first 0:20 rows are taken from binomialCoefsTo20 and |
| 114 | +! the remaining rows are generated by recursion. |
| 115 | + |
| 116 | +! This should return the coefs for any n in Table 24.1 in Abramowitz and Stegun |
| 117 | +! "Handbook of Mathematical Functions" for a given n. The recursion |
| 118 | +! relationship in given in Abramowitz and Stegun Eq. 3.1.4 (chapter 3, page 10 |
| 119 | +! in my version) |
| 120 | + |
| 121 | + Implicit NONE |
| 122 | + |
| 123 | + Integer, Intent(IN) :: n |
| 124 | + Real(DP), Intent(INOUT) :: bincoefs(0:n) |
| 125 | + |
| 126 | + Integer :: i, j |
| 127 | + |
| 128 | + Real(DP) :: lastRow(0:n) |
| 129 | + |
| 130 | +! Initialize bincof array to zero and bincoef(0,0) to 1 |
| 131 | + |
| 132 | + bincoefs = 0.0_DP |
| 133 | + bincoefs(0) = 1.0_DP |
| 134 | + |
| 135 | +! For n <= 20 get values from table |
| 136 | + |
| 137 | + If (n <= 20) Then |
| 138 | + |
| 139 | + bincoefs(0:n) = binomialCoefsTo20(n,0:n) |
| 140 | + |
| 141 | + Else |
| 142 | + |
| 143 | +! Load lower triangle of binceofs array with row coefficients using the fact |
| 144 | +! that the value of a coefficient in the nth row is the sum of the two |
| 145 | +! values above and to the left of it in the preceeding (n-1) row. |
| 146 | + |
| 147 | + lastRow = 0.0_DP |
| 148 | + lastRow(0:20) = binomialCoefsTo20(20,:) |
| 149 | + |
| 150 | + Do i=21,n |
| 151 | + Do j=1,i |
| 152 | + bincoefs(j) = lastRow(j) + lastRow(j-1) |
| 153 | + EndDo |
| 154 | + bincoefs(0) = 1.0_DP |
| 155 | + bincoefs(i) = 1.0_DP |
| 156 | + lastRow(0:i) = bincoefs(0:i) |
| 157 | + EndDo |
| 158 | + EndIf |
| 159 | + |
| 160 | + End Subroutine binomialCoefsArray |
| 161 | + |
| 162 | + Pure Subroutine getBinomialCoefsTable(n, bincoefs) |
| 163 | + |
| 164 | +! Generate binomial coefficient table for a given n . bincoefs |
| 165 | +! must be dimensioned (0:n, 0:n), For n <= 20 the |
| 166 | +! binonialCoefsTo20 is used to create the table. For n>20 |
| 167 | +! the first 0:20 rows are taken from binomialCoefsTo20 and |
| 168 | +! the remaining rows are generated by recursion. |
| 169 | + |
| 170 | +! This should reproduce Table 24.1 in Abramowitz and Stegun |
| 171 | +! "Handbook of Mathematical Functions" for a given n. The recursion |
| 172 | +! relationship in given in Abramowitz and Stegun Eq. 3.1.4 (chapter 3, page 10 |
| 173 | +! in my version) |
| 174 | + |
| 175 | + Implicit NONE |
| 176 | + |
| 177 | + Integer, Intent(IN) :: n |
| 178 | + Real(DP), Intent(INOUT) :: bincoefs(0: , 0:) |
| 179 | + |
| 180 | + Integer :: i, j |
| 181 | + |
| 182 | +! Initialize bincof array to zero and bincoef(0,0) to 1 |
| 183 | + |
| 184 | + bincoefs = 0.0_DP |
| 185 | + bincoefs(0,0) = 1.0_DP |
| 186 | + |
| 187 | +! For n <= 20 get values from table |
| 188 | + |
| 189 | + If (n <= 20) Then |
| 190 | + Do i=1,n |
| 191 | + bincoefs(i,0:i) = binomialCoefsTo20(i,0:i) |
| 192 | + EndDo |
| 193 | + Else |
| 194 | + |
| 195 | +! Load lower triangle of binceofs array with row coefficients using the fact |
| 196 | +! that the value of a coefficient in the nth row is the sum of the two |
| 197 | +! values above and to the left of it in the preceeding (n-1) row. |
| 198 | + |
| 199 | + bincoefs(0:20, 0:20) = binomialCoefsTo20(:,:) |
| 200 | + Do i=20,n-1 |
| 201 | + Do j=1,i |
| 202 | + bincoefs(i+1,j) = bincoefs(i,j) + bincoefs(i,j-1) |
| 203 | + EndDo |
| 204 | + bincoefs(i+1,0) = 1.0_DP |
| 205 | + bincoefs(i+1,i+1) = 1.0_DP |
| 206 | + EndDo |
| 207 | + EndIf |
| 208 | + |
| 209 | + End Subroutine getBinomialCoefsTable |
| 210 | + |
| 211 | + Pure Function binomialCoef(n, j) Result(bc) |
| 212 | + |
| 213 | +! Get the single binomial coefficient defined at (n,j) |
| 214 | + |
| 215 | + Implicit NONE |
| 216 | + |
| 217 | + Integer, Intent(IN) :: n, j |
| 218 | + Real(DP) :: bc |
| 219 | + |
| 220 | + Real(DP) :: bincoefs(0:n) |
| 221 | + |
| 222 | + If (j>n .OR. j<0) Then |
| 223 | + bc = 0.0_DP |
| 224 | + RETURN |
| 225 | + EndIf |
| 226 | + |
| 227 | + If (j==n) Then |
| 228 | + bc = 1.0_DP |
| 229 | + ElseIf (j>n) Then |
| 230 | + bc = 0.0_DP |
| 231 | + ElseIf (j==0) Then |
| 232 | + bc = 1.0_DP |
| 233 | + ElseIf(n<=20) Then |
| 234 | + bc = binomialCoefsTo20(n,j) |
| 235 | + Else |
| 236 | + Call binomialCoefsArray(n, bincoefs) |
| 237 | + bc = bincoefs(j) |
| 238 | + EndIf |
| 239 | + |
| 240 | + End Function binomialCoef |
| 241 | + |
| 242 | +End Module binomial_coefficients |
0 commit comments