/******************************************************************************
 * The MIT License (MIT)
 * 
 * Copyright (c) 2014 Baldur Karlsson
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 * 
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 ******************************************************************************/

#include "common/common.h"

#include <math.h>

///////////////////////////////////////////////////////////////////////////
// Grisu2 implementation (slightly simpler than Grisu3) for converting
// doubles to strings
//
// Sources:
// Based on Florian Loitsch 2010 "Printing Floating-Point Numbers Quickly
//                                and Accurately with Integers"
//     http://florian.loitsch.com/publications/dtoa-pldi2010.pdf
//     https://github.com/floitsch/double-conversion (BSD licensed)
//
// Also implementations by Milo Yip and night-shift used as reference
//     https://github.com/miloyip/dtoa-benchmark
//     https://github.com/night-shift/fpconv

struct diy_fp
{
	diy_fp() {}
	diy_fp(uint64_t mant, int exponent) : mantissa(mant), exp(exponent) {}
	uint64_t mantissa;
	int exp;

	// q in the paper, bits in the mantissa of the fixed point
	// approximation
	static const int bitsq = 64;
};

// subtract from Florian paper
diy_fp operator -(const diy_fp &x, const diy_fp &y)
{
	// assume same exponent
	return diy_fp(x.mantissa-y.mantissa, x.exp);
}

// multiply from Florian paper
diy_fp operator *(const diy_fp &x, const diy_fp &y)
{
	// _a = upper 32 bits, _b = lower 32 bits
	uint64_t xa = x.mantissa >> 32, xb = x.mantissa & 0xFFFFFFFF;
	uint64_t ya = y.mantissa >> 32, yb = y.mantissa & 0xFFFFFFFF;

	// perform each pair of multiplies
	uint64_t upper = xa*ya;
	uint64_t lower = xb*yb;
	uint64_t cross1 = xb*ya;
	uint64_t cross2 = xa*yb;

	uint64_t tmp = (lower>>32) + (cross1&0xFFFFFFFF) + (cross2&0xFFFFFFFF);
	tmp += 1U << 31; // Round up

	// note - exponent is no longer normalised
	return diy_fp(upper + (cross1>>32) + (cross2>>32) + (tmp>>32), x.exp + y.exp + 64);
}

