Converting single-precision floating point numbers to double-precision for division

Question

Working as a High-Performance-Computing guy, we tend to default to single-precision floating point numbers (float or real) whenever possible. This is because you can perform more operations per second if each operation is individually faster to perform.

One of the more senior people I work with, however, always insists that (when accuracy is required) you should temporarily convert your single-precision data to double-precision in order to perform division. That is:

float a, b;
float ans = ((double)a)/((double)b);

or

real :: a, b, ans
ans = real(dble(a)/dble(b))

depending on the language you're working in. In my opinion, this looks really ugly, and to be honest I don't even know if the answer held in ans will be more accurate than if you had simply written ans = a/b in single-point precision.

Can someone tell me whether converting your numbers prior to arithmetic, specifically for performing division, will actually result in a more accurate answer? Would this be a language/compiler specific question, or would this be up to IEEE? With what number values would this accuracy improvement be most noticeable?

Any enlightening comments/answers would be much appreciated.

Answer 1

float ans = ((double)a)/((double)b);

This article demonstrates that ans is always the same as would be computed by a single-precision division for IEEE 754 arithmetics and FLT_EVAL_METHOD=0.

When FLT_EVAL_METHOD=1, the same property is also trivially true.

When FLT_EVAL_METHOD=2, I am not sure. It is possible that one might interpret the rules as meaning that the long double computation of a/b must first be rounded to double, then to float. In this case, it can be less accurate than directly rounding from long double to float (the latter produces the correctly rounded results, whereas the former could fail to do so in extremely rare cases, unless another theorem like Figueroa's applies and shows that this never happens).

Long story short, for any modern, reasonable floating-point computing platform (*), it is superstition that float ans = ((double)a)/((double)b); has any benefits. You should ask the senior people you refer to in your question to exhibit one pair a, b of values for which the result is different, not to mention more accurate. Surely if they insist that this is better it should be no trouble for them to provide one single pair of values for which it makes a difference.

(*) remember to use -fexcess-precision=standard with GCC to preserve your sanity

Answer 2

This depends greatly on what platform is being used.

An 80x86 (or a 1980s-era 8087) using non-SSE instructions performs all its arithmetic using 80-bit precision (long double or real*10). It is the "store" instruction which moves results from the numeric processor to memory which loses precision.

Unless it is a really bone-headed compiler, maximum precision should occur from

float a = something, b = something_else;
float ans = a/b;

since to perform the division, the single precision operands will be extended precision after loading and the result will be extended precision.

If you were doing something more intricate and wanted to maintain maximum precision, don't store intermediate results in smaller-sized variables:

float a, b, c, d;

float prod_ad = a * d;
float prod_bc = b * c;
float sum_both = prod_ad + prod_bc;   // less accurate

That gives a less precise result than doing it all at once since most compilers will produce code which keeps all the intermediates values at extended precision:

float a, b, c, d;

float sum_both = a * d + b * c;   // more accurate

Building on Eugeniu Rosca's example program:

#include "stdio.h"
void main(void)
{
    float a=73;
    float b=19;

    long double a1 = a;
    long double b1 = b;

    float ans1 = (a*a*a/b/b/b);
    float ans2 = ((double)a*(double)a*(double)a/(double)b/(double)b/(double)b);
    float ans3 = a1*a1*a1/b1/b1/b1;
    long double ans4 = a1*a1*a1/b1/b1/b1;

    printf ("plain:  %.20g\n", ans1);
    printf ("cast:   %.20g\n", ans2);
    printf ("native: %.20g\n", ans3);
    printf ("full:   %.20Lg\n", ans4);
}

provides, no matter the optimization level

plain:  56.716281890869140625
cast:   56.71628570556640625
native: 56.71628570556640625
full:   56.716285172765709289

This is showing that for trivial operations, there isn't much difference. However, changing the constants to be more of a precision challenge:

float a=0.333333333333333333333333;
float b=0.1;

provides

plain:  37.03704071044921875
cast:   37.037036895751953125
native: 37.037036895751953125
full:   37.037038692721614131

where the precision difference is displaying a more pronounced effect.

Answer 3

Yes, converting to double precision will give you better accuracy (or, shall I say, precision) in division. One could say that this is up to IEEE, but only because IEEE defines the formats and standards. doubles are inherently more precise than floats, with storage of numbers as well as division.

To answer your last question, this would be most noticeable with large a and small (less than 1) b, because then you end up with a very large quotient, in the range at which all floating point numbers are less granular.

Answer 4

Running this on x86 (GCC 4.9.3):

#include "stdio.h"
int main(int arc, char **argv)
{
    float a=73;
    float b=19;

    float ans1 = (a*a*a/b/b/b);
    float ans2 = ((double)a*(double)a*(double)a/(double)b/(double)b/(double)b);
    printf("plain: %f\n", ans1);
    printf("cast:  %f\n", ans2);
    return 0;
}

outputs:

plain: 56.716282
cast:  56.716286

The same operations in a Windows calculator return:

56.716285172765709287068085726782

Clearly, the second result has greater accuracy.