#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <x86intrin.h>
#include <math.h>

#define ERROR_MAX 100

// dumb funcs to test...won't agree at -0 but don't care
float func0(float x) { return x>=0 ? 1.f : -1.f; }
float func1(float x) { return copysign(1.f,x); }

//-------------------------------

#if !defined(_MSC_VER)
#define FP32_MIN     0x1.0p-126f       // smallest normal
#define FP32_MIN_DN  0x1.0p-149f       // smallest denormal
#else
#include <float.h>
#define FP32_MIN     FLT_MIN
#define FP32_MIN_DN  FLT_TRUE_MIN
#endif

uint32_t ecount;

float error_table[ERROR_MAX];

void reset_test()
{
  ecount = 0;
}

void test_via_integer()
{
  union { float f; uint32_t i; } cnt;
  cnt.f = 0;

  do {
    float f0 = func0(cnt.f);
    float f1 = func1(cnt.f);
    cnt.i++;

    if (f0 == f1) continue;

    error_table[ecount++] = cnt.f;
    if (ecount == ERROR_MAX) return;

  } while(cnt.i != 0x7f800000);
}

void test_via_ulp()
{
  float x   = 0;
  float ulp = FP32_MIN_DN;
  float e   = 2.f*FP32_MIN;
  float ddf;

  do {
    while(x != e) {
      float f0 = func0(x);
      float f1 = func1(x);

      x += ulp;

      if (f0 == f1) continue;

      error_table[ecount++] = x;
      if (ecount == ERROR_MAX) return;
    }

    ulp += ulp;
    e   += e;

  } while(x <= 0x1.0p127f);
}


void run_tests()
{
  uint64_t s0,s1;

  printf("integer walk\n");
  for(uint32_t t=0; t<10; t++) {
    s0 = _rdtsc();
    test_via_integer();
    s1 = _rdtsc();
    printf("%f\n",1.0*(s1-s0));
  }

  printf("ulp walk\n");
  for(uint32_t t=0; t<10; t++) {
    s0 = _rdtsc();
    test_via_ulp();
    s1 = _rdtsc();
    printf("%f\n",1.0*(s1-s0));
  }
}


// walk through all the positive ordered binary32, both by integer & ulp increment
void sanity_check()
{
  uint32_t i=0;
  union { float f; uint32_t i; } cnt;
  cnt.f = 0;

  float x   = 0;
  float ulp = FP32_MIN_DN;
  float e   = 2.f*FP32_MIN;
  float ddf;

  do {
    while(x != e) {
      i++;
      cnt.i++;
      x += ulp;

      if (x == cnt.f) continue;

      printf("broken --- %08x %a %a | %a\n", i, cnt.f, x,ulp);
      return;
    }

    ulp += ulp;
    e   += e;

  } while(x <= 0x1.0p127f);
}


void main(int argc, char** argv)
{
  //sanity_check();
  run_tests();
}