schubfach.hpp

Go to the documentation of this file.

 
#ifndef __DCONV_SCHUBFACH_HPP__
#define __DCONV_SCHUBFACH_HPP__
 
// dconv.
#include <dconv/schubfachpow.hpp>
#include <dconv/format.hpp>
 
// C.
#include <cstring>
#include <cstdint>
 
namespace dconv
{
    namespace schubfach
    {
        constexpr int floorLog10Pow2 (int q, bool asymmetric) noexcept
        {
            return (q * 315653 - (asymmetric ? 131008 : 0)) >> 20;
        }
 
        constexpr int floorLog2Pow10 (int k) noexcept
        {
            return (k * 3483294) >> 20;
        }
 
        inline uint64_t roundToOdd (const details::Power& g, uint64_t cp) noexcept
        {
            uint64_t middle, unused;
            details::umul128 (g.lo, cp, middle, unused);
 
            uint64_t high, low;
            details::umul128 (g.hi, cp, high, low);
 
            const uint64_t z = low + middle;
            high += (z < low);
            return high | (z != 0);
        }
 
        inline uint64_t shortest (double value, int& k) noexcept
        {
            uint64_t bits;
            memcpy (&bits, &value, sizeof (double));
 
            const uint64_t fraction = bits & 0x000FFFFFFFFFFFFF;
            const int biased = static_cast<int> (bits >> 52);
            const uint64_t c = (biased != 0) ? (fraction | 0x0010000000000000) : fraction;
            const int q = (biased != 0) ? biased - 1075 : -1074;
 
            if ((static_cast<unsigned> (-q) <= 52) && ((c & ((static_cast<uint64_t> (1) << -q) - 1)) == 0))
            {
                k = 0;
                return c >> -q;
            }
 
            const bool asymmetric = (fraction == 0) && (biased > 1);
            const uint64_t out = c & 1;
 
            const uint64_t cb = c << 2;
            const uint64_t cbl = asymmetric ? cb - 1 : cb - 2;
            const uint64_t cbr = cb + 2;
 
            k = floorLog10Pow2 (q, asymmetric);
            const int h = q + floorLog2Pow10 (-k) + 3;
            const details::Power& g = details::schubfachpow[292 - k];
 
            const uint64_t vb = roundToOdd (g, cb << h);
            const uint64_t vbl = roundToOdd (g, cbl << h);
            const uint64_t vbr = roundToOdd (g, cbr << h);
 
            const uint64_t s = vb >> 2;
            const uint64_t sp10 = (s / 10) * 10;
            const uint64_t tp10 = sp10 + 10;
            const bool upin = vbl + out <= (sp10 << 2);
            const bool wpin = (tp10 << 2) + out <= vbr;
            if (upin != wpin)
            {
                uint64_t digits = upin ? s / 10 : s / 10 + 1;
                ++k;
                if (digits < 1000000000000000)
                {
                    digits *= 10;
                    --k;
                }
                return digits;
            }
 
            const uint64_t t = s + 1;
            const bool uin = vbl + out <= (s << 2);
            const bool win = (t << 2) + out <= vbr;
            if (uin != win)
            {
                return uin ? s : t;
            }
 
            const int64_t cmp = static_cast<int64_t> (vb - ((s + t) << 1));
            return ((cmp < 0) || ((cmp == 0) && ((s & 1) == 0))) ? s : t;
        }
 
        inline void generate (char* buffer, int& length, int& k, double value) noexcept
        {
            const uint64_t digits = shortest (value, k);
            length = details::countDigits (digits);
            const int zeros = details::writeDigits (buffer, digits, length);
            length -= zeros;
            k += zeros;
        }
 
        inline char* dtoa (char* buffer, double value) noexcept
        {
            uint64_t bits;
            memcpy (&bits, &value, sizeof (double));
            bool is_negative = (bits >> 63) != 0;
 
            *buffer = '-';
            buffer += is_negative;
            value = is_negative ? -value : value;
 
            if (value == 0.0)
            {
                memcpy (buffer, "0.0", 3);
                return buffer + 3;
            }
 
            int k = 0;
            const uint64_t digits = shortest (value, k);
            return details::format (buffer, digits, details::countDigits (digits), k);
        }
    }
}
 
#endif