CODEWITHSUNDEEP
c++c++11floating-pointnumericieee-754
Frank
Frank

Reputation: 4417

How to make my C++ library compliant with IEEE 754

I have written a small library (in C++11) to compute with quaternions, but I quickly realized there are many precision/accuracy issues given that quite a few floating point operations are involved. That led me to looking at std::complex and IEEE 754. std::complex does a fair amount of work dealing with NaNs in particular. Is that the only thing I need to worry about to become "IEEE 754 compliant"? More generally, is there a "recipe" for how to turn an otherwise naive numerical library into one that would be "IEEE 754 compliant"?

Note: the problem is not about whether the compiler is IEEE 754, but about whether I should take special steps in my own algorithms to satisfy IEEE 754. After all, I am in the same position as whoever wrote std::complex, and they did take extra steps, e.g. for NaNs.

Upvotes: 1

Views: 685

Answers (3)

PavDub
PavDub

Reputation: 181

[@Frank : comments below my last answer Aug 3 at 7:38]

When multiplying comlex numbers, NaN can be a result of - using the notation of XCode iplementation at uploaded by @Frank - PasteBin:

__x is NaN if any of __a, __b, __c, __d is Nan or
__ac == +INF && __bd == +INF or
__ac == -INF && __bd == -INF or
__a == (+/-)INF && __c == 0 or *vice versa* or
__b == (+/-)INF && __d == 0 or *vice versa* or
    
__y is NaN if any of __a, __b, __c, __d is Nan or
__ad == -INF && __bc == +INF or
__ad == +INF && __bc == -INF or
__a == (+/-)INF && __d == 0 or *vice versa* or
__b == (+/-)INF && __c == 0 or *vice versa* or

Then they (XCode STL authors) play around with these to allow for cases such as

(INF + INFi)*(c + di)
(INF + INFi)*(0 + 0i)
etc...

The trick (most likely) is to set the INFs- and NaNs-containing values to either (+/-)1 or (+/-)0 so that the final recalculation

if (__recalc)
{
    __x = _Tp(INFINITY) * (__a * __c - __b * __d);
    __y = _Tp(INFINITY) * (__a * __d + __b * __c);
}

results in

INF * (+/-)INF = (+/-)INF or
INF * (+/-)0   = NaN

as appropriate.

This is the

how you translate them to the lib user

ps.> I do not have time to go so deep into the problem so that I could judge their implementation but it is most likely reasonable. (this Math forum link is helpful but I am of no help to readily translate it into the [a+bi] complex numers notation) Anyway, MinGW 4.9.1. 32bit Win implementation leaves this task open and the user should treat the corner cases. You can adopt the MinGW strategy at first or you can take that of XCode if you know how to treat the corner cases of the quaterions algebra properly. In any case this is no longer much about How to make my C++ library compliant with IEEE 754 ...

Upvotes: 0

PavDub
PavDub

Reputation: 181

Making the arithmetic IEEE-compliant is not your task but that of the compiler.

Your compiler either is or is not IEEE-compliant (@Paul) and you can do nothing much about it. Specifically, you cannot write a non-IEEE compliant program using an IEEE-compliant compiler double arithmetic (please don't take me at my word here :) ). This is just similar as you can hardly write a non-ISO/IEC 14882:2011-comliant program with an ISO/IEC 14882:2011-comliant compiler.

This however does not secure that the program gives an expected output and it is free of weird operations, including UB. And that is similar to where NaNs and INFs come into play in IEEE 754. They are there to give a “reasonable answer” to not well defined mathematical operations such 0/0 or out of range operations such as 1/0.

The standard does not care of how you treat these answers in your library (i.e. how you translate them to the lib user). The standard only specifies a way of how using these “reasonable answers” in yet another mathematical operation gives yet another “reasonable answer”.

You may be interested in http://perso.ens-lyon.fr/jean-michel.muller/goldberg.pdf (rather tough reading but worth the effort at least for gripping the idea)

Ps I checked my lib < complex > implementation (MinGW 4.9.1. 32bit) and did not notice any NaNs magic there.

Upvotes: 0

Paul Evans
Paul Evans

Reputation: 27577

You can simply check your library's Header <limits> constant numeric_limits::is_iec559 (IEC 559 is the same as IEEE 754) to see whether or not your C++ compiler already supports IEEE 754.

Upvotes: 1

Related Questions