1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
|
/**CFile****************************************************************
FileName [utilDouble.h]
SystemName [ABC: Logic synthesis and verification system.]
PackageName []
Synopsis [Double floating point number implementation.]
Author [Alan Mishchenko, Bruno Schmitt]
Affiliation [UC Berkeley / UFRGS]
Date [Ver. 1.0. Started - February 11, 2017.]
Revision []
***********************************************************************/
#ifndef ABC__sat__Xdbl__Xdbl_h
#define ABC__sat__Xdbl__Xdbl_h
#include "misc/util/abc_global.h"
ABC_NAMESPACE_HEADER_START
////////////////////////////////////////////////////////////////////////
/// STRUCTURE DEFINITIONS ///
////////////////////////////////////////////////////////////////////////
/*
The xdbl floating-point number is represented as a 64-bit unsigned int.
The number is (2^Exp)*Mnt, where Exp is a 16-bit exponent and Mnt is a
48-bit mantissa. The decimal point is located between the MSB of Mnt,
which is always 1, and the remaining 15 digits of Mnt.
Currently, only positive numbers are represented.
The range of possible values is [1.0; 2^(2^16-1)*1.111111111111111]
that is, the smallest possible number is 1.0 and the largest possible
number is 2^(---16 ones---).(1.---47 ones---)
Comparison of numbers can be done by comparing the underlying unsigned ints.
Only addition, multiplication, and division by 2^n are currently implemented.
*/
typedef word xdbl;
static inline word Xdbl_Exp( xdbl a ) { return a >> 48; }
static inline word Xdbl_Mnt( xdbl a ) { return (a << 16) >> 16; }
static inline xdbl Xdbl_Create( word Exp, word Mnt ) { assert(!(Exp>>16) && (Mnt>>47)==(word)1); return (Exp<<48) | Mnt; }
static inline xdbl Xdbl_Const1() { return Xdbl_Create( (word)0, (word)1 << 47 ); }
static inline xdbl Xdbl_Const2() { return Xdbl_Create( (word)1, (word)1 << 47 ); }
static inline xdbl Xdbl_Const3() { return Xdbl_Create( (word)1, (word)3 << 46 ); }
static inline xdbl Xdbl_Const12() { return Xdbl_Create( (word)3, (word)3 << 46 ); }
static inline xdbl Xdbl_Const1point5() { return Xdbl_Create( (word)0, (word)3 << 46 ); }
static inline xdbl Xdbl_Const2point5() { return Xdbl_Create( (word)1, (word)5 << 45 ); }
static inline xdbl Xdbl_Maximum() { return ~(word)0; }
static inline double Xdbl_ToDouble( xdbl a ) { assert(Xdbl_Exp(a) < 1023); return Abc_Word2Dbl(((Xdbl_Exp(a) + 1023) << 52) | (((a<<17)>>17) << 5)); }
static inline xdbl Xdbl_FromDouble( double a ) { word A = Abc_Dbl2Word(a); assert(a >= 1.0); return Xdbl_Create((A >> 52)-1023, (((word)1) << 47) | ((A << 12) >> 17)); }
////////////////////////////////////////////////////////////////////////
/// FUNCTION DEFINITIONS ///
////////////////////////////////////////////////////////////////////////
/**Function*************************************************************
Synopsis [Adding two floating-point numbers.]
Description []
SideEffects []
SeeAlso []
***********************************************************************/
static inline xdbl Xdbl_Add( xdbl a, xdbl b )
{
word Exp, Mnt;
if ( a < b ) a ^= b, b ^= a, a ^= b;
assert( a >= b );
Mnt = Xdbl_Mnt(a) + (Xdbl_Mnt(b) >> (Xdbl_Exp(a) - Xdbl_Exp(b)));
Exp = Xdbl_Exp(a);
if ( Mnt >> 48 ) // new MSB is created
Exp++, Mnt >>= 1;
if ( Exp >> 16 ) // overflow
return Xdbl_Maximum();
return Xdbl_Create( Exp, Mnt );
}
/**Function*************************************************************
Synopsis [Multiplying two floating-point numbers.]
Description []
SideEffects []
SeeAlso []
***********************************************************************/
static inline xdbl Xdbl_Mul( xdbl a, xdbl b )
{
word Exp, Mnt, MntA, MntB, MntAh, MntBh, MntAl, MntBl;
if ( a < b ) a ^= b, b ^= a, a ^= b;
assert( a >= b );
MntA = Xdbl_Mnt(a);
MntB = Xdbl_Mnt(b);
MntAh = MntA>>32;
MntBh = MntB>>32;
MntAl = (MntA<<32)>>32;
MntBl = (MntB<<32)>>32;
Mnt = ((MntAh * MntBh) << 17) + ((MntAl * MntBl) >> 47) + ((MntAl * MntBh) >> 15) + ((MntAh * MntBl) >> 15);
Exp = Xdbl_Exp(a) + Xdbl_Exp(b);
if ( Mnt >> 48 ) // new MSB is created
Exp++, Mnt >>= 1;
if ( Exp >> 16 ) // overflow
return Xdbl_Maximum();
return Xdbl_Create( Exp, Mnt );
}
/**Function*************************************************************
Synopsis [Dividing floating point number by a degree of 2.]
Description []
SideEffects []
SeeAlso []
***********************************************************************/
static inline xdbl Xdbl_Div( xdbl a, unsigned Deg2 )
{
if ( Xdbl_Exp(a) >= (word)Deg2 )
return Xdbl_Create( Xdbl_Exp(a) - Deg2, Xdbl_Mnt(a) );
return Xdbl_Const1(); // underflow
}
/**Function*************************************************************
Synopsis [Testing procedure.]
Description [Helpful link https://www.h-schmidt.net/FloatConverter/IEEE754.html]
SideEffects []
SeeAlso []
***********************************************************************/
static inline void Xdbl_Test()
{
xdbl c1 = Xdbl_Const1();
xdbl c2 = Xdbl_Const2();
xdbl c3 = Xdbl_Const3();
xdbl c12 = Xdbl_Const12();
xdbl c1p5 = Xdbl_Const1point5();
xdbl c2p5 = Xdbl_Const2point5();
xdbl c1_ = Xdbl_FromDouble(1.0);
xdbl c2_ = Xdbl_FromDouble(2.0);
xdbl c3_ = Xdbl_FromDouble(3.0);
xdbl c12_ = Xdbl_FromDouble(12.0);
xdbl c1p5_ = Xdbl_FromDouble(1.5);
xdbl c2p5_ = Xdbl_FromDouble(2.5);
xdbl sum1 = Xdbl_Add(c1, c1p5);
xdbl mul1 = Xdbl_Mul(c2, c1p5);
xdbl sum2 = Xdbl_Add(c1p5, c2p5);
xdbl mul2 = Xdbl_Mul(c1p5, c2p5);
xdbl a = Xdbl_FromDouble(1.2929725);
xdbl b = Xdbl_FromDouble(10.28828287);
xdbl ab = Xdbl_Mul(a, b);
xdbl ten100 = Xdbl_FromDouble( 1e100 );
xdbl ten100_ = ABC_CONST(0x014c924d692ca61b);
assert( ten100 == ten100_ );
// float f1 = Xdbl_ToDouble(c1);
// Extra_PrintBinary( stdout, (int *)&c1, 32 ); printf( "\n" );
// Extra_PrintBinary( stdout, (int *)&f1, 32 ); printf( "\n" );
printf( "1 = %lf\n", Xdbl_ToDouble(c1) );
printf( "2 = %lf\n", Xdbl_ToDouble(c2) );
printf( "3 = %lf\n", Xdbl_ToDouble(c3) );
printf( "12 = %lf\n", Xdbl_ToDouble(c12) );
printf( "1.5 = %lf\n", Xdbl_ToDouble(c1p5) );
printf( "2.5 = %lf\n", Xdbl_ToDouble(c2p5) );
printf( "Converted 1 = %lf\n", Xdbl_ToDouble(c1_) );
printf( "Converted 2 = %lf\n", Xdbl_ToDouble(c2_) );
printf( "Converted 3 = %lf\n", Xdbl_ToDouble(c3_) );
printf( "Converted 12 = %lf\n", Xdbl_ToDouble(c12_) );
printf( "Converted 1.5 = %lf\n", Xdbl_ToDouble(c1p5_) );
printf( "Converted 2.5 = %lf\n", Xdbl_ToDouble(c2p5_) );
printf( "1.0 + 1.5 = %lf\n", Xdbl_ToDouble(sum1) );
printf( "2.0 * 1.5 = %lf\n", Xdbl_ToDouble(mul1) );
printf( "1.5 + 2.5 = %lf\n", Xdbl_ToDouble(sum2) );
printf( "1.5 * 2.5 = %lf\n", Xdbl_ToDouble(mul2) );
printf( "12 / 2^2 = %lf\n", Xdbl_ToDouble(Xdbl_Div(c12, 2)) );
printf( "12 / 2^2 = %lf\n", Xdbl_ToDouble(Xdbl_Div(c12, 2)) );
printf( "%.16lf * %.16lf = %.16lf (%.16lf)\n", Xdbl_ToDouble(a), Xdbl_ToDouble(b), Xdbl_ToDouble(ab), 1.2929725 * 10.28828287 );
assert( sum1 == c2p5 );
assert( mul1 == c3 );
}
ABC_NAMESPACE_HEADER_END
#endif
|
