#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]);
}
}