Skip to content

Commit b122cb2

Browse files
committed
adding a double-double command line utility
1 parent c4752c9 commit b122cb2

File tree

5 files changed

+131
-20
lines changed

5 files changed

+131
-20
lines changed

include/universal/number/dd/dd_impl.hpp

+21-9
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,9 @@ namespace sw { namespace universal {
3131

3232
// fwd references to free functions used in to_digits()
3333
dd operator*(const dd& lhs, const dd&);
34-
dd pown(dd const&, int);
34+
std::ostream& operator<<(std::ostream&, const dd&);
35+
dd pown(const dd&, int);
36+
dd frexp(const dd&, int*);
3537

3638
// dd is an unevaluated pair of IEEE-754 doubles that provides a (1,11,106) floating-point triple
3739
class dd {
@@ -811,6 +813,16 @@ inline std::string to_pair(const dd& v, int precision = 17) {
811813
return s.str();
812814
}
813815

816+
inline std::string to_triple(const dd& v, int precision = 17) {
817+
std::stringstream s;
818+
bool isneg = v.isneg();
819+
int scale = v.scale();
820+
int exponent;
821+
dd fraction = frexp(v, &exponent);
822+
s << '(' << (isneg ? '1' : '0') << ", " << scale << ", " << std::setprecision(precision) << fraction << ')';
823+
return s.str();
824+
}
825+
814826
inline std::string to_binary(const dd& number, bool bNibbleMarker = false) {
815827
std::stringstream s;
816828
constexpr int nrLimbs = 2;
@@ -863,7 +875,7 @@ inline dd abs(dd a) {
863875
return dd(hi, lo);
864876
}
865877

866-
inline dd ceil(dd const& a)
878+
inline dd ceil(const dd& a)
867879
{
868880
if (a.isnan()) return a;
869881

@@ -878,7 +890,7 @@ inline dd ceil(dd const& a)
878890
return dd(hi, lo);
879891
}
880892

881-
inline dd floor(dd const& a) {
893+
inline dd floor(const dd& a) {
882894
if (a.isnan()) return a;
883895

884896
double hi = std::floor(a.high());
@@ -974,7 +986,7 @@ inline dd mul_pwr2(const dd& a, double b) {
974986
// quad-double operators
975987

976988
// quad-double + double-double
977-
void qd_add(double const a[4], dd const& b, double s[4]) {
989+
void qd_add(double const a[4], const dd& b, double s[4]) {
978990
double t[5];
979991
s[0] = two_sum(a[0], b.high(), t[0]); // s0 - O( 1 ); t0 - O( e )
980992
s[1] = two_sum(a[1], b.low(), t[1]); // s1 - O( e ); t1 - O( e^2 )
@@ -991,7 +1003,7 @@ void qd_add(double const a[4], dd const& b, double s[4]) {
9911003
}
9921004

9931005
// quad-double = double-double * double-double
994-
void qd_mul(dd const& a, dd const& b, double p[4]) {
1006+
void qd_mul(const dd& a, const dd& b, double p[4]) {
9951007
double p4, p5, p6, p7;
9961008

9971009
// powers of e - 0, 1, 1, 1, 2, 2, 2, 3
@@ -1025,15 +1037,15 @@ void qd_mul(dd const& a, dd const& b, double p[4]) {
10251037
}
10261038
}
10271039

1028-
inline dd fma(dd const& a, dd const& b, dd const& c) {
1040+
inline dd fma(const dd& a, const dd& b, const dd& c) {
10291041
double p[4];
10301042
qd_mul(a, b, p);
10311043
qd_add(p, c, p);
10321044
p[0] = two_sum(p[0], p[1] + p[2] + p[3], p[1]);
10331045
return dd(p[0], p[1]);
10341046
}
10351047

1036-
inline dd sqr(dd const& a) {
1048+
inline dd sqr(const dd& a) {
10371049
if (a.isnan()) return a;
10381050

10391051
double p2, p1 = two_sqr(a.high(), p2);
@@ -1044,7 +1056,7 @@ inline dd sqr(dd const& a) {
10441056
return dd(s1, s2);
10451057
}
10461058

1047-
inline dd reciprocal(dd const& a) {
1059+
inline dd reciprocal(const dd& a) {
10481060
if (a.iszero()) return dd(SpecificValue::infpos);
10491061

10501062
if (a.isinf()) return dd(0.0);
@@ -1065,7 +1077,7 @@ inline dd reciprocal(dd const& a) {
10651077
}
10661078
}
10671079

1068-
inline dd pown(dd const& a, int n) {
1080+
inline dd pown(const dd& a, int n) {
10691081
if (a.isnan()) return a;
10701082

10711083
int N = (n < 0) ? -n : n;

include/universal/number/dd/numeric_limits.hpp

+3-7
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,6 @@ class numeric_limits< sw::universal::dd > {
1515
using DoubleDouble = sw::universal::dd;
1616
static constexpr bool is_specialized = true;
1717
static constexpr DoubleDouble min() { // return minimum value
18-
// return DoubleDouble(sw::universal::SpecificValue::minpos);
1918
return DoubleDouble(radix * (numeric_limits< double >::min() / numeric_limits< double >::epsilon()));
2019
}
2120
static constexpr DoubleDouble max() { // return maximum value
@@ -27,26 +26,23 @@ class numeric_limits< sw::universal::dd > {
2726
return (-(max)());
2827
}
2928
static constexpr DoubleDouble epsilon() { // return smallest effective increment from 1.0
30-
return numeric_limits< double >::epsilon() * numeric_limits< double >::epsilon() / radix;
29+
constexpr double epsilon{ std::numeric_limits< double >::epsilon() };
30+
return (epsilon * epsilon) * 0.5;
3131
}
3232
static constexpr DoubleDouble round_error() { // return largest rounding error
3333
return DoubleDouble(1.0 / radix);
3434
}
3535
static constexpr DoubleDouble denorm_min() { // return minimum denormalized value
36-
// return DoubleDouble(sw::universal::SpecificValue::minpos);
37-
return 0.0;
36+
return DoubleDouble(std::numeric_limits<double>::denorm_min());
3837
}
3938
static constexpr DoubleDouble infinity() { // return positive infinity
4039
return DoubleDouble(sw::universal::SpecificValue::infpos);
41-
//return numeric_limits< double >::infinity();
4240
}
4341
static constexpr DoubleDouble quiet_NaN() { // return non-signaling NaN
4442
return DoubleDouble(sw::universal::SpecificValue::qnan);
45-
//return numeric_limits< double >::quiet_NaN();
4643
}
4744
static constexpr DoubleDouble signaling_NaN() { // return signaling NaN
4845
return DoubleDouble(sw::universal::SpecificValue::snan);
49-
//return numeric_limits< double >::signaling_NaN();
5046
}
5147

5248
static constexpr int digits = 2 * std::numeric_limits<double>::digits;

static/dd/api/api.cpp

+28
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,34 @@ namespace sw {
6161
ostr << str << '\n';
6262
}
6363

64+
void construct_largest_doubledouble() {
65+
using Scalar = dd;
66+
67+
double firstLimb = std::numeric_limits<double>::max();
68+
dd a = std::numeric_limits<Scalar>::max();
69+
std::cout << std::setprecision(32) << a << '\n';
70+
int expOfFirstLimb = scale(a);
71+
std::cout << to_binary(expOfFirstLimb) << " : " << expOfFirstLimb << '\n';
72+
// second limb exponent
73+
int expOfSecondLimb = expOfFirstLimb - std::log10(1ull << 53);
74+
std::cout << "exponent of the first limb : " << expOfFirstLimb << '\n';
75+
std::cout << "exponent of the second limb : " << expOfSecondLimb << '\n';
76+
// construct the second limb
77+
double secondLimb = std::ldexp(1.0, expOfSecondLimb);
78+
std::cout << "1.0 " << to_binary(1.0) << '\n';
79+
std::cout << "first limb " << to_binary(firstLimb) << '\n';
80+
std::cout << "second limb " << to_binary(secondLimb) << '\n';
81+
82+
dd aa(firstLimb, secondLimb);
83+
std::cout << std::setprecision(16) << firstLimb << '\n';
84+
std::cout << std::setprecision(16) << aa << '\n';
85+
std::cout << std::setprecision(32) << aa << '\n';
86+
87+
dd b = ulp(std::numeric_limits<double>::max());
88+
dd c = a + b;
89+
std::cout << c << '\n';
90+
}
91+
6492
}
6593
}
6694

tools/cmd/doubledouble.cpp

+74
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
// doubledouble.cpp: components of a double-double: cli to show the sign/scale/limb components of a double-double floating-point
2+
//
3+
// Copyright (C) 2017 Stillwater Supercomputing, Inc.
4+
// SPDX-License-Identifier: MIT
5+
//
6+
// This file is part of the universal numbers project, which is released under an MIT Open Source license.
7+
#include <universal/utility/directives.hpp>
8+
#include <universal/utility/bit_cast.hpp>
9+
#include <limits>
10+
#include <universal/number/dd/dd.hpp>
11+
#include <universal/common/number_traits_reports.hpp>
12+
13+
// ShowRepresentations prints the different output formats for the long double type
14+
template<typename Scalar>
15+
void ShowRepresentations(std::ostream& ostr, const Scalar& f) {
16+
using namespace sw::universal;
17+
auto oldprec = ostr.precision(); // save stream state
18+
19+
constexpr int max_digits10 = std::numeric_limits<Scalar>::max_digits10; // floating-point attribute for printing scientific format
20+
21+
Scalar v(f); // convert to target cfloat
22+
ostr << "scientific : " << std::setprecision(max_digits10) << v << '\n';
23+
ostr << "triple form : " << to_triple(v) << '\n';
24+
ostr << "binary form : " << to_binary(v, true) << '\n';
25+
ostr << "color coded : " << color_print(v) << '\n';
26+
27+
ostr << std::setprecision(oldprec);
28+
}
29+
30+
// receive a float and print the components of a long double representation
31+
int main(int argc, char** argv)
32+
try {
33+
using namespace sw::universal;
34+
using Scalar = dd;
35+
36+
if (argc != 2) {
37+
std::cerr << "doubledouble: components of a double-double floating-point\n";
38+
std::cerr << "Show the sign/scale/fraction components of an double-double.\n";
39+
std::cerr << "Usage: doubledouble fp_value_string\n";
40+
std::cerr << "Example: doubledouble 0.03124999\n";
41+
ShowRepresentations<Scalar>(std::cerr, 0.03124999);
42+
43+
std::cout << "Number Traits of a double-double\n";
44+
numberTraits<Scalar>(std::cout);
45+
46+
std::cout << "largest normal number\n";
47+
std::cout << to_binary(std::numeric_limits<Scalar>::max()) << '\n';
48+
std::cout << "smallest normal number\n";
49+
std::cout << to_binary(std::numeric_limits<Scalar>::min()) << '\n';
50+
std::cout << "smallest denormalized number\n";
51+
std::cout << to_binary(std::numeric_limits<Scalar>::denorm_min()) << '\n';
52+
53+
constexpr Scalar epsilon{ std::numeric_limits< Scalar >::epsilon() };
54+
std::cout << "epsilon : " << epsilon << '\n';
55+
std::cout << to_binary(epsilon) << '\n';
56+
57+
std::cout.flush();
58+
return EXIT_SUCCESS; // signal successful completion for ctest
59+
}
60+
61+
dd doubledouble(argv[1]);
62+
ShowRepresentations<Scalar>(std::cout, doubledouble);
63+
64+
std::cout.flush();
65+
return EXIT_SUCCESS;
66+
}
67+
catch (const char* const msg) {
68+
std::cerr << msg << std::endl;
69+
return EXIT_FAILURE;
70+
}
71+
catch (...) {
72+
std::cerr << "Caught unknown exception" << std::endl;
73+
return EXIT_FAILURE;
74+
}

tools/cmd/quaddouble.cpp

+5-4
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212

1313
// ShowRepresentations prints the different output formats for the quad-double type
1414
template<typename Scalar>
15-
void ShowRepresentations(std::ostream& ostr, sw::universal::qd f) {
15+
void ShowRepresentations(std::ostream& ostr, const Scalar& f) {
1616
using namespace sw::universal;
1717
auto defaultPrecision = ostr.precision(); // save stream state
1818

@@ -43,12 +43,12 @@ try {
4343

4444
if (argc != 2) {
4545
std::cerr << "quaddouble: components of a quad-double floating-point\n";
46-
std::cerr << "Show the sign/scale/limbs components of a quad-double.\n";
46+
std::cerr << "Show the sign/scale/fraction components of a quad-double.\n";
4747
std::cerr << "Usage: quaddouble fp_value_string\n";
4848
std::cerr << "Example: quaddouble 0.03124999\n";
4949
ShowRepresentations<Scalar>(std::cerr, 0.03124999);
5050

51-
std::cout << "Number Traits of quad-double\n";
51+
std::cout << "Number Traits of a quad-double\n";
5252
numberTraits<Scalar>(std::cout);
5353

5454
std::cout << "largest normal number\n";
@@ -58,9 +58,10 @@ try {
5858
std::cout << "smallest denormalized number\n";
5959
std::cout << to_binary(std::numeric_limits<Scalar>::denorm_min()) << '\n';
6060

61-
Scalar epsilon{ std::numeric_limits< Scalar >::epsilon() };
61+
constexpr Scalar epsilon{ std::numeric_limits< Scalar >::epsilon() };
6262
std::cout << "epsilon : " << epsilon << '\n';
6363
std::cout << to_binary(epsilon) << '\n';
64+
6465
std::cout.flush();
6566
return EXIT_SUCCESS; // signal successful completion for ctest
6667
}

0 commit comments

Comments
 (0)