/* ibmbench - time IBM System/360 float -> IEEE conversion on obsidian (Claude Opus 5.5, 2026-10-17).
 * Converts the samples of a raw SEG-Y file in memory, repeatedly; prints the best rate.
 * Usage: ibmbench file.sgy [reps]   (ns fixed at 3000, 240-byte trace headers, 3600-byte reel header) */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

static double now(void) { struct timeval t; gettimeofday(&t, 0); return t.tv_sec + t.tv_usec * 1e-6; }

/* the straightforward way: sign, base-16 exponent, 24-bit fraction */
static float ibm_ldexp(unsigned int w)
{
	int e = (int)((w >> 24) & 0x7f) - 64;
	double f = (double)(w & 0x00ffffff) / 16777216.0;
	double v = ldexp(f, 4 * e);
	return (float)((w & 0x80000000u) ? -v : v);
}

/* bit manipulation: normalise the hex fraction into an IEEE single (as SU's conv does, roughly) */
static float ibm_bits(unsigned int w)
{
	unsigned int fr = w & 0x00ffffff, s = w & 0x80000000u;
	int e = (int)((w >> 22) & 0x1fc) - 256 + 127 - 1;   /* 4*(exp-64), IEEE bias, minus 1 */
	union { unsigned int u; float f; } r;
	if (fr == 0) return 0.0f;
	while (!(fr & 0x00800000)) { fr <<= 1; e--; }
	if (e <= 0) return 0.0f;
	if (e >= 255) e = 254;
	r.u = s | ((unsigned int)e << 23) | (fr & 0x007fffff);
	return r.f;
}

int main(int argc, char **argv)
{
	FILE *fp; long n, i, ntr, j; unsigned int *buf, *w; float *out;
	int reps = argc > 2 ? atoi(argv[2]) : 7, r, k;
	double t0, best[2] = {1e9, 1e9}, chk = 0;
	if (argc < 2 || !(fp = fopen(argv[1], "rb"))) { fprintf(stderr, "usage\n"); return 1; }
	fseek(fp, 0, SEEK_END); n = ftell(fp); fseek(fp, 3600, SEEK_SET);
	ntr = (n - 3600) / (240 + 3000 * 4);
	buf = malloc((size_t)ntr * 3000 * 4); out = malloc((size_t)ntr * 3000 * 4);
	w = malloc(240 + 3000 * 4);
	for (i = 0; i < ntr; i++) {           /* keep samples only (AIX is big-endian: no swap) */
		fread(w, 1, 240 + 3000 * 4, fp);
		for (j = 0; j < 3000; j++) buf[i * 3000 + j] = w[60 + j];
	}
	n = ntr * 3000;
	for (r = 0; r < reps; r++) for (k = 0; k < 2; k++) {
		t0 = now();
		if (k == 0) for (i = 0; i < n; i++) out[i] = ibm_ldexp(buf[i]);
		else        for (i = 0; i < n; i++) out[i] = ibm_bits(buf[i]);
		t0 = now() - t0; if (t0 < best[k]) best[k] = t0;
		chk += out[n / 2];
	}
	for (i = 0; i < 1000; i++) if (ibm_ldexp(buf[i * 7]) != ibm_bits(buf[i * 7]) && buf[i * 7] & 0xffffff)
		{ fprintf(stderr, "mismatch at %ld: %g %g\n", i * 7, ibm_ldexp(buf[i * 7]), ibm_bits(buf[i * 7])); break; }
	printf("%ld samples; ldexp: %.3f s = %.2f M/s; bits: %.3f s = %.2f M/s (chk %g)\n",
	       n, best[0], n / best[0] / 1e6, best[1], n / best[1] / 1e6, chk);
	return 0;
}