# Software floating-point multiplication

I wrote a floating-point multiplication function as an excercise. The program compares its result to the usual hardware multiplication result and for this purpose I use unspecified behavior, but the function itself should be fine. I do get warnings about XOR-ing 1-bit ints but I have no idea why.

I noticed there aren't that many comments in the multiplication function so I wonder what I could put there.

#include <math.h>
#include <assert.h>
#include <memory.h>
#include <limits.h>
#include <fenv.h>
#include <float.h>

#include <mutex>
#include <iostream>
#include <random>

int msb( uint64_t v )
{
if(v == 0) return 0;
return 63-__builtin_clzll(v);
}

void roundedShift(uint64_t* a, unsigned int n) {
if(n >= sizeof(uint64_t)*CHAR_BIT) {
// we don't care about rounding if n == 64 because mantissa product is always less than 2^46
*a = 0;
}
uint64_t roundoff = *a - ((*a >> n) << n);
*a >>= n;
if(roundoff*2 > (1ull<<n)
||
(roundoff*2 == (1ull<<n) && (*a)%2 == 1)) {
++(*a);
}
}

struct SoftFloat {
unsigned int mantissa : 23;
unsigned int exponent : 8;
unsigned int sign : 1;

static const SoftFloat nan;
SoftFloat operator *(SoftFloat right) const {
if(exponent == 255 || right.exponent == 255) {
return specialMultiplication(right);
}
uint64_t fullLeftMantissa = (mantissa)+(exponent!=0?(1<<23):0);
uint64_t fullRightMantissa = (right.mantissa)+(right.exponent!=0?(1<<23):0);
uint64_t fullNewMantissa = fullLeftMantissa * fullRightMantissa;
// experiment have shown that operating on biased exponents is faster than
// computing unbiased exponent and then adding bias
short leftNormalBiasedExponent = exponent != 0 ? exponent: 1;
short rightNormalBiasedExponent = right.exponent != 0 ? right.exponent : 1;
short newNormalBiasedExponent = leftNormalBiasedExponent +
rightNormalBiasedExponent - 127 - 23;

int totalShift = 0;
if(newNormalBiasedExponent < 1) {
int diff = -newNormalBiasedExponent + 1;
newNormalBiasedExponent = 1;
totalShift += diff;
}
int implicitBit = msb(fullNewMantissa);
int shift = implicitBit - totalShift - 23;
if(shift >= 0) {
newNormalBiasedExponent += shift;
totalShift += shift;
fullNewMantissa &= ~(1ll << implicitBit);
} else {
newNormalBiasedExponent = 0;
}
roundedShift(&fullNewMantissa, totalShift);
if(fullNewMantissa == (1 << 23)) {
++newNormalBiasedExponent;
fullNewMantissa = 0;
}
if(newNormalBiasedExponent >= 255) {
newNormalBiasedExponent = 255;
fullNewMantissa = 0;
}
return SoftFloat{(uint)fullNewMantissa, (uint)newNormalBiasedExponent, sign ^ right.sign};
// i don't really understan why this return expression gives warning
// narrowing conversion of ‘(((int)((const SoftFloat*)this)->SoftFloat::sign) ^ ((int)right.SoftFloat::sign))’ from ‘int’ to ‘unsigned int’
}
bool isNan() const {
return exponent == 255 && mantissa != 0;
// This is same as (*reinterpret_cast<const unsigned int*>(this) << 1) > 0b11111111000000000000000000000000u,
// but the experiments have shown it doesn't make any difference for speed on my machine.
// Current version does not invoke any UB
}
bool isZero() const {
return exponent == 0 && mantissa == 0;
}

bool isRepresentationEqual(const SoftFloat& right) {
return !memcmp(this, &right, sizeof(SoftFloat));
}

// Correctness testing relies on specific layout of floats and of fields inside SoftFloat, but the multiplication function itself does not
float toHardFloat() const {
float res;
memcpy(&res, this, sizeof(SoftFloat));
return res;
}
static SoftFloat fromOrdinalNumber(const unsigned int& a) {
SoftFloat res;
memcpy(&res, &a, sizeof(SoftFloat));
return res;
}
static SoftFloat fromHardFloat(const float& a) {
SoftFloat res;
memcpy(&res, &a, sizeof(SoftFloat));
return res;
}
private:
SoftFloat specialMultiplication(SoftFloat right) const {
// precondition - at least one of *this, right is either inf, -inf or nan
if(isNan() || right.isNan()) {
return nan;
}
if(isZero() || right.isZero()) {
return nan;
}
return {0, 255, sign ^ right.sign};
}
};
const SoftFloat SoftFloat::nan = {1, 255, 0};

