Declared double precision Function in F90 and return type mismatch error

If one is computing many binomial coefficients, then it might be reasonable to use this expression where factorial(:) is a precomputed array of factorial values. The problem with this approach, whether computed with integers or with floating point numbers, is the danger of unnecessary overflow. For int32, the largest value that can be stored is 12!, and for int64 the largest value is 20!. For real32, the largest value that can be stored without overflowing the exponent is 34!, and in real64 the largest value is 170!.

With the expressions given above in this thread, binomial values with larger values of n can be computed correctly without overflow. This is because the values are computed with both the multiplications in the numerator and the divisions in the denominator done together, thus keeping the intermediate values as small as possible.

I’m not sure these expressions are the result of clever programming. The expression above in this thread is one of the traditional ways to compute these values by hand. It is based on the identity n.choose.r=((n-1).choose.(r-1))*n/r. When computed with integers, you want to evaluate the numerator products first before dividing by the denominators because you don’t need to worry about remainders with that operation order. I guess the floating point values should be computed that way for the same reason, although the floating point computations will eventually have roundoff errors for large n regardless. Doing the floating point divisions first allows some values to be computed without overflow, but at the expense of introducing roundoff errors, so there are tradeoffs involved with either expression order.

I guess it should be mentioned in this context that the binomial coefficients can be computed as Pascal’s triangle using at most one addition per value (no multiplications, no divisions). If the range of n is small enough, then the whole triangle can be computed and stored, so only an array lookup is needed for each value. There are some applications where this is indeed the optimal approach for computing and accessing the values.

1 Like