// Test my cube root routine for all 32-bit floats to see if it fails
// badly for some inputs.  I think this should take about 20 minutes
// on my palmtop.  ^\ (SIGQUIT) produces a status report saying how
// far it's gotten.  ^C (SIGINT) gives up early.

#include <unistd.h>
#include <string.h>
#include <stdio.h>
#include <signal.h>
#include <stdint.h>
#include <math.h>
#include <assert.h>


volatile uint32_t counter;

// convert to non-nul-terminated bytes in a signal-handler-safe way
static char *
u32toa(uint32_t n, char *s, int digits)
{
  uint32_t m = n/10;
  if (m) {
    s = u32toa(m, s, digits + 1);
    if (2 == digits % 3) *s++ = '_';
  }
  *s = '0' + n % 10;
  return s + 1;
}

// signal handler to report status
static void
status_handler(int sig)
{
  static char buf[16];
  char *end = u32toa(counter, buf, 0);
  *end++ = '\n';
  write(1, buf, end - buf);
}

volatile int stop = 0;
static void
sigint_handler(int sig)
{
  signal(SIGINT, SIG_DFL); // just die with a second ^C if hung
  stop = 1;
  status_handler(0);
}

float cuberoot(float y);  // from floatcbrt.c

int
main(int argc, char **argv)
{
  double worst_err = 0, worst_x = 0/0.0;
  uint32_t worst = 0xdeadbeef, last_checkpoint = 65536;
  assert(sizeof(float) == sizeof counter);

  signal(SIGQUIT, status_handler);
  signal(SIGINT, sigint_handler);
  
  do {
    float y;
    uint32_t n = counter;
    // https://stackoverflow.com/questions/4328342/float-bits-and-strict-aliasing
    memcpy(&y, &n, sizeof y);
    float x = cuberoot(y);
 
   if (!isnan(x) && !isnan(y)) {
      double err = fabs(y - x*x*x);

      if (err > worst_err) {
        worst_err = err;
        worst_x = x;
        worst = counter;
      }
    }
 
    counter++;

    if (counter / 2u >= last_checkpoint || stop || !counter) {
      last_checkpoint = counter;
      float worst_y, last_y;
      memcpy(&worst_y, &worst, sizeof worst);
      memcpy(&last_y, &last_checkpoint, sizeof last_y);
      printf("At %d (%g), worst error: got ∛%g (0x%08x, %A) ≈ %g : err %g\n",
             counter, last_y, worst_y, worst, worst_y, worst_x, worst_err);
    }
  } while (counter && !stop);

  return 0;
}