Skip to content

Commit c49c548

Browse files
RavenwaterRavenwater
committed
first principle construction of double-double and quad-double values
1 parent e39a314 commit c49c548

File tree

4 files changed

+87
-18
lines changed

4 files changed

+87
-18
lines changed

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

+1-1
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@
3030

3131
namespace sw { namespace universal {
3232

33+
// this is debug infrastructure: TODO: remove when decimal conversion is solved reliably
3334
constexpr bool bTraceDecimalConversion = false;
3435
constexpr bool bTraceDecimalRounding = false;
3536
std::ostream& operator<<(std::ostream& ostr, const std::vector<char>& s) {
@@ -668,7 +669,6 @@ class dd {
668669
if (iszero()) {
669670
exponent = 0;
670671
for (int i = 0; i < precision; ++i) s[i] = '0';
671-
s[precision] = 0; // termination null
672672
return;
673673
}
674674

‎include/universal/number/qd/qd_impl.hpp

+9-13
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#include <iomanip>
1818
#include <limits>
1919
#include <cmath>
20+
#include <vector>
2021

2122
// supporting types and functions
2223
#include <universal/native/ieee754.hpp>
@@ -672,7 +673,7 @@ class qd {
672673

673674
int nrDigitsForFixedFormat = nrDigits;
674675
if (fixed)
675-
nrDigitsForFixedFormat = std::max(60, nrDigits); // can be much longer than the max accuracy for double-double
676+
nrDigitsForFixedFormat = std::max(120, nrDigits); // can be much longer than the max accuracy for quad-double
676677

677678
// a number in the range of [0.5, 1.0) to be printed with zero precision
678679
// must be rounded up to 1 to print correctly
@@ -682,28 +683,28 @@ class qd {
682683
}
683684

684685
if (fixed && nrDigits <= 0) {
685-
// process values with negative exponents (powerOfTenScale < 0)
686+
// process values that are near zero
686687
s += '0';
687688
if (precision > 0) {
688689
s += '.';
689690
s.append(static_cast<unsigned int>(precision), '0');
690691
}
691692
}
692693
else {
693-
char* t;
694+
std::vector<char> t;
694695

695696
if (fixed) {
696-
t = new char[static_cast<size_t>(nrDigitsForFixedFormat + 1)];
697+
t.resize(nrDigitsForFixedFormat + 1);
697698
to_digits(t, e, nrDigitsForFixedFormat);
698699
}
699700
else {
700-
t = new char[static_cast<size_t>(nrDigits + 1)];
701+
t.resize(nrDigits + 1);
701702
to_digits(t, e, nrDigits);
702703
}
703704

704705
if (fixed) {
705706
// round the decimal string
706-
round_string(t, nrDigits, &integerDigits);
707+
round_string(t, nrDigits+1, &integerDigits);
707708

708709
if (integerDigits > 0) {
709710
int i;
@@ -727,7 +728,6 @@ class qd {
727728
s += t[i];
728729

729730
}
730-
delete[] t;
731731
}
732732
}
733733

@@ -869,7 +869,7 @@ class qd {
869869
/// functional helpers
870870

871871
// precondition: string s must be all digits
872-
void round_string(char* s, int precision, int* decimalPoint) const {
872+
void round_string(std::vector<char>& s, int precision, int* decimalPoint) const {
873873
int nrDigits = precision;
874874
// round decimal string and propagate carry
875875
int lastDigit = nrDigits - 1;
@@ -891,8 +891,6 @@ class qd {
891891
(*decimalPoint)++; // increment decimal point
892892
++precision;
893893
}
894-
895-
s[precision] = 0; // aqd termination null
896894
}
897895

898896
void append_exponent(std::string& str, int e) const {
@@ -923,17 +921,15 @@ class qd {
923921
/// <param name="s"></param>
924922
/// <param name="exponent"></param>
925923
/// <param name="precision"></param>
926-
void to_digits(char* s, int& exponent, int precision) const {
924+
void to_digits(std::vector<char>& s, int& exponent, int precision) const {
927925
constexpr qd _one(1.0), _ten(10.0);
928926
constexpr double _log2(0.301029995663981);
929927
double hi = x[0];
930928
//double lo = x[1];
931929

932930
if (iszero()) {
933-
std::cout << "I am zero\n";
934931
exponent = 0;
935932
for (int i = 0; i < precision; ++i) s[i] = '0';
936-
s[precision] = 0; // termination null
937933
return;
938934
}
939935

‎static/dd/api/experiments.cpp

+1
Original file line numberDiff line numberDiff line change
@@ -103,6 +103,7 @@ try {
103103
double exponent = -std::ceil(std::abs(std::log10(low)));
104104
std::cout << "exponent : " << exponent << '\n';
105105

106+
// now let's walk that bit down to the ULP
106107
std::cout << std::setprecision(32);
107108
for (int i = 0; i < 54; ++i) {
108109
low = (std::pow(2.0, -53.0 - double(i)));

‎static/qd/api/experiments.cpp

+76-4
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,16 @@ namespace sw {
7575
return ostr << v.v;
7676
}
7777

78+
79+
std::string centered(const std::string& label, unsigned columnWidth) {
80+
unsigned length = static_cast<unsigned>(label.length());
81+
if (columnWidth < length) return label;
82+
83+
unsigned padding = columnWidth - length;
84+
unsigned leftPadding = (padding >> 1);
85+
unsigned rightPadding = padding - leftPadding;
86+
return std::string(leftPadding, ' ') + label + std::string(rightPadding, ' ');
87+
}
7888
}
7989
}
8090

@@ -96,13 +106,75 @@ try {
96106
// dd = high + lo
97107
// = 1*2^0 + 1*2^-53
98108
// = 1.0e00 + 1.0elog10(2^-53)
99-
ReportValue(std::pow(2.0, 0.0), "2^0");
100-
ReportValue(std::pow(2.0, -53.0), "2^-53");
101-
std::cout << std::log10(std::pow(2.0, -53.0)) << '\n';
102-
double exponent = std::ceil(std::log10(std::pow(2.0, -53.0)));
109+
double high{ std::pow(2.0, 0.0) };
110+
ReportValue(high, "2^0");
111+
double low{ std::pow(2.0, -53.0) };
112+
ReportValue(low, "2^-53");
113+
std::cout << std::log10(low) << '\n';
114+
double exponent = -std::ceil(std::abs(std::log10(low)));
103115
std::cout << "exponent : " << exponent << '\n';
116+
117+
// now let's walk that bit down to the ULP
118+
double x0{ 1.0 };
119+
double x1{ 0.0 };
120+
double x2{ 0.0 };
121+
double x3{ 0.0 };
122+
int precisionForRange = 16;
123+
std::cout << std::setprecision(precisionForRange);
124+
x0 = 1.0;
125+
qd a(x0, x1, x2, x3);
126+
std::cout << centered("quad-double", precisionForRange + 6) << " : ";
127+
std::cout << centered("binary form of x0", 68) << " : ";
128+
std::cout << centered("real value of x0", 15) << '\n';
129+
std::cout << a << " : " << to_binary(x0) << " : " << x0 << '\n';
130+
for (int i = 1; i < 53; ++i) {
131+
x0 = 1.0 + (std::pow(2.0, - double(i)));
132+
qd a(x0, x1, x2, x3);
133+
std::cout << a << " : " << to_binary(x0) << " : " << std::setprecision(7) << x0 << std::setprecision(precisionForRange) << '\n';
134+
}
135+
// x0 is 1.0 + eps() at this point
136+
// std::cout << to_binary(x0) << '\n';
137+
std::cout << to_binary(qd(x0, x1, x2, x3)) << '\n';
138+
x0 = 1.0;
139+
precisionForRange = 32;
140+
std::cout << std::setprecision(precisionForRange);
141+
std::cout << centered("quad-double", precisionForRange + 6) << " : ";
142+
std::cout << centered("binary form of x1", 68) << " : ";
143+
std::cout << centered("real value of x1", 15) << '\n';
144+
for (int i = 0; i < 54; ++i) {
145+
x1 = (std::pow(2.0, -53.0 - double(i)));
146+
qd a(x0, x1, x2, x3);
147+
std::cout << a << " : " << to_binary(x1) << " : " << std::setprecision(7) << x1 << std::setprecision(precisionForRange) << '\n';
148+
}
149+
std::cout << to_binary(qd(x0, x1, x2, x3)) << '\n';
150+
x1 = 0.0;
151+
precisionForRange = 48;
152+
std::cout << std::setprecision(precisionForRange);
153+
std::cout << centered("quad-double", precisionForRange + 6) << " : ";
154+
std::cout << centered("binary form of x2", 68) << " : ";
155+
std::cout << centered("real value of x2", 15) << '\n';
156+
for (int i = 0; i < 54; ++i) {
157+
x2 = (std::pow(2.0, -106.0 - double(i)));
158+
qd a(x0, x1, x2, x3);
159+
std::cout << a << " : " << to_binary(x2) << " : " << std::setprecision(7) << x2 << std::setprecision(precisionForRange) << '\n';
160+
}
161+
std::cout << to_binary(qd(x0, x1, x2, x3)) << '\n';
162+
x2 = 0.0;
163+
precisionForRange = 64;
164+
std::cout << std::setprecision(precisionForRange);
165+
std::cout << centered("quad-double", precisionForRange + 6) << " : ";
166+
std::cout << centered("binary form of x3", 68) << " : ";
167+
std::cout << centered("real value of x3", 15) << '\n';
168+
for (int i = 0; i < 54; ++i) {
169+
x3 = (std::pow(2.0, -159.0 - double(i)));
170+
qd a(x0, x1, x2, x3);
171+
std::cout << a << " : " << to_binary(x3) << " : " << std::setprecision(7) << x3 << std::setprecision(precisionForRange) << '\n';
172+
}
173+
std::cout << to_binary(qd(x0, x1, x2, x3)) << '\n';
174+
std::cout << std::setprecision(defaultPrecision);
104175
}
105176

177+
return 0;
106178
{
107179
// what is the difference between ostream fmt scientific/fixed
108180

0 commit comments

Comments
 (0)