optimization/10644: float calculation exceeds representation max, only when optimizer is turned on

[bird at bainet dot com]

>Number:         10644
>Category:       optimization
>Synopsis:       float calculation exceeds representation max, only when optimizer is turned on
>Confidential:   no
>Severity:       non-critical
>Priority:       medium
>Responsible:    unassigned
>State:          open
>Class:          wrong-code
>Submitter-Id:   net
>Arrival-Date:   Tue May 06 16:46:00 UTC 2003
>Closed-Date:
>Last-Modified:
>Originator:     Rebecca Bird
>Release:        gcc version 3.2.3 (also 3.2.1 and 3.2)
>Organization:
>Environment:
Red Hat Linux 7.2, kernel 2.4.18, on Pentium IV
>Description:
#include <iostream>			// ostream
#include <limits>		      	// numeric_limits<T>::max
using namespace std; 

// -------------------------------------------------------------------
template<class T>			    // T = double or float
static T factorial(T value)
{
    T fact = 1.0;
    while (value > 1.0) 
        fact *= value--;	
    return fact;
}

// -------------------------------------------------------------------
template<class T>			// T = double or float
static void test(T)
{
    cout << "representation max = " << numeric_limits<T>::max() << "\n\n";

    cout << "factorial(34) = "<< factorial(T(34)) << endl;	
    cout << "factorial(35) = " << factorial(T(35)) << endl;	
    cout << "factorial(170) = " << factorial(T(170)) << endl;	
    cout << "factorial(171) = " << factorial(T(171)) << "\n\n";	
}

// -------------------------------------------------------------------
int main()
{
    cout << "factorial limits for floats:\n";
    test((float)0);
    cout << "factorial limits for doubles:\n";
    test((double)0);
    return 0;
}
>How-To-Repeat:
The code explores the range of the factorial function for floats vs doubles.  When run without optimization:
    g++-3.2  factorial.cpp -o test
the code correctly outputs "inf" for arguments > 34 for float instantiation, and "inf" for arguments > 170 for double.

When any optimization level is used, say:
    g++-3.2  -O3 factorial.cpp -o test
the double output remains the same, but factorial(float 35) and factorial(float 170) now give incorrect results.  Specifically, the float instantiation returns values 1.03331e+40 for 35! and 7.25742e+306 for 170! which are more than the max float = 2.95233e+38.

>Fix:
If I add code to factorial such as:
    T fact = 1.0;
    ostringstream ist;		
    ist <<  fact;	
the program produces correct answers even with optimization turned on.  However, this slows down the processing time for the factorial routine (and shouldn't be necessary).
>Release-Note:
>Audit-Trail:
>Unformatted:
----gnatsweb-attachment----
Content-Type: text/plain; name="factorial.cpp"
Content-Disposition: inline; filename="factorial.cpp"

// Show limits of factorial routine for floats vs doubles

#include <iostream>			// ostream
#include <limits>			// numeric_limits<T>::max
using namespace std; 

// -------------------------------------------------------------------
template<class T>			// T = double or float
static T factorial(T value)
{
	T fact = 1.0;
	while (value > 1.0) 
		fact *= value--;	
	return fact;
}

// -------------------------------------------------------------------
template<class T>			// T = double or float
static void test(T)
{
	cout << "representation max = " << numeric_limits<T>::max() << "\n\n";

	cout << "factorial(34) = "<< factorial(T(34)) << endl;	
	cout << "factorial(35) = " << factorial(T(35)) << endl;	
	cout << "factorial(170) = " << factorial(T(170)) << endl;	
	cout << "factorial(171) = " << factorial(T(171)) << "\n\n";	
}

// -------------------------------------------------------------------
int main()
{
	cout << "factorial limits for floats:\n";
	test((float)0);
	cout << "factorial limits for doubles:\n";
	test((double)0);
	return 0;
}