1 files changed, 133 insertions, 0 deletions
diff --git a/src/sat/satoko/utils/sdbl.h b/src/sat/satoko/utils/sdbl.h
new file mode 100644
index 00000000..9f90ba02
--- /dev/null
+++ b/src/sat/satoko/utils/sdbl.h
@@ -0,0 +1,133 @@
+//===--- sdbl.h -------------------------------------------------------------===
+//
+//                     satoko: Satisfiability solver
+//
+// This file is distributed under the BSD 2-Clause License.
+// See LICENSE for details.
+//
+//===------------------------------------------------------------------------===
+// by Alan Mishchenko
+#ifndef satoko__utils__sdbl_h
+#define satoko__utils__sdbl_h
+
+#include "misc/util/abc_global.h"
+ABC_NAMESPACE_HEADER_START
+/*
+    The sdbl_t floating-point number is represented as a 64-bit unsigned int.
+    The number is (2^expt)*mnt, where expt 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 sdbl_t;
+
+static sdbl_t SDBL_CONST1 = ABC_CONST(0x0000800000000000);
+static sdbl_t SDBL_MAX = ~(sdbl_t)(0);
+
+union ui64_dbl { word ui64; double dbl; };
+
+static inline word sdbl_exp(sdbl_t a) { return a >> 48;         }
+static inline word sdbl_mnt(sdbl_t a) { return (a << 16) >> 16; }
+
+static inline double sdbl2double(sdbl_t a) {
+    union ui64_dbl temp;
+    assert(sdbl_exp(a) < 1023);
+    temp.ui64 = ((sdbl_exp(a) + 1023) << 52) | (((a << 17) >> 17) << 5);
+    return temp.dbl;
+}
+
+static inline sdbl_t double2sdbl(double value)
+{
+    union ui64_dbl temp;
+    sdbl_t expt, mnt;
+    assert(value >= 1.0);
+    temp.dbl = value;
+    expt = (temp.ui64 >> 52) - 1023;
+    mnt = SDBL_CONST1 | ((temp.ui64 << 12) >> 17);
+    return (expt << 48) + mnt;
+}
+
+static inline sdbl_t sdbl_add(sdbl_t a, sdbl_t b)
+{
+    sdbl_t expt, mnt;
+    if (a < b) {
+        a ^= b;
+        b ^= a;
+        a ^= b;
+    }
+    assert(a >= b);
+    expt = sdbl_exp(a);
+    mnt = sdbl_mnt(a) + (sdbl_mnt(b) >> (sdbl_exp(a) - sdbl_exp(b)));
+    /* Check for carry */
+    if (mnt >> 48) {
+        expt++;
+        mnt >>= 1;
+    }
+    if (expt >> 16) /* overflow */
+        return SDBL_MAX;
+    return (expt << 48) + mnt;
+}
+
+static inline sdbl_t sdbl_mult(sdbl_t a, sdbl_t b)
+{
+    sdbl_t expt, mnt;
+    sdbl_t a_mnt, a_mnt_hi, a_mnt_lo;
+    sdbl_t b_mnt, b_mnt_hi, b_mnt_lo;
+    if (a < b) {
+        a ^= b;
+        b ^= a;
+        a ^= b;
+    }
+    assert( a >= b );
+    a_mnt  = sdbl_mnt(a);
+    b_mnt  = sdbl_mnt(b);
+    a_mnt_hi = a_mnt>>32;
+    b_mnt_hi = b_mnt>>32;
+    a_mnt_lo = (a_mnt<<32)>>32;
+    b_mnt_lo = (b_mnt<<32)>>32;
+    mnt = ((a_mnt_hi * b_mnt_hi) << 17) +
+          ((a_mnt_lo * b_mnt_lo) >> 47) +
+          ((a_mnt_lo * b_mnt_hi) >> 15) +
+          ((a_mnt_hi * b_mnt_lo) >> 15);
+    expt = sdbl_exp(a) + sdbl_exp(b);
+    /* Check for carry */
+    if (mnt >> 48) {
+        expt++;
+        mnt >>= 1;
+    }
+    if (expt >> 16) /* overflow */
+        return SDBL_MAX;
+    return (expt << 48) + mnt;
+}
+
+static inline sdbl_t sdbl_div(sdbl_t a, unsigned deg2)
+{
+    if (sdbl_exp(a) >= (word)deg2)
+        return ((sdbl_exp(a) - deg2) << 48) + sdbl_mnt(a);
+    return SDBL_CONST1;
+}
+
+static inline void sdbl_test()
+{
+    sdbl_t ten100_ = ABC_CONST(0x014c924d692ca61b);
+    printf("%f\n", sdbl2double(ten100_));
+    printf("%016lX\n", double2sdbl(1 /0.95));
+    printf("%016lX\n", SDBL_CONST1);
+    printf("%f\n", sdbl2double(SDBL_CONST1));
+    printf("%f\n", sdbl2double(ABC_CONST(0x000086BCA1AF286B)));
+
+}
+
+ABC_NAMESPACE_HEADER_END
+
+#endif /* satoko__utils__sdbl_h */