Skip to content

Commit b7b1ab4

Browse files
RavenwaterRavenwater
committed
double-double ulp function
1 parent 2ddd4b3 commit b7b1ab4

File tree

2 files changed

+119
-8
lines changed

2 files changed

+119
-8
lines changed

‎include/universal/number/dd/dd_impl.hpp

+85-7
Original file line numberDiff line numberDiff line change
@@ -849,8 +849,84 @@ inline std::string to_triple(const dd& v, int precision = 17) {
849849
return s.str();
850850
}
851851

852-
inline std::string to_binary(const dd& number, bool bNibbleMarker = false) {
852+
inline std::string to_binary(const dd& number, bool nibbleMarker = false) {
853853
std::stringstream s;
854+
855+
double_decoder decoder;
856+
decoder.d = number.high();
857+
int highExponent = static_cast<int>(decoder.parts.exponent) - ieee754_parameter<double>::bias;
858+
859+
s << "0b";
860+
// print sign bit
861+
s << (decoder.parts.sign ? '1' : '0') << '.';
862+
863+
// print exponent bits
864+
{
865+
uint64_t mask = 0x400ull;
866+
for (int bit = 10; bit >= 0; --bit) {
867+
s << ((decoder.parts.exponent & mask) ? '1' : '0');
868+
if (nibbleMarker && bit != 0 && (bit % 4) == 0) s << '\'';
869+
mask >>= 1;
870+
}
871+
}
872+
873+
s << '.';
874+
875+
// print hi fraction bits
876+
uint64_t mask = (1ull << 51);
877+
for (int bit = 51; bit >= 0; --bit) {
878+
s << ((decoder.parts.fraction & mask) ? '1' : '0');
879+
if (nibbleMarker && bit != 0 && ((bit+1) % 4) == 0) s << '\'';
880+
mask >>= 1;
881+
}
882+
883+
// print lo fraction bits
884+
decoder.d = number.low();
885+
// high limb low limb
886+
// 52 51 ..... 3210 52 51 ...... 3210
887+
// h. ffff ffff ...... ffff ffff h. ffff ffff ...... ffff ffff
888+
// 105 104 53 52 51 ...... 3210 dd_bit
889+
// | <--- exponent is exp(hi) - 53
890+
// h. ffff ffff ...... ffff ffff 0. 0000 000h. ffff ffff ...... ffff ffff
891+
// | <----- exponent would be exp(hi) - 61
892+
// h. ffff ffff ...... ffff ffff 0. 0000 0000 ...... 000h. ffff ffff ...... ffff ffff
893+
// | <----- exponent would be exp(hi) - 102
894+
// h. ffff ffff ...... ffff ffff 0. 0000 0000 ...... 0000 000h. ffff ffff ...... ffff ffff
895+
// | <----- exponent would be exp(hi) - 106
896+
// the low segment is always in normal form
897+
int lowExponent = static_cast<int>(decoder.parts.exponent) - ieee754_parameter<double>::bias;
898+
899+
assert(highExponent >= lowExponent + 53 && "exponent of lower limb is not-aligned");
900+
901+
// enumerate in the bit offset space of the double-double
902+
// that means, the first bit of the second limb is bit (105 - 53) == 52 and it cycles down to 0
903+
// representing 2^-53 through 2^-106 relative to the MSB of the high limb
904+
int offset = highExponent - 53 - lowExponent - 1;
905+
mask = (1ull << 51);
906+
s << '|'; // visual delineation between the two limbs
907+
for (int ddbit = 52; ddbit >= 0; --ddbit) {
908+
if (offset == 0) {
909+
s << (decoder.d == 0.0 ? '0' : '1'); // show hidden bit when not-zero
910+
}
911+
else if (offset > 0) {
912+
// we have to introduce a leading zero as the hidden bit is positioned at a lower ddbit offset
913+
s << '0';
914+
}
915+
else {
916+
// we have reached the fraction bits
917+
s << ((decoder.parts.fraction & mask) ? '1' : '0');
918+
mask >>= 1;
919+
}
920+
if (nibbleMarker && ddbit != 0 && (ddbit % 4) == 0) s << '\'';
921+
--offset;
922+
}
923+
924+
return s.str();
925+
}
926+
927+
inline std::string to_components(const dd& number, bool nibbleMarker = false) {
928+
std::stringstream s;
929+
s << std::setprecision(16);
854930
constexpr int nrLimbs = 2;
855931
for (int i = 0; i < nrLimbs; ++i) {
856932
double_decoder decoder;
@@ -862,27 +938,27 @@ inline std::string to_binary(const dd& number, bool bNibbleMarker = false) {
862938
// print sign bit
863939
s << (decoder.parts.sign ? '1' : '0') << '.';
864940

865-
// print exponent bits
941+
// print the segment's exponent bits
866942
{
867943
uint64_t mask = 0x400;
868944
for (int bit = 10; bit >= 0; --bit) {
869945
s << ((decoder.parts.exponent & mask) ? '1' : '0');
870-
if (bNibbleMarker && bit != 0 && (bit % 4) == 0) s << '\'';
946+
if (nibbleMarker && bit != 0 && (bit % 4) == 0) s << '\'';
871947
mask >>= 1;
872948
}
873949
}
874950

875951
s << '.';
876952

877-
// print hi fraction bits
953+
// print the segment's fraction bits
878954
uint64_t mask = (uint64_t(1) << 51);
879955
for (int bit = 51; bit >= 0; --bit) {
880956
s << ((decoder.parts.fraction & mask) ? '1' : '0');
881-
if (bNibbleMarker && bit != 0 && (bit % 4) == 0) s << '\'';
957+
if (nibbleMarker && bit != 0 && (bit % 4) == 0) s << '\'';
882958
mask >>= 1;
883959
}
884960

885-
if (i < 1) s << '\n';
961+
s << " : " << number[i] << " : binary scale " << scale(number[i]) << '\n';
886962
}
887963

888964
return s.str();
@@ -896,9 +972,11 @@ inline dd ulp(const dd& a) {
896972
double nlo;
897973
if (lo == 0.0) {
898974
nlo = std::numeric_limits<double>::epsilon() / 2.0;
975+
int binaryExponent = scale(hi) - 53;
976+
nlo /= std::pow(2.0, -binaryExponent);
899977
}
900978
else {
901-
nlo = std::nextafter(lo, INFINITY);
979+
nlo = (hi < 0.0 ? std::nextafter(lo, -INFINITY) : std::nextafter(lo, +INFINITY));
902980
}
903981
dd n(hi, nlo);
904982

‎static/dd/api/experiments.cpp

+34-1
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ namespace sw {
4242
}
4343

4444
// specialize ReportValue for double-double (dd)
45-
void ReportValue(const dd& a, const std::string& label = "", unsigned labelWidth = 20, unsigned precision = 32) {
45+
void ReportValue_(const dd& a, const std::string& label = "", unsigned labelWidth = 20, unsigned precision = 32) {
4646
auto defaultPrecision = std::cout.precision();
4747
std::cout << std::setprecision(precision);
4848
std::cout << std::setw(labelWidth) << label << " : " << a << '\n';
@@ -192,6 +192,39 @@ try {
192192
ReportValue(nlo, "new low");
193193
ReportValue(n, "n");
194194
ReportValue(n - a, "n - a");
195+
196+
std::cout << '\n';
197+
for (int i = 0; i < 10; ++i) {
198+
double a(1ull << i);
199+
double ulpAtI = ulp(a);
200+
std::string label = "ulpAt<double>(2^" + std::to_string(i) + ")";
201+
ReportValue(ulpAtI, label);
202+
}
203+
std::cout << '\n';
204+
for (int i = 0; i < 5; ++i) {
205+
dd a(1ull << i);
206+
dd ulpAtI = ulp(a);
207+
std::string label = "ulpAt<dd>(2^" + std::to_string(i) + ")";
208+
ReportValue(ulpAtI, label, 20, 32);
209+
}
210+
std::cout << std::right << std::setw(20) << "......." << " :\n" << std::internal;
211+
for (int i = 53; i < 64; ++i) {
212+
dd a(1ull << i);
213+
dd ulpAtI = ulp(a);
214+
std::string label = "ulpAt<dd>(2^" + std::to_string(i) + ")";
215+
ReportValue(ulpAtI, label, 20, 32);
216+
}
217+
std::cout << '\n';
218+
std::cout << " with a non-zero low segment\n";
219+
for (int i = 0; i < 5; ++i) {
220+
dd a(1ull << i);
221+
dd ulpAtI = ulp(a);
222+
ulpAtI += a;
223+
std::string label = "ulpAt<dd>(2^" + std::to_string(i) + "+ulp)";
224+
ReportValue(ulpAtI, label, 20, 32);
225+
// std::cout << to_components(ulpAtI) << std::setprecision(32) << ulpAtI << std::setprecision(defaultPrecision) << '\n';
226+
// std::cout << to_binary(ulpAtI, true) << '\n';
227+
}
195228
}
196229

197230
return 0;

0 commit comments

Comments
 (0)