stod.c

#include <stdint.h>
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <stdio.h>
#include <string.h>

#include <time.h>

#if defined(_MSC_VER)

#include <intrin.h>

#define stod_noinline __declspec(noinline)
#define stod_forceinline __forceinline

static uint32_t nbits(uint64_t v)
{
	if (v == 0) return 0;
	unsigned long index = 0;
	_BitScanReverse64(&index, (unsigned long long)v);
	index++;
	return (uint32_t)index;
}

#define divhi128(hi, div, p_rem) _udiv128(hi, 0, div, p_rem)

#else

#define stod_noinline __attribute__((noinline))
#define stod_forceinline inline __attribute__((always_inline))

static uint32_t nbits(uint64_t v)
{
	if (v == 0) return 0;
	return 64 - __builtin_clzll(v);
}

static stod_forceinline
uint64_t divhi128(uint64_t hi, uint64_t div, uint64_t *p_rem)
{
	unsigned __int128 hi128 = (unsigned __int128)hi << 64u;
	uint64_t a = (uint64_t)(hi128 / div);
	uint64_t b = (uint64_t)(hi128 % div);
	*p_rem = b;
	return a;
}

#endif

static const uint64_t pow10[] = {
	UINT64_C(1),
	UINT64_C(10),
	UINT64_C(100),
	UINT64_C(1000),
	UINT64_C(10000),
	UINT64_C(100000),
	UINT64_C(1000000),
	UINT64_C(10000000),
	UINT64_C(100000000),
	UINT64_C(1000000000),
	UINT64_C(10000000000),
	UINT64_C(100000000000),
	UINT64_C(1000000000000),
	UINT64_C(10000000000000),
	UINT64_C(100000000000000),
	UINT64_C(1000000000000000),
	UINT64_C(10000000000000000),
	UINT64_C(100000000000000000),
	UINT64_C(1000000000000000000),
	UINT64_C(10000000000000000000),
};

stod_noinline double parse_double(const char *str)
{
	uint64_t integer = 0;
	uint64_t decimals = 0;
	uint32_t n_decimals = 0;
	bool negative = false;

	if (*str == '-') {
		negative = true;
		str++;
	}
	while (*str && *str != '.') {
		integer = integer * 10 + (uint64_t)(*str++ - '0');
	}
	if (*str == '.') {
		str++;
		while (*str) {
			decimals = decimals * 10 + (uint64_t)(*str++ - '0');
			n_decimals++;
		}
	}

	if (!decimals) {
		return (negative ? -1.0 : 1.0) * (double)integer;
	}

	uint64_t divisor = pow10[n_decimals];

	uint64_t b_int = integer;
	uint32_t n_int = nbits(b_int);

	uint64_t rem_hi, rem_lo;
	uint64_t b_hi = divhi128(decimals, divisor, &rem_hi);
	uint32_t n_hi = nbits(b_hi);

	int32_t exponent;
	uint64_t mantissa;
	bool nonzero_tail;
	if (b_int) {
		mantissa = b_int << (64u - n_int) | (b_hi >> n_int);
		nonzero_tail = (b_hi << (64u - n_int) | rem_hi) != 0;
		exponent = (int32_t)n_int - 1;
	} else if (n_hi >= 54) {
		mantissa = b_hi << (64u - n_hi);
		nonzero_tail = rem_hi != 0;
		exponent = (int32_t)n_hi - 65;
	} else {
		uint64_t b_lo = divhi128(rem_hi, divisor, &rem_lo);
		mantissa = b_hi << (64u - n_hi) | (b_lo >> n_hi);
		nonzero_tail = (b_lo << (64u - n_hi) | rem_lo) != 0;
		exponent = (int32_t)n_hi - 65;
	}

	bool r_odd = mantissa & (1 << 11u);
	bool r_round = mantissa & (1 << 10u);
	bool r_tail = (mantissa & (1 << 10u)) - 1 != 0 || nonzero_tail;
	
	uint64_t round = (r_round && (r_odd || r_tail)) ? 1u : 0u;

	uint64_t bits
		= (uint64_t)negative << 63u
		| (uint64_t)(exponent + 1023) << 52u
		| ((mantissa >> 11u) & ~(UINT64_C(1) << 52u));
	bits += round;

	double result;
	memcpy(&result, &bits, 8);
	return result;
}

static stod_forceinline bool inc(char *buf, size_t len)
{
	char *s = buf + (len - 2);
	for (;;) {
		if (*s == '9') {
			*s = '0';
			if (s == buf) break;
			--s;
		} else if (*s == '.') {
			if (s == buf) break;
			--s;
		} else {
			++*s;
			break;
		}
	}
	return true;
}

stod_noinline void bench_a(double *dst, size_t count)
{
	char num[] = "1.000000000000001";

	for (size_t i = 0; i < count; i++) {
		dst[i] = parse_double(num);
		if (!inc(num, sizeof(num))) break;
	}
}

stod_noinline void bench_b(double *dst, size_t count)
{
	char num[] = "1.000000000000001";

	for (size_t i = 0; i < count; i++) {
		dst[i] = strtod(num, NULL);
		if (!inc(num, sizeof(num))) break;
	}
}

int main()
{
#if 0
	size_t bench_size = 64*1024*1024;
	double *result = (double*)malloc(bench_size * sizeof(double));
	memset(result, 0, bench_size);

	clock_t a_begin = clock();
	bench_a(result, bench_size);
	clock_t a_end = clock();

	printf("a: %.2fs\n", (double)(a_end - a_begin) / (double)CLOCKS_PER_SEC);

	clock_t b_begin = clock();
	bench_b(result, bench_size);
	clock_t b_end = clock();

	printf("b: %.2fs\n", (double)(b_end - b_begin) / (double)CLOCKS_PER_SEC);
#endif

	//char num[] = "0.200000000000001";
	char num[] = "0000.000000";
	// char num[] = "0.000000000014479";

	do {
		double a = parse_double(num);
		double b = strtod(num, NULL);
		if (a != b) {
			printf("%s -> %.20f, %.20f\n", num, a, b);
		}
	} while (inc(num, sizeof(num)));

	return 1;
}