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
|
#pragma once
#include <cmath>
#include <functional>
#include <vector>
#include <iostream>
typedef void i0;
typedef int8_t i8;
typedef uint8_t u8;
typedef int16_t i16;
typedef uint16_t u16;
typedef int32_t i32;
typedef uint32_t u32;
typedef int64_t i64;
typedef uint64_t u64;
typedef ssize_t isz;
typedef size_t usz;
typedef float f32;
typedef double f64;
#if __SIZEOF_LONG_DOUBLE__ == 16
typedef long double f128;
#else
typedef long double _f64;
#endif
#if defined(_MSC_VER) || !defined(__cpp_multidimensional_subscript) || __cplusplus < 202110L
#warning "can not use multidimentional subscript operator: falling back to `operator()`"
#undef MDSUBSCRIPT
#else
#define MDSUBSCRIPT
#endif
#ifdef NO_MDSUBSCRIPT
#undef MDSUBSCRIPT
#endif
#ifdef MDSUBSCRIPT
#define SUBSCR_OPN [
#define SUBSCR_CLS ]
#define SUBSCR_OPRTR operator[]
#else
#define SUBSCR_OPN (
#define SUBSCR_CLS )
#define SUBSCR_OPRTR operator()
#endif
template <typename T>
using function_t = std::function<T(const std::vector<T> &)>;
namespace prak {
/// stolen from [cppreference.com](https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon)
/// Compares 2 floating-point values up to `ulps` ULPS (units in the last place)
template <class T>
std::enable_if_t<not std::numeric_limits<T>::is_integer, bool>
fequal(T x, T y, std::size_t ulps = 1)
{
// Since `epsilon()` is the gap size (ULP, unit in the last place)
// of floating-point numbers in interval [1, 2), we can scale it to
// the gap size in interval [2^e, 2^{e+1}), where `e` is the exponent
// of `x` and `y`.
// If `x` and `y` have different gap sizes (which means they have
// different exponents), we take the smaller one. Taking the bigger
// one is also reasonable, I guess.
const T m = std::min(std::fabs(x), std::fabs(y));
// Subnormal numbers have fixed exponent, which is `min_exponent - 1`.
const int exp = m < std::numeric_limits<T>::min()
? std::numeric_limits<T>::min_exponent - 1
: std::ilogb(m);
// We consider `x` and `y` equal if the difference between them is
// within `n` ULPs.
return std::fabs(x - y) <= ulps * std::ldexp(std::numeric_limits<T>::epsilon(), exp);
}
/// prints a vector
template <typename T>
void printv(const T &vec) {
std::cout << "std::vector { ";
for (const auto &x : vec) {
std::cout << x << ' ';
}
std::cout << '}' << std::endl;
}
/// An allocator that aligns memory to 64 bytes. Needed for AVX instructions.
/// C++17 required for std::align_val_t
template <class T>
struct align_alloc : public std::allocator<T> {
constexpr T* allocate( std::size_t n ) {
return static_cast<T*>(::operator new(sizeof (T) * n, std::align_val_t{64}));
}
constexpr void deallocate( T* p, std::size_t n ) {
::operator delete(p, std::align_val_t{64});
}
};;
/// alias prak::vector that is the same as std::vector, but uses aligned allocator
template <typename T>
using vector = std::vector<T, align_alloc<T>>;
/// prak value / pair value: a value with an error
template <typename T>
struct pvalue { T val, err; };
template <typename T>
struct pvalue<T> operator*(const struct pvalue<T> &v, T a) {
return pvalue<T>{v.val * a, v.err * a};
}
template <typename T>
std::ostream &operator<<(std::ostream &os, const struct pvalue<T> &p) {
/* return os << "value {" << p.val << "±" << p.err << "}"; */
return os << p.val << "±" << p.err;
}
} // namespace prak
|
