GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 93
Functions: -% 0 / 0 / 0
Branches: -% 0 / 0 / 0
num/quaternion_hurwitz.hpp
Line Branch Exec Source
1 #pragma once
2
3 #include <utility>
4 #include <array>
5 #include <tuple>
6 #include <ostream>
7
8 namespace wala {
9
10 template <typename num = int>
11 struct hurwitz_quaternion {
12 // we store the doubled quaternion
13 num s,x,y,z;
14 ✗ hurwitz_quaternion() : s(0), x(0), y(0), z(0) {}
15 ✗ hurwitz_quaternion(num v) : s(2*v), x(0), y(0), z(0) {}
16 ✗ hurwitz_quaternion(num s_, num x_, num y_, num z_) : s(2*s_), x(2*x_), y(2*y_), z(2*z_) {}
17 struct doubled_coords_tag {};
18 ✗ hurwitz_quaternion(doubled_coords_tag, num s_, num x_, num y_, num z_) : s(s_), x(x_), y(y_), z(z_) {
19 ✗ assert((s & 1) == (x & 1) && (s & 1) == (y & 1) && (s & 1) == (z & 1));
20 }
21 ✗ friend std::ostream& operator << (std::ostream& o, const hurwitz_quaternion& q) {
22 ✗ o << double(q.s)/2;
23 {
24 ✗ std::ios_base::fmtflags f(o.flags());
25 ✗ o << std::showpos << double(q.x)/2 << "i" << double(q.y)/2 << "j" << double(q.z)/2 << "k";
26 ✗ o.flags(f);
27 }
28 ✗ return o;
29 }
30
31 ✗ explicit operator bool() const {
32 ✗ return s || x || y || z;
33 }
34
35 ✗ friend bool operator == (const hurwitz_quaternion& a, const hurwitz_quaternion& b) {
36 ✗ return std::tie(a.s,a.x,a.y,a.z) == std::tie(b.s,b.x,b.y,b.z);
37 }
38 ✗ friend bool operator != (const hurwitz_quaternion& a, const hurwitz_quaternion& b) { return !(a == b); }
39
40 ✗ num real_doubled() const {
41 ✗ return s;
42 }
43 ✗ num real() const {
44 ✗ assert(!(s & 1));
45 ✗ return s >> 1;
46 }
47 ✗ std::array<num, 3> imag_doubled() const {
48 ✗ return {x, y, z};
49 }
50 ✗ std::array<num, 3> imag() const {
51 ✗ assert(!(s & 1));
52 ✗ return {x>>1, y>>1, z>>1};
53 }
54 ✗ std::array<num, 4> coords_doubled() const {
55 ✗ return {s, x, y, z};
56 }
57 ✗ std::array<num, 4> coords() const {
58 ✗ assert(!(s & 1));
59 ✗ return {s>>1, x>>1, y>>1, z>>1};
60 }
61
62 ✗ friend num norm(const hurwitz_quaternion& q) {
63 ✗ return (q.s * q.s + q.x * q.x + q.y * q.y + q.z * q.z) >> 2;
64 }
65 ✗ friend hurwitz_quaternion conj(const hurwitz_quaternion& q) {
66 ✗ return hurwitz_quaternion(doubled_coords_tag{}, q.s, -q.x, -q.y, -q.z);
67 }
68
69 ✗ friend hurwitz_quaternion operator + (const hurwitz_quaternion& q) {
70 ✗ return hurwitz_quaternion(doubled_coords_tag{}, +q.s, +q.x, +q.y, +q.z);
71 }
72 ✗ friend hurwitz_quaternion operator - (const hurwitz_quaternion& q) {
73 ✗ return hurwitz_quaternion(doubled_coords_tag{}, -q.s, -q.x, -q.y, -q.z);
74 }
75
76 ✗ hurwitz_quaternion& operator += (const hurwitz_quaternion& o) {
77 ✗ s += o.s;
78 ✗ x += o.x;
79 ✗ y += o.y;
80 ✗ z += o.z;
81 ✗ return *this;
82 }
83 ✗ friend hurwitz_quaternion operator + (const hurwitz_quaternion& a, const hurwitz_quaternion& b) {
84 ✗ return hurwitz_quaternion(doubled_coords_tag{}, a.s + b.s, a.x + b.x, a.y + b.y, a.z + b.z);
85 }
86 ✗ hurwitz_quaternion& operator -= (const hurwitz_quaternion& o) {
87 ✗ s -= o.s;
88 ✗ x -= o.x;
89 ✗ y -= o.y;
90 ✗ z -= o.z;
91 ✗ return *this;
92 }
93 ✗ friend hurwitz_quaternion operator - (const hurwitz_quaternion& a, const hurwitz_quaternion& b) {
94 ✗ return hurwitz_quaternion(doubled_coords_tag{}, a.s - b.s, a.x - b.x, a.y - b.y, a.z - b.z);
95 }
96
97 ✗ friend hurwitz_quaternion operator * (const num& a, const hurwitz_quaternion& q) {
98 ✗ return hurwitz_quaternion(doubled_coords_tag{}, a*q.s, a*q.x, a*q.y, a*q.z);
99 }
100 ✗ friend hurwitz_quaternion operator * (const hurwitz_quaternion& q, const num& a) {
101 ✗ return hurwitz_quaternion(doubled_coords_tag{}, q.s*a, q.x*a, q.y*a, q.z*a);
102 }
103 ✗ hurwitz_quaternion& operator *= (const num& a) {
104 ✗ s *= a;
105 ✗ x *= a;
106 ✗ y *= a;
107 ✗ z *= a;
108 ✗ return *this;
109 }
110
111 ✗ friend hurwitz_quaternion operator * (const hurwitz_quaternion& a, const hurwitz_quaternion& b) {
112 ✗ return hurwitz_quaternion(
113 doubled_coords_tag{},
114 (a.s * b.s - a.x * b.x - a.y * b.y - a.z * b.z) >> 1,
115 (a.s * b.x + a.x * b.s + a.y * b.z - a.z * b.y) >> 1,
116 (a.s * b.y + a.y * b.s + a.z * b.x - a.x * b.z) >> 1,
117 (a.s * b.z + a.z * b.s + a.x * b.y - a.y * b.x) >> 1
118 );
119 }
120 ✗ hurwitz_quaternion& operator *= (const hurwitz_quaternion& o) {
121 ✗ return *this = *this * o;
122 }
123
124 struct div_t {
125 hurwitz_quaternion quot, rem;
126 };
127 // a = b * quot + rem
128 ✗ friend div_t right_div(const hurwitz_quaternion& a, const hurwitz_quaternion& b) {
129 ✗ hurwitz_quaternion numer = conj(b) * a;
130 ✗ num denom = norm(b);
131
132 ✗ auto floor_div = [](num u, num v) -> num {
133 ✗ if ((u^v) >= 0) {
134 ✗ return u/v;
135 } else {
136 ✗ auto res = std::div(u, v);
137 ✗ return res.quot - bool(res.rem);
138 }
139 };
140 ✗ num s = floor_div(numer.s, denom);
141 ✗ num x = floor_div(numer.x, denom);
142 ✗ num y = floor_div(numer.y, denom);
143 ✗ num z = floor_div(numer.z, denom);
144
145 ✗ hurwitz_quaternion q_odd(doubled_coords_tag{}, s | 1, x | 1, y | 1, z | 1);
146 ✗ hurwitz_quaternion r_odd = a - b * q_odd;
147 ✗ hurwitz_quaternion q_even(doubled_coords_tag{}, (s+1)&~num(1), (x+1)&~num(1), (y+1)&~num(1), (z+1)&~num(1));
148 ✗ hurwitz_quaternion r_even = a - b * q_even;
149 ✗ div_t res = norm(r_odd) < norm(r_even) ? div_t{q_odd, r_odd} : div_t{q_even, r_even};
150 ✗ assert(norm(res.rem) < norm(b));
151 ✗ return res;
152 }
153
154 // a = ga', b = gb'
155 ✗ friend hurwitz_quaternion right_gcd(hurwitz_quaternion a, hurwitz_quaternion b) {
156 ✗ while (a) {
157 ✗ b = right_div(b, a).rem;
158 ✗ std::swap(a, b);
159 }
160 ✗ return b;
161 }
162 };
163
164 } // namespace wala
165