Detect FPU rounding mode on a GPU

Question

I was delving into multi-precision arithmetics, and there is a nice fast class of algorithms, described in Jonathan Richard Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", 1997, Discrete & Computational Geometry, pages: 305–363. However, these algorithms rely on the FPU using round-to-even tiebreaking.

On CPU, it would be easy, one would just check or set the FPU state word and would be sure. However, there is no such instruction (yet?) for GPU programming.

That is why I was wondering if there is a dependable way of detecting (not setting) the rounding mode on an unknown FPU, perhaps by calculating a couple of tests and looking at the bit patterns of the resulting floating point numbers?

EDIT:

Just to summarize, the accepted code indeed seems to work, you can try:

#include <stdio.h>
#include <stdlib.h>
#include <float.h> // _controlfp()
#include <stdint.h>

int is_round_to_nearest()
{
    union {
        double f;
        uint64_t n;
    } special;

    special.n = 0 | (((uint64_t)(-0x100 + 1023) & 0x7ff) << 52) | 0; // no sign, 1.0 mantissa is expressed as zeroes, the 1 is implicit
    //const double special.f = atof("0x1.0p-100");

    if( 1.0 + special.f !=  1.0)
        return 0;
    if( 1.0 - special.f !=  1.0)
        return 0;
    if(-1.0 + special.f != -1.0)
        return 0;
    if(-1.0 - special.f != -1.0)
        return 0;
    return 1;
}

void main()
{
    printf("default : %d\n", is_round_to_nearest());
    _controlfp(_RC_CHOP, _MCW_RC);
    printf("_RC_CHOP : %d\n", is_round_to_nearest());
    _controlfp(_RC_UP, _MCW_RC);
    printf("_RC_UP : %d\n", is_round_to_nearest());
    _controlfp(_RC_DOWN, _MCW_RC);
    printf("_RC_DOWN : %d\n", is_round_to_nearest());
    _controlfp(_RC_NEAR, _MCW_RC);
    printf("_RC_NEAR : %d\n", is_round_to_nearest());
}

And that gives the following answers:

default : 1
_RC_CHOP : 0
_RC_UP : 0
_RC_DOWN : 0
_RC_NEAR : 1

Note that I was unable to set round to nearest, ties away from zero mode on my machine. In Visual Studio, floating point mode needs to be set to strict (/fp:strict), otherwise this will not work in release mode (all modes identified as nearest).

The following code seems to work even in release, even with the default or fast (/fp:precise, /fp:fast) rounding modes, but there is still no guarantee how will your compiler optimize the code:

int is_round_to_nearest()
{
    union {
        double f;
        uint64_t n;
    } special;

    special.n = 0 | (((uint64_t)(-0x100 + 1023) & 0x7ff) << 52) | 0; // no sign, 1.0 mantissa is expressed as zeroes, the 1 is implicit
    //const double special.f = atof("0x1.0p-100");

    volatile double v;
    v = 1.0; v += special.f;
    if(v !=  1.0)
        return 0;
    v = 1.0; v -= special.f;
    if(v !=  1.0)
        return 0;
    v = -1.0; v += special.f;
    if(v != -1.0)
        return 0;
    v = -1.0; v -= special.f;
    if(v != -1.0)
        return 0;
    return 1;
}

Answer 1

This C code tells you that you are either in round-to-nearest-even or using a strange floating-point architecture indeed:

int is_round_to_nearest(void)
{
  if ( 1.0 + 0x1.0p-100 !=  1.0) return 0;
  if ( 1.0 - 0x1.0p-100 !=  1.0) return 0;
  if (-1.0 + 0x1.0p-100 != -1.0) return 0;
  if (-1.0 - 0x1.0p-100 != -1.0) return 0;
  return 1;
}

You can add an f suffix to all twelve floating-point constants above to obtain a single-precision function.