Comparison between "braindead" double-to-ascii and sprintf

braindead.c

#include "stdio.h"
#include "stdint.h"
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define LOG2_10 3.32192809488736

typedef enum { false, true } bool;

typedef union {
	double d;
	uint64_t l;
} ieee_t;

typedef struct {
	uint64_t bcd;
	int power;
	bool exact;
	ieee_t original;
} bcd_t;

ieee_t random_ieee() {
	ieee_t ieee;
	ieee.l = 
		(((uint64_t) rand() <<  0) & 0x000000000000FFFFL) | 
		(((uint64_t) rand() << 16) & 0x00000000FFFF0000L) | 
		(((uint64_t) rand() << 32) & 0x0000FFFF00000000L) |
		(((uint64_t) rand() << 48) & 0xFFFF000000000000L);
	return ieee;
}

bcd_t dbltobcd_sprintf(ieee_t input) {
	char rep[25];
	sprintf(rep, "%+1.*ex", DBL_DIG, input.d);
	uint64_t bcd = 0;
	bcd |= (((uint64_t)(rep[1]-'0')) << 60);
	for (int i=0; i<15; i++) {
		bcd |= (((uint64_t)(rep[i+3]-'0')) << ((14-i)*4));
	}
	int power = 0;
	for (int i=20; rep[i] != 'x'; i++) {
		power = power*10 + (rep[i]-'0');
	}
	if (rep[19] == '-') power *= -1;
	power -= 15;

	bcd_t result;
	result.bcd = bcd;
	result.power = power;
	result.exact = true;
	result.original = input;
	return result;
}

bcd_t dbltobcd_braindead(ieee_t input) {
	double n = input.d;
	if (n < 0) n *= -1;

	int exponent = ((input.l & 0x7ff0000000000000L) >> 52) - 0x3ff;
	int approximate_power = ((double)exponent) / LOG2_10;
	int delta = 15 - approximate_power;

	if (delta >= 0) {
		int i = delta;
		// 1e22 is the largest exactly representable power of ten
		for (; i >= 22; i -= 22) n *= 1e22;
		for (; i >= 9; i -= 9) n *= 1000000000;
		for (; i >= 3; i -= 3) n *= 1000;
		for (; i >= 1; i -= 1) n *= 10;
	} else {
		int i = delta;
		// 1e22 is the largest exactly representable power of ten
		for (; i <= -22; i += 22) n /= 1e22;
		for (; i <= -9; i += 9) n /= 1000000000;
		for (; i <= -3; i += 3) n /= 1000;
		for (; i <= -1; i += 1) n /= 10;
	}

	uint64_t bcd = 0;
	uint64_t rec = llround(n);
	int num_digits = 0;
	while (rec != 0) {
		uint64_t digit = rec % 10;
		rec /= 10;
		bcd = (bcd >> 4) | (digit << 60);
		num_digits++;
	}
	int power = num_digits - delta - 16;

	bcd_t result;
	result.bcd = bcd;
	result.power = power;
	result.exact = false;
	result.original = input;
	return result;
}

int bcdcmp(bcd_t bcd1, bcd_t bcd2) {
	if (bcd1.power > bcd2.power) return 1;
	else if (bcd1.power < bcd2.power) return -1;
	else if (bcd1.bcd > bcd2.bcd) return 1;
	else if (bcd1.bcd < bcd2.bcd) return -1;
	else return 0;
}

uint64_t bcd2long(bcd_t bcd) {
	uint64_t result = 0L;
	for (int i=15; i>=0; i--) {
		result = result*10 + ((bcd.bcd >> (4*i)) & 0xf);
	}
	return result;
}

uint64_t bcd_difference(bcd_t bcd1, bcd_t bcd2) {
	// Make bcd1 smaller than bcd2
	if (bcdcmp(bcd1, bcd2) > 0) {
		bcd_t temp = bcd2;
		bcd2 = bcd1;
		bcd1 = temp;
	}

	uint64_t l1 = bcd2long(bcd1);
	uint64_t l2 = bcd2long(bcd2);
	for (int d=bcd2.power-bcd1.power; d>0; d++) {
		l2 *= 10;
	}

	return l2 - l1;
}

uint64_t bcd_unsafe_round(bcd_t bcd, int num_digits) {
	int first_digit = (bcd.bcd>>((16-num_digits-1)*4)) & 0xf;
	bcd.bcd >>= (16-num_digits)*4;
	uint64_t result = bcd2long(bcd);
	if (first_digit >= 5) result++;
	return result;
}

