float to int conversion going wrong (even though the float is already an int)

Question

I was writing a little function to calculate the binomial coefficiant using the tgamma function provided by c++. tgamma returns float values, but I wanted to return an integer. Please take a look at this example program comparing three ways of converting the float back to an int:

#include <iostream>
#include <cmath>

int BinCoeffnear(int n,int k){
    return std::nearbyint( std::tgamma(n+1) / (std::tgamma(k+1)*std::tgamma(n-k+1)) );
}

int BinCoeffcast(int n,int k){
    return static_cast<int>( std::tgamma(n+1) / (std::tgamma(k+1)*std::tgamma(n-k+1)) );
}

int BinCoeff(int n,int k){
    return (int) std::tgamma(n+1) / (std::tgamma(k+1)*std::tgamma(n-k+1));
}

int main()
{
    int n = 7;
    int k = 2;
    std::cout << "Correct: " << std::tgamma(7+1) / (std::tgamma(2+1)*std::tgamma(7-2+1)); //returns 21
    std::cout << " BinCoeff: " << BinCoeff(n,k); //returns 20 
    std::cout << " StaticCast: " << BinCoeffcast(n,k); //returns 20
    std::cout << " nearby int: " << BinCoeffnear(n,k); //returns 21

    return 0;   
}

why is it, that even though the calculation returns a float equal to 21, 'normal' conversion fails and only nearbyint returns the correct value. What is the nicest way to implement this?

EDIT: according to c++ documentation here tgamma(int) returns a double.

Answer 1

Floating-point numbers have rounding errors associated with them. Here is a good article on the subject: What Every Computer Scientist Should Know About Floating-Point Arithmetic.

In your case the floating-point number holds a value very close but less than 21. Rules for implicit floating–integral conversions say:

The fractional part is truncated, that is, the fractional part is discarded.

Whereas std::nearbyint:

Rounds the floating-point argument arg to an integer value in floating-point format, using the current rounding mode.

In this case the floating-point number will be exactly 21 and the following implicit conversion would return 21.

The first cout outputs 21 because of rounding that happens in cout by default. See std::setprecition.

Here's a live example.

---

What is the nicest way to implement this?

Use the exact integer factorial function that takes and returns unsigned int instead of tgamma.

Answer 2

The remark by @nos is on point. Note that the first line

std::cout << "Correct: " << 
 std::tgamma(7+1) / (std::tgamma(2+1)*std::tgamma(7-2+1));

Prints a double value and does not perform a floating point to integer conversion.

The result of your calculation in floating point is indeed less than 21, yet this double precision value is printed by cout as 21.

On my machine (x86_64, gnu libc, g++ 4.8, optimization level 0) setting cout.precision(18) makes the results explicit.

Correct: 20.9999999999999964  BinCoeff: 20 StaticCast: 20 nearby int: 21

In this case practical to replace integer operations with floating point operations, but one has to keep in mind that the result must be integer. The intention is to use std::round.

The problem with std::nearbyint is that depending on the rounding mode it may produce different results.

std::fesetround(FE_DOWNWARD);
std::cout << " nearby int: " << BinCoeffnear(n,k); 

would return 20.

So with std::round the BinCoeff function might look like

int BinCoeffRound(int n,int k){
            return static_cast<int>(
                            std::round(
                                    std::tgamma(n+1) /
                                       (std::tgamma(k+1)*std::tgamma(n-k+1)) 
                                  ));
}

Answer 3

From this std::tgamma reference:

If arg is a natural number, std::tgamma(arg) is the factorial of arg-1. Many implementations calculate the exact integer-domain factorial if the argument is a sufficiently small integer.

It seems that the compiler you're using is doing that, calculating the factorial of 7 for the expression std::tgamma(7+1).

The result might differ between compilers, and also between optimization levels. As demonstrated by Jonas there is a big difference between optimized and unoptimized builds.

Answer 4

the problem is on handling the floats. floats cant 2 as 2 but as 1.99999 something like that. So converting to int will drop out the decimal part.

So instead of converting to int immediately first round it to by calling the ceil function w/c declared in cmath or math.h.

this code will return all 21

#include <iostream>
#include <cmath>

int BinCoeffnear(int n,int k){
    return std::nearbyint( std::tgamma(n+1) / (std::tgamma(k+1)*std::tgamma(n-k+1)) );
}
int BinCoeffcast(int n,int k){
    return static_cast<int>( ceil(std::tgamma(n+1) / (std::tgamma(k+1)*std::tgamma(n-k+1))) );
}
int BinCoeff(int n,int k){
    return (int) ceil(std::tgamma(n+1) / (std::tgamma(k+1)*std::tgamma(n-k+1)));
}
int main()
{
    int n = 7;
    int k = 2;
    std::cout << "Correct: " << (std::tgamma(7+1) / (std::tgamma(2+1)*std::tgamma(7-2+1))); //returns 21
    std::cout << " BinCoeff: " << BinCoeff(n,k); //returns 20 
    std::cout << " StaticCast: " << BinCoeffcast(n,k); //returns 20
    std::cout << " nearby int: " << BinCoeffnear(n,k); //returns 21
    std::cout << "\n" << (int)(2.9995) << "\n";
}