std::mutex coutMutex;
void checkOneMultiplication(SoftFloat left, SoftFloat right) {
SoftFloat softProduct = left*right;
SoftFloat hardProduct = SoftFloat::fromHardFloat(
left.toHardFloat() *
right.toHardFloat());
bool equalRepresentation = softProduct.isRepresentationEqual(hardProduct);
bool bothNan = softProduct.isNan() &&
hardProduct.isNan();
if(!(equalRepresentation || bothNan)) {
std::lock_guard<std::mutex> lock(coutMutex);
std::cerr << "failed\n";
std::cerr << "left operand\t" << left.mantissa << "  " << left.exponent << "  " << left.sign << "  " << left.toHardFloat() << std::endl;
std::cerr << "right operand\t" << right.mantissa << "  " << right.exponent << "  " << right.sign << "  " << right.toHardFloat() << std::endl;
std::cerr << "soft product\t" << softProduct.mantissa << "  " << softProduct.exponent << "  " << softProduct.sign << "  " << softProduct.toHardFloat() << std::endl;
std::cerr << "hard product\t" << hardProduct.mantissa << "  " << hardProduct.exponent << "  " << hardProduct.sign << "  " << hardProduct.toHardFloat() << std::endl;
abort();
}
}

void checkRange(unsigned int start, unsigned int end, int threadNumber) {
// this loop checks multiplication with every one of 2^28 possible floats in given range of representaions. Right operand is pseudorandom but deterministic
std::uniform_int_distribution<unsigned int> distr;
for(unsigned int representationNumber = start; representationNumber != end; representationNumber++) {
// can't use < in condition because last block ends with 0
int rigntOperandRepresentationNumber = distr(gen);
SoftFloat leftOperand = SoftFloat::fromOrdinalNumber(representationNumber);
SoftFloat rightOperand = SoftFloat::fromOrdinalNumber(rigntOperandRepresentationNumber);
checkOneMultiplication(leftOperand, rightOperand);
if(representationNumber%0x10000000 == 0 && representationNumber > start) {
std::lock_guard<std::mutex> lock(coutMutex);
std::cout << "thread " << threadNumber << " checked 0x" << std::hex << representationNumber-start << " multiplications" << std::endl;
}
}
}

int main(int argc, char *argv[])
{
static_assert(sizeof(SoftFloat) == sizeof(float));
static_assert(alignof(SoftFloat) == alignof(float), "");
static_assert(sizeof(SoftFloat) == sizeof(int), "");
static_assert(__BYTE_ORDER__ == LITTLE_ENDIAN, "");
// Unfortunately we can't statically check further details of layout of structures

assert(SoftFloat::fromHardFloat(std::numeric_limits<float>::denorm_min()).isRepresentationEqual(SoftFloat{1, 0, 0}));
assert(SoftFloat::fromHardFloat(FLT_MIN).isRepresentationEqual(SoftFloat{0, 1, 0}));

float nonConstant = -1;
assert(SoftFloat::fromHardFloat(nonConstant*0).isRepresentationEqual(SoftFloat{0, 0, 1}));
// Test can't run if any of the above asserts fail

fesetround(FE_TONEAREST); // assuming mantissa is rounded to even when there are two nearest

const unsigned int blockSize = UINT32_MAX/4+1;

for(int i = 1; i < threadCount; i++) {
}
checkRange(0, blockSize, 0);

SoftFloat specialCases[] = {{0, 0, 0}, // 0
{0, 0, 1}, // -0
{0, 255, 0}, //inf
{0, 255, 1}, //-inf
SoftFloat::nan,
SoftFloat::fromHardFloat(std::nanf("")),
SoftFloat::fromHardFloat(FLT_MIN),
SoftFloat::fromHardFloat(1),
SoftFloat::fromHardFloat(FLT_MAX),
SoftFloat::fromHardFloat(FLT_EPSILON),
SoftFloat::fromHardFloat(std::numeric_limits<float>::denorm_min()),
SoftFloat::fromHardFloat(-FLT_MAX),
SoftFloat::fromHardFloat(FLT_MIN_EXP),
SoftFloat::fromHardFloat(FLT_MIN_10_EXP),
SoftFloat::fromHardFloat(FLT_MAX_EXP),
SoftFloat::fromHardFloat(FLT_MAX_10_EXP)};
int numberOfSpecialCases = sizeof(specialCases)/sizeof(SoftFloat);
for(int i = 0; i < numberOfSpecialCases; i++) {
for(int j = 0; j < numberOfSpecialCases; j++) {
checkOneMultiplication(specialCases[i], specialCases[j]);
}
}
for(int i = 1; i < threadCount; i++) {
}

std::cout << "checked 0x" << UINT32_MAX + numberOfSpecialCases*numberOfSpecialCases << " multiplications";
return 0;
}


