C[omp]ute

From: andrew cooke <andrew@...>

Date: Tue, 1 Jul 2025 15:27:44 -0400

I am writing some (audio) synthesis code for a CPU without fast
division.  Internally, amplitude is described with int16_t, but for
some processes (eg wavefolding) it's useful to switch to floats.
Without fast division this switch could be expensive.  So I started to
prematurely optimize...

I realised that one way to convert between the two was to directly
write to the mantissa part of an IEEE float - the float then
"automatically" has the correct value without any explicit conversion.
The only tricky technical details seemed to be handling the hidden bit
and getting the exponent correct (avoiding denormalized values).

But I ended up spending too much time on the details of type punning
(that C trick of writing and reading from a union using different
types).

First, it turns out that type punning is not supported in C++, but
seems to work in g++ (and C++23 has start_lifetime_as which can, I think,
be trivially used to replace the type punning in my final code below).

Second, the code at https://stackoverflow.com/a/15685301 DOES NOT
WORK.  Maybe it works in C - I don't know, I haven't tried it.  But in
C++ it fails.  The reason is described at
https://en.cppreference.com/w/cpp/language/bit_field.html - second
example, with the comment "6 bits for b2 - doesn't fit into the 1st
byte => starts a 2nd".

(If you have no idea what I'm talking about
https://en.wikipedia.org/wiki/IEEE_754 might help).

So you end up having a union of a float with an uint32_t and doing the
separation into the different components manually.  And that does
work!  Although I still have no idea if it's faster that division.

Here's the code.

header:

  class IEEEFloat {

  public:

    IEEEFloat(double f);
    IEEEFloat(float f);
    IEEEFloat(int v);
    IEEEFloat(int16_t v);
    IEEEFloat(uint32_t m, uint32_t e, uint32_t s);

    float f();
    uint32_t m();
    uint32_t e();
    uint32_t s();
    int16_t sample();
    void dump(std::ostream& c);

  private:

    const uint32_t hidden = 1 << 23;
    const uint32_t mask = hidden - 1;

    typedef union {
      float f;
      uint32_t u;
    } float_cast;

    float_cast fc;

  };


  float sample2float(int16_t s);

  int16_t float2sample(float f);


and implementation:

IEEEFloat::IEEEFloat(double v) : IEEEFloat(static_cast<float>(v)) {};

IEEEFloat::IEEEFloat(float v) : fc({.f=v}) {};

IEEEFloat::IEEEFloat(uint32_t m, uint32_t e, uint32_t s) : fc({.u=0}) {
  fc.u = ((s & 1) << 31) | (e & 255) << 23 | (m & mask);
}

IEEEFloat::IEEEFloat(int v) : IEEEFloat(static_cast<int16_t>(v)) {};

IEEEFloat::IEEEFloat(int16_t v) : fc({.u=0}) {
  if (v != 0) {
    uint32_t s = static_cast<uint32_t>(v < 0) << 31;
    uint32_t e = 127;
    uint32_t m = static_cast<uint32_t>(abs(v)) << 8;
    while (! (m & hidden)) {
      m = m << 1;
      e--;
    }
    fc.u = s | (e << 23) | (m & mask);
  }
}

float IEEEFloat::f() {
  return fc.f;
}

uint32_t IEEEFloat::s() {
  return (fc.u >> 31) & 1;
}

uint32_t IEEEFloat::e() {
  return (fc.u >> 23) & 255;
}

uint32_t IEEEFloat::m() {
  return fc.u & mask;
}

int16_t IEEEFloat::sample() {
  if (e()) {
    int16_t v = static_cast<int16_t>((m() | hidden) >> (8 + 127 - e()));
    if (s()) v = -v;
    return v;
  } else {
    return 0;
  }
}

void IEEEFloat::dump(std::ostream& c) {
  c << fc.f << " (" << std::hex << m() << ", " << std::dec << e() << ", " << s() << ")" << std::endl;
}


float sample2float(int16_t s) {
  return IEEEFloat(s).f();
}

int16_t float2sample(float f) {
  return IEEEFloat(std::max(-0.999969f, std::min(0.999969f, f))).sample();
}



TEST_CASE("IEEEFloat") {
  CHECK(sample2float(sample_max) == doctest::Approx(1).epsilon(0.001));
  CHECK(sample2float(sample_max/2) == doctest::Approx(0.5).epsilon(0.001));
  CHECK(sample2float(0) == doctest::Approx(0).epsilon(0.001));
  CHECK(sample2float(sample_min/2) == doctest::Approx(-0.5).epsilon(0.001));
  CHECK(sample2float(sample_min) == doctest::Approx(-1).epsilon(0.001));
  CHECK(float2sample(1.0) == doctest::Approx(sample_max).epsilon(0.001));
  CHECK(float2sample(0.5) == doctest::Approx(sample_max/2).epsilon(0.001));
  CHECK(float2sample(0.0) == doctest::Approx(0).epsilon(0.001));
  CHECK(float2sample(-0.5) == doctest::Approx(sample_min/2).epsilon(0.001));
  CHECK(float2sample(-1.0) == doctest::Approx(sample_min).epsilon(0.001));
  for (int32_t s = sample_min; s <= sample_max; s += 1234) {
    CHECK(float2sample(sample2float(s)) == doctest::Approx(s).epsilon(0.001));
  }
  for (float f = -1; f <= 1; f += 0.123) {
    CHECK(sample2float(float2sample(f)) == doctest::Approx(f).epsilon(0.001));
  }
}


where sample_max and sample_min are 32767 and -32767 respectively.

Andrew