Class diyfp

Source/Utils/json.hpp:14388–14504

            }

            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;