Class diyfp

Source/external/json.hpp:16840–16956

}

struct diyfp // f * 2^e
{
	static constexpr int kPrecision = 64; // = q

	std::uint64_t f = 0;
	int e = 0;

	constexpr diyfp(std::uint64_t f_, int e_) noexcept : f(f_), e(e_) {}

	/*!
	@brief returns x - y
	@pre x.e == y.e and x.f >= y.f
	*/
	static diyfp sub(const diyfp& x, const diyfp& y) noexcept
	{
		JSON_ASSERT(x.e == y.e);
		JSON_ASSERT(x.f >= y.f);

		return {x.f - y.f, x.e};
	}

	/*!
	@brief returns x * y
	@note The result is rounded. (Only the upper q bits are returned.)
	*/
	static diyfp mul(const diyfp& x, const diyfp& y) noexcept
	{
		static_assert(kPrecision == 64, "internal error");

		// Computes:
		//  f = round((x.f * y.f) / 2^q)
		//  e = x.e + y.e + q

		// Emulate the 64-bit * 64-bit multiplication:
		//
		// p = u * v
		//   = (u_lo + 2^32 u_hi) (v_lo + 2^32 v_hi)
		//   = (u_lo v_lo		 ) + 2^32 ((u_lo v_hi		 ) + (u_hi v_lo		 )) + 2^64 (u_hi v_hi		 )
		//   = (p0				) + 2^32 ((p1				) + (p2				)) + 2^64 (p3				)
		//   = (p0_lo + 2^32 p0_hi) + 2^32 ((p1_lo + 2^32 p1_hi) + (p2_lo + 2^32 p2_hi)) + 2^64 (p3				)
		//   = (p0_lo			 ) + 2^32 (p0_hi + p1_lo + p2_lo					  ) + 2^64 (p1_hi + p2_hi + p3)
		//   = (p0_lo			 ) + 2^32 (Q										  ) + 2^64 (H				 )
		//   = (p0_lo			 ) + 2^32 (Q_lo + 2^32 Q_hi						   ) + 2^64 (H				 )
		//
		// (Since Q might be larger than 2^32 - 1)
		//
		//   = (p0_lo + 2^32 Q_lo) + 2^64 (Q_hi + H)
		//
		// (Q_hi + H does not overflow a 64-bit int)
		//
		//   = p_lo + 2^64 p_hi

		const std::uint64_t u_lo = x.f & 0xFFFFFFFFu;
		const std::uint64_t u_hi = x.f >> 32u;
		const std::uint64_t v_lo = y.f & 0xFFFFFFFFu;
		const std::uint64_t v_hi = y.f >> 32u;

		const std::uint64_t p0 = u_lo * v_lo;