cpp-is-magic/examples/example_3.cpp at main · MatthieuHernandez/cpp-is-magic · GitHub

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
#include <compare>
#include <cstdint>
#include <numbers>

template <std::floating_point CT>
class more_precise {
    using T = typename std::remove_const<CT>::type;

  public:
    constexpr more_precise(const T& v = {}) {
        value_1 = v;
        split_values();
    }

    constexpr T get() const { return value_1 + value_2; }

    constexpr more_precise(const more_precise<T>&& other) { this->operator=(other); }

    constexpr more_precise(const more_precise<T>& other) { this->operator=(other); }

    constexpr more_precise<T>& operator=(const more_precise<T>&& other) {
        value_1 = other.value_1;
        value_2 = other.value_2;
        split_values();
        return *this;
    }

    constexpr more_precise<T>& operator=(const more_precise<T>& other) {
        value_1 = other.value_1;
        value_2 = other.value_2;
        split_values();
        return *this;
    }

    constexpr more_precise<T> operator+(const more_precise<T>& other) const {
        more_precise<T> result;
        result.value_1 = value_1 + other.value_1;
        result.value_2 = value_2 + other.value_2;
        return result;
    }

    constexpr more_precise<T> operator-(const more_precise<T>& other) const {
        more_precise<T> result;
        result.value_1 = value_1 - other.value_1;
        result.value_2 = value_2 - other.value_2;
        return result;
    }

    constexpr more_precise<T> operator*(const more_precise<T>& other) const {
        more_precise<T> result;
        result.value_1 = value_1 * other.value_1 + value_1 * other.value_2;
        result.value_2 = value_2 * other.value_1 + value_2 * other.value_2;
        return result;
    }

    constexpr more_precise<T> operator/(const more_precise<T>& other) const {
        more_precise<T> result;
        result.value_1 = value_1 / (other.value_1 + other.value_2);
        result.value_2 = value_2 / (other.value_1 + other.value_2);
        return result;
    }

    constexpr bool operator==(const more_precise<T>& other) const = default;

    constexpr std::partial_ordering operator<=>(const more_precise<T>& other) const {
        if (value_1 < other.value_1) {
            return std::partial_ordering::less;
        } else if (value_1 > other.value_1) {
            return std::partial_ordering::greater;
        } else if (value_1 == other.value_1 && value_2 < other.value_2) {
            return std::partial_ordering::less;
        } else if (value_1 == other.value_1 && value_2 > other.value_2) {
            return std::partial_ordering::greater;
        } else if (value_1 == other.value_1 && value_1 == other.value_2) {
            return std::partial_ordering::equivalent;
        } else {
            return std::partial_ordering::unordered;
        }
    }

  private:
    constexpr void split_values() {
        static constexpr T split_factor = 128 * sizeof(T);
        if ((value_2 < value_1 && value_1 > 0) || (value_2 > value_1 && value_1 < 0)) {
            T tmp_1 = value_1 * split_factor;
            // Equivalent to std::floor which is not yet constexpr.
            tmp_1 = static_cast<T>(static_cast<long long>(tmp_1));
            tmp_1 = tmp_1 / split_factor;
            value_2 += (value_1 - tmp_1);
            value_1 = tmp_1;
        }
    }
    T value_1{};
    T value_2{};
};

template <typename CT>
static CT constexpr square_root(CT x) {
    using T = typename std::remove_const<CT>::type;
    T guess = x;
    T previous_guess;
    T diff;
    T two = 2.0;
    for (int8_t i = 32; i > 0; i--) {
        previous_guess = guess;
        guess = (guess + x / guess) / two;
        diff = guess - previous_guess;
    }

    return guess;
}

template <typename CT>
static consteval CT gauss_legendre_algorithm() {
    using T = typename std::remove_const<CT>::type;
    T two = 2.0l;
    T four = 4.0l;
    T a = 1.0l;
    T b = 1.0l / square_root(2.0l);
    T p = 1.0l;
    T t = 1.0l / 4.0l;
    for (int8_t i = 5; i > 0; i--) {
        T a_next = (a + b) / two;
        b = square_root(a * b);
        t = t - p * (a - a_next) * (a - a_next);
        p = two * p;
        a = a_next;
    }
    return ((a + b) * (a + b) / (four * t));
}

int main() {
    auto compare_pi = [](auto t) {
        using T = decltype(t);
        constexpr T std_pi = std::numbers::pi_v<T>;
        constexpr T my_pi = gauss_legendre_algorithm<T>();
        constexpr T my_more_precise_pi = gauss_legendre_algorithm<more_precise<T>>().get();
        static_assert(my_pi != std_pi);
        static_assert(my_more_precise_pi == std_pi);
        return 0;
    };

    auto iterate_types = [&](auto... types) { (compare_pi(types), ...); };

    using long_double = long double;
    iterate_types(float{}, double{}, long_double{});
    return 0;
}