#include <gmp.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <time.h>

typedef uint32_t int_equiv_t;
typedef int_equiv_t fb_t;
typedef float flt_equiv_t;

/* shift 0b111... right by n, then left by n, so you get 0b1111000...,
   then invert that. */
#define N_BITS(n) ~((~(uint32_t)0 >> (n)) << (n))

#define TEST_LIMIT (1L << (8 * sizeof(fb_t)))

#define BAR_BITS 5
#define ND_BITS 26

#define SIGN_MASK (1 << (8 * sizeof(fb_t) - 1))
#define BAR_MASK  (N_BITS(BAR_BITS) << ND_BITS)
#define BAR_SHIFT (ND_BITS)
#define ND_MASK   (~(SIGN_MASK | BAR_MASK))

int64_t
fb_num(fb_t v)
{
        int32_t s = (v & SIGN_MASK) ? -1 : 1;
        uint32_t b = (v & BAR_MASK) >> BAR_SHIFT;
        uint32_t nd = (v & ND_MASK);
        return s * (int32_t) (nd >> b);
}

uint64_t
fb_den(fb_t v)
{
        uint32_t b = (v & BAR_MASK) >> BAR_SHIFT;
        uint32_t nd = (v & ND_MASK);
        return (nd & N_BITS(b)) | (1 << b);
}

void
fb_print(fb_t v) {
        printf("%u: %ld/%lu\n", v, fb_num(v), fb_den(v));
}

void
fb_to_gmp(mpq_t rop, fb_t fb)
{
        mpq_set_si(rop, fb_num(fb), fb_den(fb));
}

void
flt_to_gmp(mpq_t rop, flt_equiv_t flt)
{
        /* this converts floats to doubles, but that conversion is exact
           since doubles are a superset of floats */
        mpq_set_d(rop, flt);
}

#define BUCKET_COUNT 100000
#define BUCKET_SIZE_NUM 1
#define BUCKET_SIZE_DEN 20

mpq_t limits[BUCKET_COUNT];
uint64_t fb_hist[BUCKET_COUNT] = {0};
uint64_t flt_hist[BUCKET_COUNT] = {0};

int
main()
{
        mpq_t bsize, qv;

        mpq_inits(bsize, qv, NULL);
        mpq_set_ui(bsize, BUCKET_SIZE_NUM, BUCKET_SIZE_DEN);

        for (uint64_t i = 0; i < BUCKET_COUNT; i++) {
                mpq_init(limits[i]);
                mpq_set_ui(limits[i], i, 1);
                mpq_mul(limits[i], limits[i], bsize);
        }

        clock_t start = clock();
        for (uint64_t lv = 0; lv < TEST_LIMIT; lv++) {
                int_equiv_t v = (int_equiv_t) lv;
                if (v % 1000000 == 0) {
                        float complete = (float)v / TEST_LIMIT;
                        double elapsed_sec = (double) (clock() - start) / CLOCKS_PER_SEC;
                        double est_rem = elapsed_sec / complete - elapsed_sec;
                        int s = ceil(est_rem);
                        int m = s / 60;
                        s -= 60 * m;
                        int h = m / 60;
                        m -= 60 * h;
                        fprintf(stderr, "%.1f%% (%u/%ld) %02d:%02d:%02d remaining\n",
                                100. * complete, v, TEST_LIMIT, h, m, s);
                }

                fb_t fb = (fb_t) v;
                fb_to_gmp(qv, fb);
                if (mpq_sgn(qv) >= 0) {
                        mpq_div(qv, qv, bsize);
                        double bf = mpq_get_d(qv);
                        /* do this before we convert to int to avoid potential overflow */
                        if (bf >= BUCKET_COUNT) {
                                bf = BUCKET_COUNT - 1;
                        }
                        int b = floor(bf);
                        fb_hist[b]++;
                }

                flt_equiv_t flt = *((flt_equiv_t*)&v);
                int class = fpclassify(flt);
                if (__builtin_expect(class != FP_ZERO && class != FP_NORMAL, 0)) {
                        continue;
                }
                flt_to_gmp(qv, flt);
                if (mpq_sgn(qv) >= 0) {
                        mpq_div(qv, qv, bsize);
                        double bf = mpq_get_d(qv);
                        if (bf >= BUCKET_COUNT) {
                                bf = BUCKET_COUNT - 1;
                        }
                        int b = floor(bf);
                        flt_hist[b]++;
                }
        }

        printf("lim,fb,flt\n");
        for (int i = 0; i < BUCKET_COUNT; i++) {
                char *limstr = mpq_get_str(NULL, 10, limits[i]);
                printf("%s,%ld,%ld\n", limstr, fb_hist[i], flt_hist[i]);
        }
}