static diy_fp pow10cache[] = {
	diy_fp(18054884314459144840U, -1220),
	diy_fp(13451937075301367670U, -1193),
	diy_fp(10022474136428063862U, -1166),
	diy_fp(14934650266808366570U, -1140),
	diy_fp(11127181549972568877U, -1113),
	diy_fp(16580792590934885855U, -1087),
	diy_fp(12353653155963782858U, -1060),
	diy_fp(18408377700990114895U, -1034),
	diy_fp(13715310171984221708U, -1007),
	diy_fp(10218702384817765436U, -980),
	diy_fp(15227053142812498563U, -954),
	diy_fp(11345038669416679861U, -927),
	diy_fp(16905424996341287883U, -901),
	diy_fp(12595523146049147757U, -874),
	diy_fp(9384396036005875287U,  -847),
	diy_fp(13983839803942852151U, -821),
	diy_fp(10418772551374772303U, -794),
	diy_fp(15525180923007089351U, -768),
	diy_fp(11567161174868858868U, -741),
	diy_fp(17236413322193710309U, -715),
	diy_fp(12842128665889583758U, -688),
	diy_fp(9568131466127621947U,  -661),
	diy_fp(14257626930069360058U, -635),
	diy_fp(10622759856335341974U, -608),
	diy_fp(15829145694278690180U, -582),
	diy_fp(11793632577567316726U, -555),
	diy_fp(17573882009934360870U, -529),
	diy_fp(13093562431584567480U, -502),
	diy_fp(9755464219737475723U,  -475),
	diy_fp(14536774485912137811U, -449),
	diy_fp(10830740992659433045U, -422),
	diy_fp(16139061738043178685U, -396),
	diy_fp(12024538023802026127U, -369),
	diy_fp(17917957937422433684U, -343),
	diy_fp(13349918974505688015U, -316),
	diy_fp(9946464728195732843U,  -289),
	diy_fp(14821387422376473014U, -263),
	diy_fp(11042794154864902060U, -236),
	diy_fp(16455045573212060422U, -210),
	diy_fp(12259964326927110867U, -183),
	diy_fp(18268770466636286478U, -157),
	diy_fp(13611294676837538539U, -130),
	diy_fp(10141204801825835212U, -103),
	diy_fp(15111572745182864684U, -77),
	diy_fp(11258999068426240000U, -50),
	diy_fp(16777216000000000000U, -24),
	diy_fp(12500000000000000000U,   3),
	diy_fp(9313225746154785156U,   30),
	diy_fp(13877787807814456755U,  56),
	diy_fp(10339757656912845936U,  83),
	diy_fp(15407439555097886824U, 109),
	diy_fp(11479437019748901445U, 136),
	diy_fp(17105694144590052135U, 162),
	diy_fp(12744735289059618216U, 189),
	diy_fp(9495567745759798747U,  216),
	diy_fp(14149498560666738074U, 242),
	diy_fp(10542197943230523224U, 269),
	diy_fp(15709099088952724970U, 295),
	diy_fp(11704190886730495818U, 322),
	diy_fp(17440603504673385349U, 348),
	diy_fp(12994262207056124023U, 375),
	diy_fp(9681479787123295682U,  402),
	diy_fp(14426529090290212157U, 428),
	diy_fp(10748601772107342003U, 455),
	diy_fp(16016664761464807395U, 481),
	diy_fp(11933345169920330789U, 508),
	diy_fp(17782069995880619868U, 534),
	diy_fp(13248674568444952270U, 561),
	diy_fp(9871031767461413346U,  588),
	diy_fp(14708983551653345445U, 614),
	diy_fp(10959046745042015199U, 641),
	diy_fp(16330252207878254650U, 667),
	diy_fp(12166986024289022870U, 694),
	diy_fp(18130221999122236476U, 720),
	diy_fp(13508068024458167312U, 747),
	diy_fp(10064294952495520794U, 774),
	diy_fp(14996968138956309548U, 800),
	diy_fp(11173611982879273257U, 827),
	diy_fp(16649979327439178909U, 853),
	diy_fp(12405201291620119593U, 880),
	diy_fp(9242595204427927429U,  907),
	diy_fp(13772540099066387757U, 933),
	diy_fp(10261342003245940623U, 960),
	diy_fp(15290591125556738113U, 986),
	diy_fp(11392378155556871081U, 1013),
	diy_fp(16975966327722178521U, 1039),
	diy_fp(12648080533535911531U, 1066),
};

static const int firstpow10 = -348; // first cached power of 10
static const int cachestep = 8; // power of 10 steps between cache items

diy_fp find_cachedpow10(int exp, int& kout)
{
	const double inv_log2_10 = 0.30102999566398114;
	const double alpha = -60.0;

	// k calculation from the paper ceil[ (alpha - exp + q - 1) * 1/log2(10) ]
	// exponent is shifted by #bits
	int k = (int)ceil( (alpha - double(exp + diy_fp::bitsq) + diy_fp::bitsq - 1) * inv_log2_10 );

	// determine index in above array
	int idx = (-firstpow10 + k - 1) / cachestep + 1;

	// output the decimal power that corresponds to this k
	kout = (firstpow10 + idx * cachestep);

	return pow10cache[idx];
}

