Floating Point Arithmetic

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct DoubleParts {
5
    uint64_t sign;      // 0 or 1
6
    uint64_t exponent;  // 11 bits
7
    uint64_t fraction;  // 52 bits
8
};
9
10
DoubleParts dissectDouble(double x) {
11
    uint64_t bits = 0;
12
    static_assert(sizeof(double) == sizeof(uint64_t), "Unexpected double size");
13
    memcpy(&bits, &x, sizeof(bits)); // Safe bit reinterpretation
14
    DoubleParts p{};
15
    p.sign = (bits >> 63) & 0x1ull;
16
    p.exponent = (bits >> 52) & 0x7ffull; // 11 bits
17
    p.fraction = bits & 0x000fffffffffffffull; // 52 bits
18
    return p;
19
}
20
21
int main() {
22
    cout.setf(std::ios::boolalpha);
23
24
    vector<double> samples = {1.0, -0.0, 0.1, numeric_limits<double>::infinity(), -numeric_limits<double>::infinity(), 
25
                              numeric_limits<double>::quiet_NaN(), numeric_limits<double>::denorm_min()};
26
27
    for (double x : samples) {
28
        DoubleParts p = dissectDouble(x);
29
        cout << "x = " << x << "\n";
30
        cout << "  sign:      " << p.sign << "\n";
31
        cout << "  exponent:  0x" << hex << setw(3) << setfill('0') << p.exponent << dec
32
             << " (" << p.exponent << ")\n";
33
        cout << "  fraction:  0x" << hex << setw(13) << setfill('0') << p.fraction << dec << "\n";
34
35
        // Classify special cases
36
        bool isZero      = (p.exponent == 0 && p.fraction == 0);
37
        bool isSubnormal = (p.exponent == 0 && p.fraction != 0);
38
        bool isInf       = (p.exponent == 0x7ff && p.fraction == 0);
39
        bool isNaN       = (p.exponent == 0x7ff && p.fraction != 0);
40
41
        cout << "  zero?       " << isZero << "\n";
42
        cout << "  subnormal?  " << isSubnormal << "\n";
43
        cout << "  infinity?   " << isInf << "\n";
44
        cout << "  NaN?        " << isNaN << "\n";
45
46
        if (!isZero && !isSubnormal && !isInf && !isNaN) {
47
            int64_t unbiased = static_cast<int64_t>(p.exponent) - 1023; // bias for double
48
            cout << "  unbiased exponent (e - bias): " << unbiased << "\n";
49
        }
50
        cout << "\n";
51
    }
52
}
53

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
template <class T>
5
T naive_sum(const vector<T>& a) {
6
    T s = 0;
7
    for (T v : a) s += v; // Each addition rounds; errors may accumulate
8
    return s;
9
}
10
11
template <class T>
12
T kahan_sum(const vector<T>& a) {
13
    T s = 0;      // Running sum
14
    T c = 0;      // Compensation for lost low-order bits
15
    for (T v : a) {
16
        T y = v - c;      // First, subtract the compensation
17
        T t = s + y;      // Add the corrected value
18
        c = (t - s) - y;  // New compensation captures low-order roundoff
19
        s = t;
20
    }
21
    return s;
22
}
23
24
int main() {
25
    // Construct a numerically challenging sum: one big number plus many tiny ones
26
    const int n = 1000000;
27
    vector<double> a; a.reserve(n + 1);
28
    a.push_back(1e9);
29
    for (int i = 0; i < n; ++i) a.push_back(1e-6);
30
31
    double s1 = naive_sum(a);
32
    double s2 = kahan_sum(a);
33
34
    cout.setf(std::ios::fixed); cout << setprecision(10);
35
    cout << "Naive sum:  " << s1 << "\n";
36
    cout << "Kahan sum:  " << s2 << "\n";
37
    cout << "Expected (real): 1000000000.999999\n";
38
    cout << "Absolute error naive: " << fabs(s1 - 1000000000.999999) << "\n";
39
    cout << "Absolute error Kahan: " << fabs(s2 - 1000000000.999999) << "\n";
40
}
41

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Absolute/relative tolerance comparison
5
bool approx_equal(double a, double b, double rel_tol = 1e-12, double abs_tol = 1e-15) {
6
    // Handle infinities and NaNs explicitly
7
    if (std::isnan(a) || std::isnan(b)) return false; // NaN is never equal
8
    if (std::isinf(a) || std::isinf(b)) return a == b; // only equal if same infinity
9
10
    double diff = fabs(a - b);
11
    if (diff <= abs_tol) return true; // near zero
12
13
    double max_ab = max(fabs(a), fabs(b));
14
    return diff <= rel_tol * max_ab;
15
}
16
17
// ULP distance for IEEE 754 double (orders all bit patterns monotonically)
18
static inline uint64_t to_ordered_uint(double x) {
19
    uint64_t u; memcpy(&u, &x, sizeof(u));
20
    // Map negatives to high range so that increasing numeric values have increasing ordered integers
21
    return (u & (1ull<<63)) ? (~u + 1ull) : (u | (1ull<<63));
22
}
23
24
uint64_t ulp_distance(double a, double b) {
25
    if (std::isnan(a) || std::isnan(b)) return UINT64_MAX; // undefined distance
26
    uint64_t ua = to_ordered_uint(a);
27
    uint64_t ub = to_ordered_uint(b);
28
    return (ua > ub) ? (ua - ub) : (ub - ua);
29
}
30
31
int main() {
32
    double x = 1.0;
33
    double y = std::nextafter(1.0, 2.0); // adjacent representable number above 1.0
34
35
    cout.setf(std::ios::boolalpha);
36
    cout << "x = " << setprecision(17) << x << ", y = " << y << "\n";
37
    cout << "approx_equal(x, y) with tight tolerances: " << approx_equal(x, y, 1e-16, 0.0) << "\n";
38
    cout << "ULP distance(x, y) = " << ulp_distance(x, y) << " (should be 1)\n";
39
40
    double a = 1e9;
41
    double b = a + 0.1; // 0.1 may not change a if it rounds away
42
    cout << "ULP distance(1e9, 1e9+0.1) = " << ulp_distance(a, b) << "\n";
43
}
44