The program tests multiplication of every possible floating-point number with a random number. This takes a long time to run, so I made it threaded. I also used fancy C++11 RNG stuff to make it deterministic.

The function uses __builtin_clzll which is available on gcc and clang but not on MSVC, I believe.

• What (compiler, OS, standard) are you targeting? Sep 22, 2019 at 1:22
• Let's say gcc, linux, x64 Sep 23, 2019 at 19:00

Fairly good code overall.

Mostly small stuff below.

Narrowing

I do get warnings about XOR-ing 1-bit ints but I have no idea why.

Weak compiler.

Perhaps use bool sign : 1;

Naked magic numbers

Rather than 23, 127, etc, consider the C-ish

#define MANTISSA_BIT_WIDTH 23


or a C++ -ish

const int mantissa_bit_width = 23;


uint?

uint appears non-standard. Perhaps unsigned?

// return SoftFloat{(uint)fullNewMantissa, (uint)newNormalBiasedExponent, sign ^ right.sign};
return SoftFloat{(unsigned)fullNewMantissa, (unsigned)newNormalBiasedExponent, sign ^ right.sign};


Minor stuff

Portability

Although OP has "gcc, linux, x64", little changes would step toward portability without sacrificing efficient emitted code.

// if(fullNewMantissa == (1 << 23)) {
if(fullNewMantissa == (1ul << 23)) {  // int could be 16-bit

// fullNewMantissa &= ~(1ll << implicitBit);
fullNewMantissa &= ~(1ull << implicitBit);  // Why mess with signed shifts?


int vs. short

Rarely is short faster/better than int unless one has an array of the type.

Consider

// short leftNormalBiasedExponent, rightNormalBiasedExponent, newNormalBiasedExponent
int leftNormalBiasedExponent, rightNormalBiasedExponent, newNormalBiasedExponent


sizeof type vs sizeof object

Consider the clearer, less maintenance

// return !memcmp(this, &right, sizeof(SoftFloat));
return !memcmp(this, &right, sizeof *this);


roundedShift()

Unclear about roundedShift() correctness. Partly due to lack of comments, partly due to "it takes time" to analyze.

• I doubt weak compiler - the 1-bit unsigned are each promoted to int for the XOR, so narrowing to 1-bit unsigned warrants a warning. Changing the type to bool and using != in place of ^ is a good suggestion, though. Sep 30, 2019 at 10:26

Using the C compatibility headers is questionable in new code:

#include <math.h>
#include <assert.h>
#include <memory.h>
#include <limits.h>
#include <fenv.h>
#include <float.h>


I recommend using the C++ versions (<cmath> etc.) which define their identifiers in the std namespace.

We're missing an include of <stdint.h> for uint64_t (or better, include <cstdint> to define std::uint64_t). Do we really need a 64-bit type, or would std::fast_uint64_t be a better choice?

There are many instances of sizeof with a typename argument, which would be clearer with an ordinary expression argument. For example:

    SoftFloat res;
memcpy(&res, &a, sizeof(SoftFloat));
return res;


Here we can show that we're correctly passing the size of the destination argument:

    SoftFloat res;
std::memcpy(&res, &a, sizeof res);
return res;