static int gen_digits(const diy_fp &lower, const diy_fp &upper, char *digits, int &kout)
{
	diy_fp delta = upper-lower;

	// generate 1.0 to the desired exponent so we can split integer from decimal part
	diy_fp one(uint64_t(1) << -upper.exp, upper.exp);

	// mask off integer and decimal parts
	uint64_t intpart = upper.mantissa >> -one.exp;
	uint64_t decpart = upper.mantissa & (one.mantissa - 1);

	// len is current number of digits produced
	int len = 0;
	// kappa is an exponent shift, to account for if we don't produce exactly the number
	// of digits to reach the decimal place, and there should be extra 0s beyond the produced
	// digits. (or negative if there should be preceeding 0s)
	int kappa = 10;
	uint32_t div = 1000000000; // highest possible pow10 in 32bits = 10^9

	// handle integer component before decimal separator
	while(kappa > 0)
	{
		// get digit at current power of ten
		uint64_t digit = intpart / div;

		// don't include preceeding 0 digits (so either include if
		// digit is non-0, or if we've started including digits ie.
		// len > 0)
		if(digit || len) digits[len++] = '0' + char(digit);

		// remove this pow10 from the int for future iterations
		intpart %= div; kappa--; div /= 10;

		// this is our termination condition, when we've produced the number.
		// delta is the difference between upper and lower, and the left side
		// is the current remainder after the currently generated digits have
		// been removed. If that is small enough that we've produced the number,
		// exit and increment kout to account for the extra exponential
		if( (intpart << -one.exp) + decpart <= delta.mantissa)
		{
			kout += kappa;
			return len;
		}
	}

	// note, after this part if we're still here, intpart is 0 as we've
	// masked off all digits, so only decpart remains.
	// Kappa has also reached 0, beyond here it decrements below 0

	// handle decimal portion after separator
	do
	{
		decpart *= 10;
		uint64_t digit = decpart >> -one.exp;

		// don't include preceeding 0s (see above - note if we've produced
		// any integer digits at all, len will be > 0)
		if(digit || len) digits[len++] = '0' + char(digit);

		// remove this pow10 from the decimal part
		decpart &= (one.mantissa-1); kappa--; delta.mantissa *= 10;

		// stop looping when decpart is lower than delta (see above for termination condition)
	} while(decpart > delta.mantissa);

	kout += kappa;

	return len;
}

int grisu2(uint64_t mantissa, int exponent, char digits[18], int &kout)
{
	// the IEEE format implicitly has a hidden 1 bit above the mantissa for all normalised
	// numbers
	const uint64_t hiddenbit = 0x0010000000000000;

	// exponent is shifted by a further 52 because input exponent is assuming mantissa
	// is 1.2345678...e exp (fraction)
	// but grisu2 treats number as
	//    12345678...e exp-52 (whole number)
	diy_fp w = diy_fp(mantissa|hiddenbit, exponent - 52);
	if(exponent == -1023) w.exp = 1 - (1023 + 52); // subnormal exponent

	// we know the w input comes from a double, so is only using the lower 52 bits at
	// most. We can safely multiply by 2 (and cancel by lowering exponent to match), then
	// add 1 to get the upper value
	diy_fp upper;
	upper.mantissa = (w.mantissa<<1) + 1;
	upper.exp = w.exp-1;

	diy_fp lower;
	if(mantissa == 0)
	{
		// if mantissa is 0 we are going to underflow, so shift by 2
		// to maintain precision/normalised value
		lower.mantissa = (w.mantissa<<2) - 1;
		lower.exp = w.exp-2;
	}
	else
	{
		lower.mantissa = (w.mantissa<<1) - 1;
		lower.exp = w.exp-1;
	}

	// normalise upper - shift the mantissa until the top bit is set, and decrement the
	// exponent each time to keep the number the same (1.2e5 and 12e4 are equivalent representations)
	while((upper.mantissa & 0x8000000000000000) == 0)
	{
		upper.mantissa <<= 1;
		upper.exp--;
	}
	// lower needs to be the same exponent as upper so we can calculate delta by upper-lower, so
	// it is not normalised, but shifted to upper's exponent
	lower.mantissa <<= (lower.exp - upper.exp);
	lower.exp = upper.exp;

	int k = 0;
	diy_fp ck = find_cachedpow10(upper.exp, k);

	lower = lower * ck;
	upper = upper * ck;

	// squeeze the range in by 1 ULP
	lower.mantissa++; upper.mantissa--;

	// set our initial exponent. This will be shifted
	// if we have any preceeding or trailing 0s to get the
	// final exponent
	kout = -k;

	return gen_digits(lower, upper, digits, kout);
}
