src/nt/gl/gcd_mpi/lib.hpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #pragma once | ||
| 2 | |||
| 3 | #include "../../../math/ds/mpi/lib.hpp" | ||
| 4 | #include "../../../math/qpow/basic/lib.hpp" | ||
| 5 | #include "../gcd/lib.hpp" | ||
| 6 | |||
| 7 | namespace tifa_libs { | ||
| 8 | namespace gcd_mpi_impl_ { | ||
| 9 | // @return {x, y} s.t. a = 2^x 5^y | ||
| 10 | 13 | CEXP ptti shrink(u32 a) NE { | |
| 11 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
|
13 | assert(a > 0); |
| 12 | 13 | csint x = std::countr_zero(a); | |
| 13 | 13 | a >>= x; | |
| 14 | // clang-format off | ||
| 15 |
12/18✓ Branch 0 taken 9 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
|
13 | return {x, a == 1 ? 0 : a == 5 ? 1 : a == 25 ? 2 : a == 125 ? 3 : a == 625 ? 4 : a == 3125 ? 5 : a == 15625 ? 6 : a == 78125 ? 7 : a == 390625 ? 8 : 9}; |
| 16 | // clang-format on | ||
| 17 | } | ||
| 18 | 20 | CEXP ptti shrink(mpi& a) NE { | |
| 19 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 20 times.
|
20 | if (assert(!a.is_neg()); a.empty()) return {0, 0}; |
| 20 | 20 | ptti res{0, 0}; | |
| 21 | while (true) { | ||
| 22 | 33 | cu32 g = gcd(mpi::D, a[0]); | |
| 23 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 13 times.
|
33 | if (g == 1) break; |
| 24 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1 time.
|
13 | if (g != mpi::D) a *= mpi::D / g; |
| 25 | 13 | a.erase(begin(a)); | |
| 26 | 13 | const auto s = shrink(g); | |
| 27 | 13 | res.first += s.first, res.second += s.second; | |
| 28 | 13 | } | |
| 29 | 20 | return res; | |
| 30 | } | ||
| 31 | } // namespace gcd_mpi_impl_ | ||
| 32 | |||
| 33 | template <bool FAST = true> | ||
| 34 | 10 | CEXP mpi gcd_mpi(mpi a, mpi b) NE { | |
| 35 | 10 | a.set_neg(false), b.set_neg(false); | |
| 36 | if CEXP (FAST) | ||
| 37 |
1/2✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
|
10 | if (max(a.size(), b.size()) <= 2) return gcd(a.to_i64(), b.to_i64()); |
| 38 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
|
10 | if (a.empty()) return b; |
| 39 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
|
10 | if (b.empty()) return a; |
| 40 | 10 | cptti s = gcd_mpi_impl_::shrink(a), t = gcd_mpi_impl_::shrink(b); | |
| 41 |
2/2✓ Branch 2 taken 3 times.
✓ Branch 3 taken 7 times.
|
10 | if (a < b) swap(a, b); |
| 42 | 63 | while (true) { | |
| 43 |
2/2✓ Branch 1 taken 10 times.
✓ Branch 2 taken 63 times.
|
73 | if (b.empty()) break; |
| 44 | if CEXP (FAST) | ||
| 45 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 63 times.
|
63 | if (a.size() <= 2) break; |
| 46 |
2/2✓ Branch 4 taken 53 times.
✓ Branch 5 taken 10 times.
|
63 | if (a = a - b; !a.empty()) |
| 47 | 53 | while (true) { | |
| 48 | 106 | cu32 g = gcd(a[0], mpi::D); | |
| 49 |
2/2✓ Branch 0 taken 53 times.
✓ Branch 1 taken 53 times.
|
106 | if (g == 1) break; |
| 50 |
1/2✓ Branch 0 taken 53 times.
✗ Branch 1 not taken.
|
53 | if (g != mpi::D) a *= mpi::D / g; |
| 51 | 53 | a.erase(begin(a)); | |
| 52 | } | ||
| 53 |
2/2✓ Branch 2 taken 31 times.
✓ Branch 3 taken 32 times.
|
63 | if (a < b) swap(a, b); |
| 54 | } | ||
| 55 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
|
10 | assert(a >= b); |
| 56 | 10 | mpi g; | |
| 57 |
1/2✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
|
10 | if CEXP (FAST) g = b.empty() ? a : gcd(a.to_i64(), b.to_i64()); |
| 58 | else g = a; | ||
| 59 | 10 | cu32 e2 = (u32)min(s.first, t.first), e5 = (u32)min(s.second, t.second); | |
| 60 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 6 times.
|
10 | if (e2) g *= qpow(mpi{2}, e2); |
| 61 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
|
10 | if (e5) g *= qpow(mpi{5}, e5); |
| 62 | 10 | return g; | |
| 63 | 10 | } | |
| 64 | |||
| 65 | } // namespace tifa_libs | ||
| 66 |