uint64_t bcd_safe_round(bcd_t bcd, int num_digits, int* num_fallbacks) {
	bool requires_fallback = true;
	int prev_digit = -1;
	for (int i=16-num_digits-1; i>=1; i--) {
		int digit = (bcd.bcd>>(i*4)) & 0xf;
		if (prev_digit == -1) {
			requires_fallback = (digit == 4 || digit == 5);
		} else if (prev_digit == 0) {
			requires_fallback = (digit == 0);
		} else if (prev_digit == 4) {
			requires_fallback = (digit == 9);
		} else if (prev_digit == 5) {
			requires_fallback = (digit == 0);
		} else if (prev_digit == 9) {
			requires_fallback = (digit == 9);
		}
		prev_digit = digit;
		if (!requires_fallback) break;
	}

	if (requires_fallback) {
		// printf("---fallback\n");
		// printf("original: %+e\n", bcd.original.d);
		// printf("long: %lld\n", bcd2long(bcd));
		(*num_fallbacks)++;
		bcd_t fallback = dbltobcd_sprintf(bcd.original);
		return bcd_unsafe_round(fallback, num_digits);
	} else {
		return bcd_unsafe_round(bcd, num_digits);
	}
}

#define DIFFS_LENGTH 26
int main() {
	printf("Note: RAND_MAX=%d\n", RAND_MAX);
	srand(time(NULL));
	int num_fallbacks = 0;
	int diffs[1000];
	for (int i=0; i<DIFFS_LENGTH; i++) diffs[i] = 0;
	for (int i=0; i<1000000; i++) {
		ieee_t ieee = random_ieee();
		if (ieee.d != ieee.d) continue; // NaN
		bcd_t bcd1 = dbltobcd_sprintf(ieee);
		bcd_t bcd2 = dbltobcd_braindead(ieee);
		uint64_t diff = bcd_difference(bcd1, bcd2);
		if (diff >= DIFFS_LENGTH - 1) {
			// Pool out-of-bounds diffs all into the highest bin
			diff = DIFFS_LENGTH - 1;
		}
		diffs[diff]++;
		uint64_t rnd1 = bcd_unsafe_round(bcd1, 12);
		// uint64_t rnd2 = bcd_unsafe_round(bcd2, 12);
		uint64_t rnd2 = bcd_safe_round(bcd2, 12, &num_fallbacks);
		if (rnd1 != rnd2) {
			printf("---\n");
			printf("ieee value: %+e\n", ieee.d);
			printf("ieee bytes: %llx\n", ieee.l);
			printf("sprintf: %llxE%d\n", bcd1.bcd, bcd1.power);
			printf("brndead: %llxE%d\n", bcd2.bcd, bcd2.power);
			printf("long1: %lld\n", bcd2long(bcd1));
			printf("long2: %lld\n", bcd2long(bcd2));
			printf("cmp: %d\n", bcdcmp(bcd1, bcd2));
			printf("differe: %lld\n", bcd_difference(bcd1, bcd2));
			printf("rnd1: %lld\n", rnd1);
			printf("rnd2: %lld\n", rnd2);
		}
	}
	for (int i=0; i<DIFFS_LENGTH; i++) {
		printf("%d: %d\n", i, diffs[i]);
	}
	printf("Total rounding fallbacks: %d\n", num_fallbacks);
	return 0;
}

/*
Difficult examples:

ieee value: -5.074791e-312
ieee bytes: 800000ef26dbb168
sprintf: 5074790912492772E-327
brndead: 5074790912500000E-327
long1: 5074790912492772
long2: 5074790912500000
cmp: -1
differe: 7228
rnd1: 507479091249
rnd2: 507479091250
---
ieee value: +8.360253e-312
ieee bytes: 189fb0c787b
sprintf: 8360253001975257E-327
brndead: 8360253002000000E-327
long1: 8360253001975257
long2: 8360253002000000
cmp: -1
differe: 24743
rnd1: 836025300198
rnd2: 836025300200
---
ieee value: -9.083660e-274
ieee bytes: 873f732164c3d408
sprintf: 9083659622814999E-289
brndead: 9083659622815010E-289
long1: 9083659622814999
long2: 9083659622815010
cmp: -1
differe: 11
rnd1: 908365962281
rnd2: 908365962282
---
ieee value: +8.646301e-311
ieee bytes: fea9babdf5a
sprintf: 8646301219605000E-326
brndead: 8646301219600000E-326
long1: 8646301219605000
long2: 8646301219600000
cmp: 1
differe: 5000
rnd1: 864630121961
rnd2: 864630121960
*/