Rocksolid Light

Welcome to RetroBBS

mail  files  register  newsreader  groups  login

Message-ID:  

Uncompensated overtime? Just Say No.


devel / comp.lang.c++ / Re: from_chars vs my parse_double<>

SubjectAuthor
* from_chars vs my parse_double<>Bonita Montero
`- Re: from_chars vs my parse_double<>Bonita Montero

1
from_chars vs my parse_double<>

<udcrpm$31mqp$2@dont-email.me>

  copy mid

https://www.rocksolidbbs.com/devel/article-flat.php?id=1352&group=comp.lang.c%2B%2B#1352

  copy link   Newsgroups: comp.lang.c++
Path: i2pn2.org!i2pn.org!eternal-september.org!news.eternal-september.org!.POSTED!not-for-mail
From: Bonita.Montero@gmail.com (Bonita Montero)
Newsgroups: comp.lang.c++
Subject: from_chars vs my parse_double<>
Date: Thu, 7 Sep 2023 17:55:36 +0200
Organization: A noiseless patient Spider
Lines: 58
Message-ID: <udcrpm$31mqp$2@dont-email.me>
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 7bit
Injection-Date: Thu, 7 Sep 2023 15:55:34 -0000 (UTC)
Injection-Info: dont-email.me; posting-host="707347a35c765578889e4f3b790e82e5";
logging-data="3201881"; mail-complaints-to="abuse@eternal-september.org"; posting-account="U2FsdGVkX18YVzbG3Ayd9BiVDeQTrnoMOdz8+0sMo0o="
User-Agent: Mozilla Thunderbird
Cancel-Lock: sha1:Yg5izO8T0ynquY5+wudp52jJ7AA=
Content-Language: de-DE
 by: Bonita Montero - Thu, 7 Sep 2023 15:55 UTC

#include <iostream>
#include <type_traits>
#include <string>
#include <chrono>
#include "char_conv.h"

using namespace std;
using namespace chrono;

int main()
{ string strValue( "-3.14159266359e-300" );
else
auto compare = [&]<bool Precise, bool Std>( char const *prefix,
uint64_t reference, bool_constant<Precise>, bool_constant<Std> )
{
double value;
if constexpr( !Std )
parse_double<Precise>( strValue.cbegin(), strValue.cend(), value );
else
from_chars( strValue.data(), strValue.data() + strValue.size(), value );
uint64_t bin = bit_cast<uint64_t>( value );
cout << prefix << hexfloat << value;
if( unsigned diffBits = 64 - countl_zero( bin ^ reference ); reference )
cout << ": " << diffBits;
cout << endl;
return bit_cast<uint64_t>( value );
};
uint64_t
dummyReference = 0,
pdImprecise = compare( "imprecise: ", dummyReference,
false_type(), false_type() ),
pdPrecise = compare( "precise vs imprecise: ", pdImprecise,
true_type(), false_type() ),
fcVsImprecise = compare( "from_chars vs. ip: ", pdImprecise,
false_type(), true_type() ),
fcVsPrecise = compare( "from_chars vs. p: ", pdPrecise,
false_type(), true_type() );
}

I implemented sth. like from_chars myself. I wanted to have a maximum,
precision with that. For the given value above the highest different
bit of my precise solution vs. from_chars is bit 11, so 12 bits of
the mantissa are more exact.
If you take a naive approach and sum up the suffix-digits multiplied
by their 10 ^ N value from left to right the internediate values you
add become smaller the further you get right and at last won't partipate
in the mantissa. So I do the math from right to left by first putting
everything into a table (100 digits, if overflow everything is put into
a thread-local vector which isn't reallocated for the next call). And
if you get further in the 10 ^ N row there are N errors of each / 10
or * 10 operation. My function has a template parameter which leads to
a computation of the 10 ^ N value for each digit by my own pow10()
function, which feeds from a table for each bit of the N-exponent,
i.e. there are two tables which store the positive and negative
eponent's values (10 ^ (2 ^ N)).
With the more precise solution I get a result where the lower 12
mantissa bits are different for the above string value.

Re: from_chars vs my parse_double<>

<uervti$1tv29$1@dont-email.me>

  copy mid

https://www.rocksolidbbs.com/devel/article-flat.php?id=1868&group=comp.lang.c%2B%2B#1868

  copy link   Newsgroups: comp.lang.c++
Path: i2pn2.org!i2pn.org!eternal-september.org!news.eternal-september.org!.POSTED!not-for-mail
From: Bonita.Montero@gmail.com (Bonita Montero)
Newsgroups: comp.lang.c++
Subject: Re: from_chars vs my parse_double<>
Date: Mon, 25 Sep 2023 14:54:11 +0200
Organization: A noiseless patient Spider
Lines: 602
Message-ID: <uervti$1tv29$1@dont-email.me>
References: <udcrpm$31mqp$2@dont-email.me>
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 7bit
Injection-Date: Mon, 25 Sep 2023 12:54:10 -0000 (UTC)
Injection-Info: dont-email.me; posting-host="e2ba2e234a6dcf9f0afcd3d90854525d";
logging-data="2030665"; mail-complaints-to="abuse@eternal-september.org"; posting-account="U2FsdGVkX1+InVxQgUSCf+A5BNuipq4jCLG41oMZkjU="
User-Agent: Mozilla Thunderbird
Cancel-Lock: sha1:IwoHDLUDpUf66MbFOXBVSs7ZSw4=
Content-Language: de-DE
In-Reply-To: <udcrpm$31mqp$2@dont-email.me>
 by: Bonita Montero - Mon, 25 Sep 2023 12:54 UTC

This is my current own from_chars variant.
It is more precise than any from_chars implementation

template<std::integral Char>
conv_result_s<Char> gen_double( std::span<Char> range, double value,
int16_t nFractDigits, bool sign, std::chars_format fmt, bool round = false )
{ static_assert(sizeof(size_t) == 4 || sizeof(size_t) == 8, "must be 64
or 32 bit platform");
using namespace std;
// absolute number of fract-digits,
// negatives are a boundary for signficant bits,
// poisitives may append further zeroes
size_t absFractDigits = nFractDigits >= 0 ? nFractDigits : -nFractDigits;
constexpr uint64_t
IMPLICIT_BIT = 1ull << 52,
MANTISSA_MASK = IMPLICIT_BIT - 1;
// binary representation of value
uint64_t const bin = bit_cast<uint64_t>( value );
// signed binary representation of value
int64_t const sBin = bin;
int exp = bin >> 52 & 0x7FF;
// error-return for values that don't fit into range
auto tooLarge = [&]( auto it ) { return conv_result_s<Char>(
errc::value_too_large, it ); };
// processing Inf / (Q)NaN
auto special = [&]( auto, string_view const &svSpecial ) ->
conv_result_s<Char>
{
auto gen = range.begin();
// range has enough space ?
if( range.size() < sign + svSpecial.length() )
return tooLarge( gen );
// prepend sign
if( sign )
*gen++ = "+-"[sBin < 0];
// append svSpecial and return
return { errc(), copy( svSpecial.begin(), svSpecial.end(), gen ) };
};
// Inf / NaN exponent ?
if( exp == 0x7FF ) [[unlikely]]
if( !(bin & MANTISSA_MASK) )
// Inf
return special( []() {}, "Inf" );
else
{
// SNaN / QNaN
static char const *nanTable[2] = { "SNaN", "QNaN" };
char const *nan = nanTable[(size_t)(bin >> 51 & 1)];
return special( []() {}, { nan, nan + 4 } );
}
// mask physical mantissa bits
uint64_t mantissa = bin & MANTISSA_MASK;
// exponent is != 0: normalized number
if( exp ) [[likely]]
// prepend implicit one
mantissa = IMPLICIT_BIT | mantissa,
// make exponent signed
exp -= 0x3FF;
else
{
if( !mantissa )
{
// +/-0.0
auto gen = range.begin();
// additional digits for exponent
size_t additional = fmt == chars_format::scientific ? 3 : 0;
// enough space for digits ?
if( range.size() >= sign + 2 + absFractDigits + additional )
// return err
return tooLarge( gen );
// prepend sign
if( sign )
*gen++ = "+-"[sBin < 0];
// leading zero, maybe separator
gen = copy_n( "0.", 1 + (bool)absFractDigits, gen );
// fill fract-digits
gen = fill_n( gen, absFractDigits, '0' );
// append exponent, if scientific
gen = copy_n( "e00", additional, gen );
return { errc(), gen };
}
// normalize mantissa
unsigned normShift = countl_zero( mantissa ) - 11;
mantissa <<= normShift;
// correct exponent according to normalization
exp = -0x3FE - normShift;
}
// reserve enough space for maximum number of digits
thread_local vector<ndi<uint8_t>>
prefixDigits( 0x400 ),
suffixDigits( 0x400 );
auto prefDigHead = prefixDigits.end();
bool skipFractZeroes = false;
// has integer-bits ?
if( exp >= 0 )
{
#if defined(_MSC_VER) && defined(_M_X64)
constexpr bool HAS128 = true;
using Word = uint64_t;
using Dword = void;
constexpr unsigned WORD_DIGITS = 19;
#else
constexpr bool HAS128 = false;
constexpr bool BITS64 = sizeof(size_t) == 8;
using Word = conditional_t<sizeof(size_t) == 8, uint32_t, uint16_t>;
using Dword = size_t;
constexpr Word LARGE_DIV = BITS64 ? 1'000'000'000 : 10'000;
constexpr unsigned WORD_DIGITS = BITS64 ? 9 : 4;
#endif
constexpr unsigned WORD_BITS = sizeof(Word) * 8;
size_t smallPrefix;
unsigned tzCnt;
uint64_t prefix;
if( exp < 52 ) [[likely]]
prefix = mantissa >> 52 - exp,
tzCnt = 0;
else
prefix = mantissa,
tzCnt = exp - 52;
// mantissa, including logical suffix, doesn't fit into a Dword ?
if( exp >= sizeof(size_t) * 8 ) [[unlikely]]
{
// prefix is a large prefix, do the large integer math;
// prefix bits, survives call but not thread termination
thread_local vector<ndi<Word>> prefixBits( (0x401 + WORD_BITS - 1) /
WORD_BITS );
// fill leading zero uint32_ts
auto prefixHead = fill_n( prefixBits.begin(), tzCnt / WORD_BITS, 0 );
// begin of sigificant bits
auto prefixSignificant = prefixHead;
// prefix shift inside uint32_t
unsigned shift = tzCnt % WORD_BITS;
// write the physical prefix shifted into the logical prefix
if constexpr( !HAS128 )
do
*prefixHead++ = (Word)(prefix << shift),
prefix >>= WORD_BITS - shift,
shift = 0;
while( prefix );
else
if( shift )
{
*prefixHead++ = (uint64_t)(prefix << shift);
if( uint64_t high = prefix >> WORD_BITS - shift; high )
*prefixHead++ = high;
}
else
*prefixHead++ = prefix;
for( ; ; )
{
// remainder of our division-row
Word remainder = 0;
// start the division from the most significant bits on
auto bitsSweep = prefixHead;
auto divWord = [&]( auto ) -> bool
{
Word div;
#if defined(_MSC_VER) && defined(_M_X64)
div = div1e19( { *--bitsSweep, remainder }, &remainder );
#else
// compute division result and remainder for sub-prefix
Dword counter = (Dword)remainder << WORD_BITS | *--bitsSweep;
div = (Word)(counter / LARGE_DIV);
remainder = (Word)(counter % LARGE_DIV);
#endif
*bitsSweep = div;
return div;
};
// divide the highest prefix half_size_t;
// reduce logical prefix if result was zero
prefixHead -= !divWord( []() {} );
// do the divides for the rest of the logical prefix
while( bitsSweep != prefixSignificant )
divWord( []() {} );
auto unroll = []<size_t ... Indices>( index_sequence<Indices ...>,
auto fn )
{
(fn.template operator ()<Indices>(), ...);
};
if( bitsSweep == prefixBits.begin() ) [[likely]]
// produce HST_DIGITS digits unrolled from the remainder
unroll( make_index_sequence<WORD_DIGITS>(),
[&]<size_t>()
{
*--prefDigHead = (uint8_t)(remainder % 10u);
remainder /= 10u;
} );
else
*--prefixSignificant = remainder,
unroll( make_index_sequence<WORD_DIGITS>(),
[&]<size_t>() { *--prefDigHead = 0; } );
size_t nWords = prefixHead - prefixBits.begin();
#if defined(_MSC_VER) && defined(_M_X64)
if( nWords == 1 ) [[unlikely]]
{
smallPrefix = prefixBits[0];
break;
}
#else
// logical mantissa fits into a Dword now ?
if( prefixHead - prefixBits.begin() <= 2 ) [[likely]]
{
// check if we have two Words (2 * size_t)
assert(prefixHead - prefixBits.begin() == 2);
// assemble small prefix
smallPrefix = (Dword)prefixBits[1] << WORD_BITS | prefixBits[0];
break;
}
#endif
}
}
else
// prefix fits in size_t
smallPrefix = (size_t)(prefix << tzCnt);
// fast path for mantissa values below 0xA00000000u
do
*--prefDigHead = (uint8_t)(smallPrefix % 10u);
while( smallPrefix /= 10u );
}
else
// prefix is 0.
*--prefDigHead = 0,
// skip leading fract zeroes if we're in scientific mode
skipFractZeroes = fmt == chars_format::scientific;
auto suffDigTail = suffixDigits.begin();
// remaining fract digts; prefix might be moved partitially to suffix
size_t
nRemFractDigits = absFractDigits + round,
nPrefixes = prefixDigits.end() - prefDigHead;
if( fmt == chars_format::scientific )
nRemFractDigits -= nRemFractDigits > nPrefixes - 1 ? nPrefixes - 1 : 0;
constexpr unsigned BITS_PER_DIGIT = 0x6A4D3D; // ceil( log( 10 ) / log(
2 ) * 2 ^ 21 )
int decExp = 0;
#if defined(NDEBUG)
constexpr bool FRACT_DEBUG = false;
#else
constexpr bool FRACT_DEBUG = true;
#endif
// compute the number of fract digits which are guaranteed to be zero;
// compute fract digits only if there are more digits than the
guaranteed zeroes
if( unsigned guaranteedZeroes = ((unsigned)(exp < -1 ? -exp - 1 : 0) <<
21) / BITS_PER_DIGIT;
exp < 52 && (FRACT_DEBUG || skipFractZeroes && nRemFractDigits ||
nRemFractDigits > guaranteedZeroes) )
{
unsigned nFractBits = 52 - exp;
uint64_t suffix60;
if( nFractBits <= 60 ) [[likely]]
suffix60 = (mantissa & (1ull << nFractBits) - 1) << 60 - nFractBits;
else
{
#if defined(__SIZEOF_INT128__)
using Word = uint64_t;
using uint128_t = unsigned __int128;
#elif defined(_MSC_VER) && defined(_M_X64) && defined(XXX)
// multiplication encapsulated in mul1e19add
using Word = uint64_t;
using uint128_t = void;
#else
using Word = uint32_t;
using uint128_t = void;
#endif
// store the words as uint64_t, do the math with 128 bits
using Dword = conditional_t<sizeof(Word) == 8, uint128_t, uint64_t>;
// have a thread-local buffer which also might fit for the next time
constexpr unsigned WORD_BITS = sizeof(Word) * 8;
// largest power of 10 that fits into a Word
constexpr Word WORD_MUL = sizeof(Word) == 8 ?
(Word)10'000'000'000'000'000'000u : 1'000'000'000;
// number of digts of a WORD_MUL remainder
constexpr unsigned WORD_DIGITS = sizeof(Word) == 8 ? 19 : 9;
// number of words that fit into an uint64_t
constexpr size_t UI64_SCALE = sizeof(uint64_t) / sizeof(Word);
// thread-local buffer for the logical mantissa to survive each call
thread_local vector<Word> fractBits( (0x3FF + 52 + WORD_BITS - 1) /
WORD_BITS + 2 );
auto
bitsHigh = fractBits.begin(),
bitsLow = fractBits.begin();
// put highest set bit to bit 63
uint64_t fractMantissa = mantissa << 11;
unsigned
// bit position of highest set bit
iHsb = nFractBits - 53,
// right shift of highest set bit in logical mantissa
shift = iHsb % WORD_BITS;
//size_t nWords = (shift + 64 - countr_zero( fractMantissa ) +
WORD_BITS - 1) / WORD_BITS;
if constexpr( sizeof(Word) == 4 )
{
bool skip = true;
// write uint32-part of fractMantissa;;
// may skip the word and also increaste bitsLow if resulting
uint32_t is zero
auto advance = [&]<int Skippable>( integral_constant<int,
Skippable>, unsigned shift )
{
uint32_t word = (uint32_t)(fractMantissa << 32 - shift);
*bitsLow = word;
skip = Skippable && skip && !word;
bitsLow += skip;
++bitsHigh;
fractMantissa >>= shift;
};
advance( integral_constant<int, 2>(), shift );
advance( integral_constant<int, 1>(), 32 );
advance( integral_constant<int, 0>(), 32 );
}
else
{
if( shift )
// we have two words, write the lower word
*bitsHigh = mantissa << 64 - shift,
// increase tail if the lower word was zero
bitsLow += !*bitsHigh++;
// write upper uint64_t
*bitsHigh++ = mantissa >> shift;
}
auto bitsZero = bitsHigh;
// fill leading logical mantissa words with zero
bitsHigh = fill_n( bitsHigh, iHsb / WORD_BITS, 0 );
for( ; ; )
{
Word overflow = 0;
// sweep from the end to produce results and overflows to the left
auto bitsSweep = bitsLow;
auto mulWord = [&]( auto ) -> bool
{
#if defined(_MSC_VER) && defined(_M_X64) && defined(XXX)
uint64_t low = mul1e19add( *bitsSweep, overflow, &overflow );
*bitsSweep++ = low;
return low;

#else
// 64 or 128 bit partitial result
Dword times10 = (Dword)*bitsSweep * WORD_MUL + overflow;
// extract overfow digit
overflow = (Word)(times10 >> WORD_BITS);
// write back lower part of result
*bitsSweep++ = (Word)times10;
return (Word)times10;
#endif
};
bitsLow += !mulWord( []() {} );
while( bitsSweep != bitsZero )
mulWord( []() {} );
// sweep didn't stop at zero word ?
if( bitsSweep == bitsHigh ) [[likely]]
{
// produce all WORD_DIGITS digits at suffDigTail
for( size_t i = WORD_DIGITS; i; )
suffDigTail[--i] = (uint8_t)(overflow % 10u),
overflow /= 10u;
// number of digits to move if necessary
size_t nDigits = WORD_DIGITS;
// index of first non-zero digit if we're skipping leading zeroes
size_t iNz = 0;
// were in scientific mode and we're currently skipping leading
zeros ?
if( skipFractZeroes )
{
// skip leading zeroes of chunk
do
if( suffDigTail[iNz] )
{
skipFractZeroes = false;
break;
}
while( ++iNz != WORD_DIGITS );
// reduce number of digits to move by number of skipped zeroes
nDigits -= iNz;
// adjust decimal exponent according to number of skipped zeroes
decExp -= (int)iNz;
}
size_t nTailSkip = nDigits > nRemFractDigits ? nDigits -
nRemFractDigits : 0;
if( (nDigits -= nTailSkip) == WORD_DIGITS ) [[likely]]
// set new tail beyond the currently produced digits
suffDigTail += WORD_DIGITS,
nRemFractDigits -= WORD_DIGITS;
else
// we either had leading zeroes skipped ot tail-sipping;
// move the digits to suffDigTail
suffDigTail = copy_n( suffDigTail + iNz, nDigits, suffDigTail ),
// reduce numer of remaining fract digits accordingly
nRemFractDigits -= nDigits;
}
else
{
// we've no overflow beyond bitsHigh;
// set the lowest zero-word to overfloe
*bitsZero++ = overflow;
// n Digits to set to zero, bounded by WORD_DIGITS
size_t nDigits = nRemFractDigits > WORD_DIGITS ? WORD_DIGITS :
nRemFractDigits;
if( !skipFractZeroes )
// we're making physical fracts
suffDigTail = fill_n( suffDigTail, nDigits, 0 ),
// reduce numer of remaining fract digits accordingly
nRemFractDigits -= nDigits;
else
// lower decimal exponent by number of skipped digits
decExp -= (int)nDigits;


Click here to read the complete article
1
server_pubkey.txt

rocksolid light 0.9.8
clearnet tor