Skip to content

Instantly share code, notes, and snippets.

@sffc
Last active February 21, 2017 20:51
Show Gist options
  • Save sffc/cc76182af46972bcf4a4104df6e45f1b to your computer and use it in GitHub Desktop.
Save sffc/cc76182af46972bcf4a4104df6e45f1b to your computer and use it in GitHub Desktop.
Comparison between "braindead" double-to-ascii and sprintf
#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
*/
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment