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