GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 66.5% 175 / 0 / 263
Functions: 87.5% 63 / 0 / 72
Branches: 86.2% 293 / 35 / 375
modnum.hpp
Line Branch Exec Source
1 #pragma once
2
3 #include <cassert>
4 #include <iostream>
5 #include <cstdint>
6 #include <concepts>
7 #include <limits>
8 #include <type_traits>
9
10 10187695 template <typename T> T mod_inv_in_range(T a, T m) {
11 // assert(0 <= a && a < m);
12 10187695 T x = a, y = m;
13 // abs coeff of a in x and y (they're always opposite sign)
14 10187695 T vx = 1, vy = 0;
15 10187695 bool swap = false;
16
8/8
unsigned int mod_inv_in_range<unsigned int>(unsigned int, unsigned int):
✓ Branch 4 → 3 taken 30939796 times.
✓ Branch 4 → 5 taken 10162400 times.
✓ Branch 10 → 7 taken 79818 times.
✓ Branch 10 → 11 taken 15804 times.
unsigned long mod_inv_in_range<unsigned long>(unsigned long, unsigned long):
✓ Branch 4 → 3 taken 14270 times.
✓ Branch 4 → 5 taken 628 times.
✓ Branch 10 → 7 taken 63092 times.
✓ Branch 10 → 11 taken 8863 times.
 41284671 while (x) {
17 31096976 T k = y / x;
18 31096976 y %= x;
19 31096976 vy += k * vx;
20 31096976 std::swap(x, y);
21 31096976 std::swap(vx, vy);
22 31096976 swap ^= 1;
23 }
24 10187695 assert(y == 1);
25
8/8
unsigned int mod_inv_in_range<unsigned int>(unsigned int, unsigned int):
✓ Branch 7 → 8 taken 9495358 times.
✓ Branch 7 → 9 taken 667042 times.
✓ Branch 13 → 14 taken 5426 times.
✓ Branch 13 → 15 taken 10378 times.
unsigned long mod_inv_in_range<unsigned long>(unsigned long, unsigned long):
✓ Branch 7 → 8 taken 476 times.
✓ Branch 7 → 9 taken 152 times.
✓ Branch 13 → 14 taken 5368 times.
✓ Branch 13 → 15 taken 3495 times.
 10187695 return swap ? vy : m - vy;
26 }
27
28 template <typename T> struct extended_gcd_result {
29 T gcd;
30 T coeff_a, coeff_b;
31 };
32 2018 template <typename T> extended_gcd_result<T> extended_gcd(T a, T b) {
33 2018 T x = a, y = b;
34 // coeff of a and b in x and y
35 2018 T ax = 1, ay = 0;
36 2018 T bx = 0, by = 1;
37
2/2
✓ Branch 14 → 9 taken 6604 times.
✓ Branch 14 → 15 taken 2018 times.
 8622 while (x) {
38
1/2
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 6604 times.
 6604 T k = y / x;
39 6604 y %= x;
40 6604 ay -= k * ax;
41 6604 by -= k * bx;
42 6604 std::swap(x, y);
43 6604 std::swap(ax, ay);
44 6604 std::swap(bx, by);
45 }
46 2018 return {y, ay, by};
47 }
48
49 ✗ template <typename T> T mod_inv(T a, T m) {
50 ✗ a %= m;
51 ✗ a = a < 0 ? a + m : a;
52 ✗ return mod_inv_in_range(a, m);
53 }
54
55 // Derives the boilerplate operator surface of a number type from its compound
56 // ops, ==, neg(), and inv().
57 // Bodies are only instantiated on use, so a type may omit some of the
58 // underlying pieces if the corresponding derived ops are never called.
59 template <typename Self>
60 struct num_ops {
61 ✗ Self operator+ () const { return static_cast<const Self&>(*this); }
62 32391534 Self operator- () const { return static_cast<const Self&>(*this).neg(); }
63
64 ✗ friend Self operator ++ (Self& a, int) { Self r = a; ++a; return r; }
65 ✗ friend Self operator -- (Self& a, int) { Self r = a; --a; return r; }
66 22912491636 friend Self operator + (const Self& a, const Self& b) { return Self(a) += b; }
67 22591061507 friend Self operator - (const Self& a, const Self& b) { return Self(a) -= b; }
68 13700966883 friend Self operator * (const Self& a, const Self& b) { return Self(a) *= b; }
69 13626853 friend Self operator / (const Self& a, const Self& b) { return Self(a) /= b; }
70
71 8201 friend bool operator != (const Self& a, const Self& b) { return !(a == b); }
72
73 ✗ friend Self neg(const Self& a) { return a.neg(); }
74 3374014 friend Self inv(const Self& a) { return a.inv(); }
75 };
76
77 // Storage and arithmetic for numbers mod Self::MOD, as a reduced
78 // representative v in [0, MOD) of unsigned type V.
79 // The type provides static MOD (of type V), reduce (value -> representative),
80 // and *=;
81 // everything else is derived here, valid for any MOD up to V's full range
82 // (sums and differences are tracked mod 2^bits, so no headroom is needed).
83 // Hooks may be overridden in the type's own body (e.g. a faster += / -=).
84 template <typename Self, typename V>
85 struct mod_ops : num_ops<Self> {
86 static_assert(std::unsigned_integral<V>);
87 V v;
88
89 struct is_reduced_tag {};
90
91 4703151036 mod_ops() : v(0) {}
92 42579229 mod_ops(V v_, is_reduced_tag) : v(v_) { assert(v < Self::MOD); }
93 1199606183 template <std::integral I> mod_ops(I x) : v(Self::reduce(x)) {}
94
95 85133791 static Self from_reduced(V v) { return Self(v, is_reduced_tag{}); }
96
97 // A negative value reduces via its nonnegative complement: x = -1 - ~x.
98 801517218 static V reduce(std::signed_integral auto x) {
99 using U = std::make_unsigned_t<decltype(x)>;
100
3/6
None:
✓ Branch 2 → 3 taken 62110 times.
✓ Branch 2 → 4 taken 406673 times.
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 6169160 times.
unsigned long mod_ops<mod_goldilocks, unsigned long>::reduce<int>(int):
✗ Branch 2 → 3 not taken.
✗ Branch 2 → 5 not taken.
 1192064713 return x < 0 ? V(Self::MOD - 1 - Self::reduce(U(~x))) : Self::reduce(U(x));
101 }
102
103 200688 explicit operator V() const { return v; }
104 44197652 std::make_signed_t<V> balanced() const {
105
6/6
None:
✓ Branch 2 → 3 taken 11556440 times.
✓ Branch 2 → 4 taken 10521129 times.
✓ Branch 6 → 7 taken 11556440 times.
✓ Branch 6 → 8 taken 10521129 times.
mod_ops<modnum<1000000007>, unsigned int>::balanced() const:
✓ Branch 6 → 7 taken 22904 times.
✓ Branch 6 → 9 taken 19610 times.
 44197652 return std::make_signed_t<V>(Self::MOD-v > v ? v : v - Self::MOD);
106 }
107
108
10/16
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 42 times.
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 20 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 42 times.
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 20 times.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 60 times.
✓ Branch 9 → 10 taken 376543 times.
✓ Branch 9 → 11 taken 7894 times.
✓ Branch 13 → 14 taken 3971648 times.
✓ Branch 13 → 15 taken 610 times.
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 29 times.
 4556122 friend bool operator == (const Self& a, const Self& b) { return a.v == b.v; }
109
16/16
✓ Branch 11 → 12 taken 16696892 times.
✓ Branch 12 → 13 taken 30644307 times.
✓ Branch 13 → 14 taken 41198909 times.
✓ Branch 14 → 15 taken 27251494 times.
✓ Branch 16 → 17 taken 1878056 times.
✓ Branch 17 → 18 taken 22302709 times.
✓ Branch 18 → 19 taken 43193692 times.
✓ Branch 19 → 20 taken 66934098 times.
✓ Branch 20 → 21 taken 44165059 times.
✓ Branch 22 → 23 taken 644060 times.
✓ Branch 23 → 24 taken 644060 times.
✓ Branch 37 → 38 taken 1919469 times.
✓ Branch 38 → 39 taken 3278649 times.
✓ Branch 39 → 40 taken 1359180 times.
✓ Branch 42 → 43 taken 1359630 times.
✓ Branch 43 → 44 taken 1359630 times.
 152496947 friend std::ostream& operator << (std::ostream& out, const Self& n) { return out << n.v; }
110 202798445 friend std::istream& operator >> (std::istream& in, Self& n) { int64_t v_; in >> v_; n = Self(v_); return in; }
111
112 ✗ Self& operator ++ () {
113 ✗ ++v;
114 ✗ if (v == Self::MOD) v = 0;
115 ✗ return self();
116 }
117 ✗ Self& operator -- () {
118 ✗ if (v == 0) v = Self::MOD;
119 ✗ --v;
120 ✗ return self();
121 }
122
17/18
None:
✓ Branch 37 → 38 taken 1919469 times.
mod_ops<modnum<2013265921>, unsigned int>::operator+=(modnum<2013265921> const&):
✓ Branch 8 → 9 taken 202208976 times.
✓ Branch 8 → 10 taken 290893518 times.
✓ Branch 21 → 22 taken 17321299 times.
✓ Branch 21 → 25 taken 4756214 times.
mod_ops<modnum<998244353>, unsigned int>::operator+=(modnum<998244353> const&):
✓ Branch 8 → 9 taken 2926750909 times.
✓ Branch 8 → 10 taken 4944177003 times.
✓ Branch 11 → 12 taken 10199479 times.
✓ Branch 11 → 13 taken 10195153 times.
✓ Branch 17 → 18 taken 15 times.
✓ Branch 17 → 36 taken 15377153 times.
✓ Branch 31 → 32 taken 15 times.
✓ Branch 60 → 61 taken 116650992 times.
✓ Branch 60 → 62 taken 2417711716 times.
✓ Branch 71 → 72 taken 5116137 times.
✓ Branch 71 → 73 taken 6587158 times.
✗ Branch 85 → 86 not taken.
✓ Branch 85 → 87 taken 11703295 times.
 19890327890 Self& operator += (const Self& o) { v = Self::sub_mod_raw(v, Self::MOD - o.v); return self(); }
123
27/27
✓ Branch 5 → 6 taken 27256182 times.
✓ Branch 5 → 7 taken 19702212 times.
✓ Branch 6 → 7 taken 2 times.
✓ Branch 6 → 8 taken 11 times.
✓ Branch 8 → 9 taken 3082738973 times.
✓ Branch 8 → 10 taken 2612144891 times.
✓ Branch 11 → 12 taken 16255453 times.
✓ Branch 11 → 13 taken 10308309 times.
✓ Branch 16 → 17 taken 2361333077 times.
✓ Branch 16 → 18 taken 2730847960 times.
✓ Branch 19 → 20 taken 2358256427 times.
✓ Branch 19 → 21 taken 2733924610 times.
✓ Branch 22 → 23 taken 2355893560 times.
✓ Branch 22 → 24 taken 2736287477 times.
✓ Branch 25 → 26 taken 2736572239 times.
✓ Branch 25 → 27 taken 2360353732 times.
✓ Branch 25 → 29 taken 17332579 times.
✓ Branch 42 → 43 taken 20396714 times.
✓ Branch 42 → 44 taken 26561680 times.
✓ Branch 45 → 46 taken 25868191 times.
✓ Branch 45 → 47 taken 21090203 times.
✓ Branch 54 → 55 taken 1180798987 times.
✓ Branch 54 → 56 taken 1353563721 times.
✓ Branch 57 → 58 taken 1363580696 times.
✓ Branch 57 → 59 taken 1170782012 times.
✓ Branch 70 → 71 taken 417264 times.
✓ Branch 70 → 72 taken 23046762 times.
 12553159019 Self& operator -= (const Self& o) { v = Self::sub_mod_raw(v, o.v); return self(); }
124 6814046 Self& operator /= (const Self& o) { return self() *= o.inv(); }
125
126 // Returns a - b mod MOD, for b in [0, MOD]; wraparound detects the underflow.
127
48/50
✓ Branch 2 → 3 taken 6011068771 times.
✓ Branch 2 → 4 taken 5255072932 times.
✓ Branch 3 → 4 taken 161157 times.
✓ Branch 3 → 5 taken 1486935 times.
✓ Branch 4 → 5 taken 23436729 times.
✓ Branch 4 → 6 taken 23438267 times.
✓ Branch 5 → 6 taken 7833249040 times.
✓ Branch 5 → 7 taken 6178706836 times.
✓ Branch 6 → 7 taken 77774034 times.
✓ Branch 6 → 8 taken 77396225 times.
✓ Branch 7 → 8 taken 12141324 times.
✓ Branch 7 → 9 taken 12141816 times.
✓ Branch 8 → 9 taken 6206173833 times.
✓ Branch 8 → 10 taken 7842990242 times.
✓ Branch 9 → 10 taken 46874996 times.
✗ Branch 9 → 11 not taken.
✓ Branch 11 → 12 taken 10199479 times.
✓ Branch 11 → 13 taken 10195153 times.
✓ Branch 13 → 14 taken 2364098538 times.
✓ Branch 13 → 15 taken 2728082499 times.
✓ Branch 14 → 15 taken 7454549 times.
✓ Branch 14 → 16 taken 7922619 times.
✓ Branch 16 → 17 taken 2366906321 times.
✓ Branch 16 → 18 taken 2738087637 times.
✓ Branch 19 → 20 taken 2358256427 times.
✓ Branch 19 → 21 taken 2733924610 times.
✓ Branch 22 → 23 taken 2357413680 times.
✓ Branch 22 → 24 taken 2752088656 times.
✓ Branch 25 → 26 taken 2731827305 times.
✓ Branch 25 → 27 taken 2360353732 times.
✓ Branch 26 → 27 taken 443749 times.
✓ Branch 26 → 28 taken 4301185 times.
✓ Branch 28 → 29 taken 10 times.
✓ Branch 28 → 30 taken 5 times.
✓ Branch 34 → 35 taken 8671485 times.
✓ Branch 34 → 36 taken 14792541 times.
✓ Branch 39 → 40 taken 20507933 times.
✓ Branch 39 → 41 taken 26450461 times.
✓ Branch 42 → 43 taken 20396714 times.
✓ Branch 42 → 44 taken 26561680 times.
✓ Branch 51 → 52 taken 1181502678 times.
✓ Branch 51 → 53 taken 1352860030 times.
✗ Branch 54 → 55 not taken.
✓ Branch 54 → 56 taken 60 times.
✓ Branch 67 → 68 taken 8388294 times.
✓ Branch 67 → 69 taken 15075732 times.
✓ Branch 68 → 69 taken 5116137 times.
✓ Branch 68 → 70 taken 6587158 times.
✓ Branch 82 → 83 taken 5116137 times.
✓ Branch 82 → 84 taken 6587158 times.
 14175795665 static V sub_mod_raw(V a, V b) { return a < b ? a - b + Self::MOD : a - b; }
128
129
6/8
mod_ops<mod_goldilocks, unsigned long>::neg() const:
✓ Branch 4 → 5 taken 340 times.
✗ Branch 4 → 8 not taken.
mod_ops<modnum<1000000007>, unsigned int>::neg() const:
✓ Branch 4 → 5 taken 2085 times.
✗ Branch 4 → 8 not taken.
mod_ops<modnum<998244353>, unsigned int>::neg() const:
✓ Branch 2 → 3 taken 28171870 times.
✓ Branch 2 → 4 taken 4214777 times.
✓ Branch 4 → 5 taken 2383 times.
✓ Branch 4 → 8 taken 79 times.
 32391534 Self neg() const { return from_reduced(v ? Self::MOD - v : 0); }
130 10187695 Self inv() const { return from_reduced(mod_inv_in_range(v, Self::MOD)); }
131
132 private:
133 8623364 Self& self() { return static_cast<Self&>(*this); }
134 };
135
136
13/15
✓ Branch 2 → 3 taken 1326178 times.
✗ Branch 2 → 5 not taken.
✓ Branch 2 → 7 taken 998 times.
✓ Branch 2 → 10 taken 1325939 times.
✓ Branch 6 → 7 taken 1330720 times.
✓ Branch 6 → 8 taken 76 times.
✓ Branch 6 → 18 taken 98 times.
✓ Branch 8 → 9 taken 146 times.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 1648003 times.
✓ Branch 22 → 23 taken 20 times.
✓ Branch 28 → 29 taken 9 times.
✓ Branch 31 → 32 taken 1919469 times.
✓ Branch 51 → 52 taken 1181502678 times.
✓ Branch 51 → 53 taken 1352860030 times.
 4629824628 template <auto MOD_> struct modnum : mod_ops<modnum<MOD_>, std::make_unsigned_t<decltype(MOD_)>> {
137 using Self = modnum;
138 static_assert(MOD_ > 0, "MOD must be positive");
139 using V = std::make_unsigned_t<decltype(MOD_)>;
140 static constexpr V MOD = V(MOD_);
141
142 using base = mod_ops<modnum, V>;
143
43/52
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 6169202 times.
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 202798465 times.
✓ Branch 4 → 5 taken 2644805 times.
✓ Branch 4 → 6 taken 42 times.
✓ Branch 5 → 6 taken 28 times.
✓ Branch 5 → 7 taken 6169180 times.
✓ Branch 6 → 7 taken 1325960 times.
✓ Branch 6 → 8 taken 52582799 times.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 3304743 times.
✓ Branch 8 → 9 taken 8 times.
✓ Branch 9 → 10 taken 376563 times.
✓ Branch 9 → 11 taken 7894 times.
✓ Branch 10 → 11 taken 43 times.
✓ Branch 11 → 12 taken 21097472 times.
✓ Branch 11 → 13 taken 16357285 times.
✓ Branch 12 → 13 taken 1696 times.
✓ Branch 12 → 14 taken 644060 times.
✓ Branch 13 → 14 taken 3972925 times.
✓ Branch 13 → 15 taken 610 times.
✓ Branch 14 → 15 taken 7454549 times.
✓ Branch 14 → 16 taken 7922619 times.
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 22077513 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 590 times.
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 10644437 times.
✓ Branch 18 → 19 taken 10618909 times.
✓ Branch 18 → 20 taken 11458604 times.
✓ Branch 21 → 22 taken 20 times.
✓ Branch 25 → 26 taken 79 times.
✓ Branch 25 → 27 taken 15 times.
✓ Branch 28 → 29 taken 10 times.
✓ Branch 28 → 30 taken 5 times.
✓ Branch 30 → 31 taken 59 times.
✓ Branch 31 → 32 taken 17 times.
✓ Branch 32 → 33 taken 59 times.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 6587187 times.
✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 29 times.
✓ Branch 65 → 66 taken 6587158 times.
✓ Branch 65 → 67 taken 5116137 times.
✓ Branch 68 → 69 taken 5116137 times.
✓ Branch 68 → 70 taken 6587158 times.
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 29 times.
✓ Branch 82 → 83 taken 5116137 times.
✓ Branch 82 → 84 taken 6587158 times.
 798096689 using base::base;
144 using base::v;
145 using base::reduce;
146
147
1/1
✓ Branch 24 → 25 taken 17 times.
 404506783 static V reduce(std::unsigned_integral auto x) { return V(x % MOD); }
148
149 22155694 explicit operator std::make_signed_t<V>() const
150 requires (MOD <= V(std::numeric_limits<std::make_signed_t<V>>::max()))
151 {
152
3/4
✓ Branch 3 → 4 taken 220 times.
✗ Branch 3 → 21 not taken.
✓ Branch 106 → 107 taken 16916 times.
✓ Branch 106 → 108 taken 21928 times.
 22155474 return std::make_signed_t<V>(v);
153 }
154
155 14742095568 Self& operator *= (const Self& o) {
156
38/41
✓ Branch 2 → 3 taken 440546128 times.
✓ Branch 2 → 4 taken 191711364 times.
✓ Branch 5 → 6 taken 2587949461 times.
✓ Branch 5 → 7 taken 3133498165 times.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 10644437 times.
✓ Branch 8 → 9 taken 15530640 times.
✓ Branch 8 → 10 taken 25673124 times.
✓ Branch 11 → 12 taken 16255453 times.
✓ Branch 11 → 13 taken 10308309 times.
✓ Branch 12 → 13 taken 321979 times.
✓ Branch 12 → 14 taken 322072 times.
✓ Branch 13 → 14 taken 1487405 times.
✓ Branch 13 → 15 taken 160526 times.
✓ Branch 15 → 16 taken 494484 times.
✓ Branch 15 → 17 taken 22687400 times.
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 10644437 times.
✓ Branch 18 → 19 taken 10618909 times.
✓ Branch 18 → 20 taken 11458604 times.
✓ Branch 25 → 26 taken 2731827305 times.
✓ Branch 25 → 27 taken 2360353732 times.
✓ Branch 28 → 29 taken 446296308 times.
✓ Branch 28 → 30 taken 448311263 times.
✓ Branch 34 → 35 taken 960158 times.
✓ Branch 34 → 36 taken 959311 times.
✓ Branch 45 → 46 taken 25868191 times.
✓ Branch 45 → 47 taken 21090203 times.
✓ Branch 49 → 50 taken 3293587 times.
✓ Branch 49 → 52 taken 3293600 times.
✓ Branch 53 → 54 taken 60 times.
✓ Branch 54 → 55 taken 1180798987 times.
✓ Branch 54 → 56 taken 1353563721 times.
✓ Branch 57 → 58 taken 1363580696 times.
✓ Branch 57 → 59 taken 1170782012 times.
✓ Branch 71 → 72 taken 5116137 times.
✓ Branch 71 → 73 taken 6587158 times.
✓ Branch 79 → 80 taken 6587158 times.
✓ Branch 79 → 81 taken 5116137 times.
✗ Branch 85 → 86 not taken.
✓ Branch 85 → 87 taken 11703295 times.
 14297827921 if constexpr (sizeof(V) <= 4) v = V(uint64_t(v) * o.v % MOD);
157 else v = V(__uint128_t(v) * o.v % MOD);
158 9110614 return *this;
159 }
160 };
161
162 struct mod_goldilocks : mod_ops<mod_goldilocks, uint64_t> {
163 using Self = mod_goldilocks;
164 static constexpr uint64_t MOD = 0xffffffff00000001ull;
165 static constexpr uint64_t EPS = -MOD;
166 // We have 2^32 is a primitive 6th root of unity.
167 // Note that omega_8 + omega_8^7 == 2^24 - 2^72 == sqrt(2)
168 // We'll pick the root so that 2^24 - 2^72 is our primitive 384th root of unity.
169 static constexpr uint64_t PRIMITIVE_ROOT = 2717;
170
171 using base = mod_ops<mod_goldilocks, uint64_t>;
172
6/7
✓ Branch 4 → 5 taken 96 times.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 48 times.
✓ Branch 8 → 9 taken 11556440 times.
✓ Branch 8 → 10 taken 10521129 times.
✓ Branch 11 → 12 taken 11556440 times.
✓ Branch 11 → 13 taken 10521129 times.
 44294161 using base::base;
173 using base::reduce;
174
3/3
✓ Branch 6 → 7 taken 166 times.
✓ Branch 6 → 8 taken 16 times.
✓ Branch 6 → 18 taken 20 times.
 74213843 mod_goldilocks() = default;
175 ✗ mod_goldilocks(__int128_t a) : base(a < 0 ? uint64_t(MOD - 1 - __uint128_t(~a) % MOD) : uint64_t(__uint128_t(a) % MOD), is_reduced_tag{}) {}
176 ✗ mod_goldilocks(__uint128_t a) : base(uint64_t(a % MOD), is_reduced_tag{}) {}
177
178 // Avoids the division: any uint64_t is within MOD of reduced.
179 85523 static uint64_t reduce(std::unsigned_integral auto x) {
180 static_assert(sizeof(x) <= 8);
181 22192939 uint64_t a = x;
182
1/6
None:
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 7172 times.
unsigned long mod_goldilocks::reduce<unsigned int>(unsigned int):
✗ Branch 2 → 3 not taken.
✗ Branch 2 → 4 not taken.
unsigned long mod_goldilocks::reduce<unsigned long>(unsigned long):
✗ Branch 2 → 3 not taken.
✗ Branch 2 → 4 not taken.
 7172 return a >= MOD ? a - MOD : a;
183 }
184
185 // returns a-b, assuming -MOD <= a-b, e.g. b <= MOD
186 1574704951 static uint64_t sub_mod_raw(uint64_t a, uint64_t b) {
187 #if defined(__x86_64__)
188 // TODO: We could try to write this using intrinsics, but GCC sometimes produces the wrong code.
189 1574704951 uint64_t res_wrapped = a;
190 1574704951 uint64_t adjustment = b;
191 1574704951 asm (
192 // AT&T syntax: SRC DST
193 "sub %[y], %[x]\n\t"
194 // Trick from plonky2 implementation:
195 // After sub, flag CF is set iff we underflowed. We want to correct by EPS == 2^32 - 1 iff C is set.
196 // sbb (subtract with borrow) computes DST <- DST - SRC - CF
197 // Thus, we can use the 32-bit form of sbb on a dummy register to load CF ? EPS : 0.
198 // Here, we'll just reuse the original register holding b.
199 "sbb %k[y], %k[y]\n\t"
200 : [x] "+r"(res_wrapped),
201 [y] "+r"(adjustment)
202 :
203 : "cc"
204 );
205 #else
206 uint64_t res_wrapped = a - b;
207 uint64_t adjustment = (res_wrapped > a) ? EPS : 0;
208 #endif
209 740657560 return res_wrapped - adjustment;
210 }
211
212 // Reduce lo + 2^64 * mi + 2^96 * hi, where hi <= MOD
213 834047391 static uint64_t reduce_u160_raw(uint64_t lo, uint32_t mi, uint64_t hi) {
214 // result = lo - hi + EPS * mi
215 // 0 <= lo <= 2^64 - 1 = MOD + EPS - 1
216 // 0 <= EPS * mi <= (2^32 - 1) * EPS = MOD - 1 - EPS
217 // 0 <= hi <= MOD
218 // -MOD <= lo - hi + EPS * mi <= 2*MOD-2
219 // so we do have some leeway
220 834047391 return sub_mod_raw(sub_mod_raw(lo, hi), MOD-(uint64_t(mi)<<32)+mi);
221 }
222
223 834047391 static uint64_t reduce_u128_raw(__uint128_t v) {
224 834047391 uint64_t hi = uint64_t(v >> 64);
225 834047391 uint64_t lo = uint64_t(v);
226 834047391 uint32_t hi_hi = uint32_t(hi >> 32);
227 834047391 uint32_t hi_lo = uint32_t(hi);
228 834047391 return reduce_u160_raw(lo, hi_lo, hi_hi);
229 }
230
231 821697866 Self& operator *= (Self o) {
232 821697866 v = reduce_u128_raw(__uint128_t(v) * __uint128_t(o.v));
233
1/2
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 22077513 times.
 821652273 return *this;
234 }
235 };
236
237 4746 template <typename T> T power(T a, long long b) {
238 4746 assert(b >= 0);
239
22/26
mod_goldilocks power<mod_goldilocks>(mod_goldilocks, long long):
✓ Branch 4 → 5 taken 17024 times.
✓ Branch 4 → 6 taken 11545 times.
✗ Branch 6 → 7 not taken.
✗ Branch 6 → 8 not taken.
✓ Branch 7 → 4 taken 28569 times.
✓ Branch 7 → 8 taken 532 times.
✓ Branch 8 → 9 taken 5632 times.
✓ Branch 8 → 10 taken 4834 times.
✗ Branch 9 → 6 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 11 → 8 taken 10466 times.
✓ Branch 11 → 12 taken 176 times.
modnum<2013265921> power<modnum<2013265921> >(modnum<2013265921>, long long):
✓ Branch 5 → 6 taken 2128 times.
✓ Branch 5 → 7 taken 8885 times.
✓ Branch 8 → 5 taken 11013 times.
✓ Branch 8 → 9 taken 904 times.
✓ Branch 8 → 10 taken 2092 times.
✓ Branch 11 → 8 taken 2464 times.
✓ Branch 11 → 12 taken 93 times.
modnum<998244353> power<modnum<998244353> >(modnum<998244353>, long long):
✓ Branch 5 → 6 taken 19551 times.
✓ Branch 5 → 7 taken 46301 times.
✓ Branch 8 → 5 taken 65852 times.
✓ Branch 8 → 9 taken 4338 times.
✓ Branch 8 → 10 taken 3580 times.
✓ Branch 11 → 8 taken 4695 times.
✓ Branch 11 → 12 taken 190 times.
 128264 T r = 1; while (b) { if (b & 1) r *= a; b >>= 1; a *= a; } return r;
240 }
241
242 template <typename U, typename V> struct pairnum : num_ops<pairnum<U, V>> {
243 using Self = pairnum;
244 U u;
245 V v;
246
247 ✗ pairnum() : u(0), v(0) {}
248 ✗ pairnum(long long val) : u(val), v(val) {}
249 ✗ pairnum(const U& u_, const V& v_) : u(u_), v(v_) {}
250
251 ✗ friend std::ostream& operator << (std::ostream& out, const Self& n) { return out << '(' << n.u << ',' << ' ' << n.v << ')'; }
252 ✗ friend std::istream& operator >> (std::istream& in, Self& n) { long long val; in >> val; n = Self(val); return in; }
253
254 ✗ friend bool operator == (const Self& a, const Self& b) { return a.u == b.u && a.v == b.v; }
255
256 ✗ Self inv() const {
257 ✗ return Self(u.inv(), v.inv());
258 }
259 ✗ Self neg() const {
260 ✗ return Self(u.neg(), v.neg());
261 }
262
263 ✗ Self& operator ++ () {
264 ✗ ++u, ++v;
265 ✗ return *this;
266 }
267 ✗ Self& operator -- () {
268 ✗ --u, --v;
269 ✗ return *this;
270 }
271
272 ✗ Self& operator += (const Self& o) {
273 ✗ u += o.u;
274 ✗ v += o.v;
275 ✗ return *this;
276 }
277 ✗ Self& operator -= (const Self& o) {
278 ✗ u -= o.u;
279 ✗ v -= o.v;
280 ✗ return *this;
281 }
282 ✗ Self& operator *= (const Self& o) {
283 ✗ u *= o.u;
284 ✗ v *= o.v;
285 ✗ return *this;
286 }
287 ✗ Self& operator /= (const Self& o) {
288 ✗ u /= o.u;
289 ✗ v /= o.v;
290 ✗ return *this;
291 }
292 };
293
294 template <typename tag> struct dynamic_modnum : mod_ops<dynamic_modnum<tag>, uint32_t> {
295 using Self = dynamic_modnum;
296 private:
297 inline static uint32_t MOD_ = 0;
298 inline static uint64_t BARRETT_M = 0;
299
300 public:
301 // Make only the const-reference public, to force the use of set_mod
302 static constexpr uint32_t const& MOD = MOD_;
303
304 using base = mod_ops<dynamic_modnum, uint32_t>;
305 using base::base;
306 using base::v;
307 using base::reduce;
308
309 // Barret reduction taken from KACTL:
310 /**
311 * Author: Simon Lindholm
312 * Date: 2020-05-30
313 * License: CC0
314 * Source: https://en.wikipedia.org/wiki/Barrett_reduction
315 * Description: Compute $a \% b$ about 5 times faster than usual, where $b$ is constant but not known at compile time.
316 * Returns a value congruent to $a \pmod b$ in the range $[0, 2b)$.
317 * Status: proven correct, stress-tested
318 * Measured as having 4 times lower latency, and 8 times higher throughput, see stress-test.
319 * Details:
320 * More precisely, it can be proven that the result equals 0 only if $a = 0$,
321 * and otherwise lies in $[1, (1 + a/2^64) * b)$.
322 */
323 ✗ static void set_mod(int mod) {
324 ✗ assert(mod > 0);
325 ✗ MOD_ = uint32_t(mod);
326 ✗ BARRETT_M = (uint64_t(-1) / MOD);
327 }
328 ✗ static uint32_t barrett_reduce_partial(uint64_t a) {
329 ✗ return uint32_t(a - uint64_t((__uint128_t(BARRETT_M) * a) >> 64) * MOD);
330 }
331 ✗ static uint32_t barrett_reduce(uint64_t a) {
332 ✗ int32_t res = int32_t(barrett_reduce_partial(a) - MOD);
333 ✗ return uint32_t((res < 0) ? res + int32_t(MOD) : res);
334 }
335
336 struct mod_reader {
337 ✗ friend std::istream& operator >> (std::istream& i, mod_reader) {
338 ✗ int mod; i >> mod;
339 ✗ Self::set_mod(mod);
340 ✗ return i;
341 }
342 };
343 ✗ static mod_reader MOD_READER() {
344 ✗ return mod_reader();
345 }
346
347 ✗ static uint32_t reduce(std::unsigned_integral auto x) {
348 static_assert(sizeof(x) <= 8);
349 ✗ return barrett_reduce(x);
350 }
351
352 ✗ explicit operator int() const { return int(v); }
353
354 ✗ Self& operator *= (const Self& o) {
355 ✗ v = barrett_reduce(uint64_t(v) * o.v);
356 ✗ return *this;
357 }
358 };
359
360 template <typename T> struct mod_constraint {
361 T v, mod;
362
363 2127 friend mod_constraint operator & (mod_constraint a, mod_constraint b) {
364
2/2
✓ Branch 2 → 3 taken 1036 times.
✓ Branch 2 → 4 taken 1091 times.
 2127 if (a.mod < b.mod) std::swap(a, b);
365
2/2
✓ Branch 4 → 5 taken 109 times.
✓ Branch 4 → 6 taken 2018 times.
 2127 if (b.mod == 1) return a;
366
367 2018 extended_gcd_result<T> egcd = extended_gcd<T>(a.mod, b.mod);
368 2018 assert(a.v % egcd.gcd == b.v % egcd.gcd);
369
370
1/2
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 2018 times.
 2018 T extra = b.v - a.v % b.mod;
371
1/2
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 2018 times.
 2018 extra /= egcd.gcd;
372
373 2018 extra *= egcd.coeff_a;
374
2/4
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 2018 times.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 2018 times.
 2018 extra %= b.mod / egcd.gcd;
375
2/2
✓ Branch 21 → 22 taken 680 times.
✓ Branch 21 → 23 taken 1338 times.
 2018 extra += (extra < 0) ? b.mod / egcd.gcd : 0;
376
377 return mod_constraint{
378 2018 a.v + extra * a.mod,
379 2018 a.mod * (b.mod / egcd.gcd)
380 2018 };
381 }
382 };
383