Optimization of modified brownian noise

test.c

// Optimization of modified brownian noise: trying to get rid of branching in
// the ‘reference’ function.

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static inline float
reference (const float a, const float b)
{
  float aplusb = a + b;
  if ((aplusb <= -1.0f) || (aplusb >= 1.0f))
    return a-b;
  return aplusb;
}

static inline float
evil (const float a, const float b)
{
  union { float f; uint32_t i; } u, ub;
  u.f = a + b;

  // Only keep the exponent, which is n+127 where n<0 if |a+b|<1.0.
  u.i &= 0xff<<23;
  // Normalize the exponent. The sign will match in both representations.
  u.i -= 127<<23;
  // Only keep the sign. Minus if |a+b|<1.0.
  u.i &= 1<<31;

  // Flip ub’s sign based on u, negating ub if |a+b|<1.0.
  ub.f = b;
  ub.i ^= u.i;

  return a - ub.f;
}

static inline float
less_evil (const float a, const float b)
{
  float c[2] = { a + b, a - b };
  return c[(c[0] <= -1.0f) || (c[0] >= 1.0f)];
}

static inline float randf (void) { return (float)rand () / (float)RAND_MAX; }

int
main (int argc, char **argv)
{
  (void)argc;

  int n = 100000000;
  float out;

  char **arg = &argv[1];
  if (! *arg) {
    fprintf (stderr, "USAGE: %s [reference|evil|less_evil|verify]...\n",
             argv[0]);
    return 1;
  }

  for (; *arg; arg++) {
    out = 0.0f;
    srand (1);

    if (! strcmp (*arg, "reference")) {
      for (int i = 0; i < n; i++) {
        float white = 2.0f*randf () - 1.0f;
        out = reference (out, 0.1f*white);
      }

    } else if (! strcmp (*arg, "evil")) {
      for (int i = 0; i < n; i++) {
        float white = 2.0f*randf () - 1.0f;
        out = evil (out, 0.1f*white);
      }

    } else if (! strcmp (*arg, "less_evil")) {
      for (int i = 0; i < n; i++) {
        float white = 2.0f*randf () - 1.0f;
        out = less_evil (out, 0.1f*white);
      }

    } else if (! strcmp (*arg, "verify")) {
      for (int i = 0; i < INT32_MAX; i++) {
        float a = 2.0f*randf () - 1.0f,
              b = 2.0f*randf () - 1.0f;

        float c0 = reference (a, b),
              c1 = evil (a, b),
              c2 = less_evil (a, b);

        if ((c0 != c1) || (c0 != c2) || (c0 < -1.0f) || (c0 > 1.0f))
          printf ("%d: a: %f, b: %f, c0: %f, "
                  "c1: %f, diff: %.10f, c2: %f, diff: %.10f\n",
                  i, a, b, c0, c1, (c1-c0), c2, (c2-c0));
      }

    } else {
      return 1;
    }
  }

  return 0;
}

// vim:set et sw=2 sts=2: