src/nt/mfsum/pi/lib.hpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #pragma once | ||
| 2 | |||
| 3 | #include "../../../math/div64/lib.hpp" | ||
| 4 | #include "../../../math/iroot/sqrt/lib.hpp" | ||
| 5 | #include "../../../util/alias/others/lib.hpp" | ||
| 6 | |||
| 7 | namespace tifa_libs { | ||
| 8 | |||
| 9 | 31 | CEXP u64 pi_min25(u64 n) NE { | |
| 10 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 30 times.
|
31 | retif_((n <= 1) [[unlikely]], 0); |
| 11 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 29 times.
|
30 | retif_((n == 2) [[unlikely]], 1); |
| 12 | 29 | cu32 lim = isqrt(n); | |
| 13 | 29 | u32 vn = (lim + 1) / 2; | |
| 14 | 87 | vecu s(vn), r(vn); | |
| 15 |
2/2✓ Branch 1 taken 2162483 times.
✓ Branch 2 taken 29 times.
|
2162512 | flt_ (u32, cx, 0, vn) s[cx] = cx; |
| 16 |
2/2✓ Branch 1 taken 2162483 times.
✓ Branch 2 taken 29 times.
|
2162512 | flt_ (u32, cx, 0, vn) r[cx] = cx * 2 + 1; |
| 17 | 29 | vecuu l(vn); | |
| 18 |
2/2✓ Branch 1 taken 2162483 times.
✓ Branch 2 taken 29 times.
|
2162512 | flt_ (u32, cx, 0, vn) l[cx] = (n / (cx * 2 + 1) - 1) / 2; |
| 19 | 29 | vecb vis(lim + 1); | |
| 20 | 29 | u32 cnt_p = 0; | |
| 21 |
2/2✓ Branch 0 taken 4277 times.
✓ Branch 1 taken 3 times.
|
4280 | for (u32 p = 3; p <= lim; p += 2) { |
| 22 |
2/2✓ Branch 2 taken 2660 times.
✓ Branch 3 taken 1617 times.
|
4277 | if (vis[p]) continue; |
| 23 | 1617 | cu32 p2 = p * p; | |
| 24 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 1591 times.
|
1617 | if ((u64)p2 * p2 > n) break; |
| 25 | 1591 | vis[p] = true; | |
| 26 |
2/2✓ Branch 2 taken 3281224 times.
✓ Branch 3 taken 1591 times.
|
3282815 | for (u32 cx = p2; cx <= lim; cx += p * 2) vis[cx] = true; |
| 27 | 1591 | u32 ns = 0; | |
| 28 |
2/2✓ Branch 0 taken 48644341 times.
✓ Branch 1 taken 1591 times.
|
48645932 | flt_ (u32, cx, 0, vn) { |
| 29 | 48644341 | cu32 cur = r[cx]; | |
| 30 |
2/2✓ Branch 2 taken 1785337 times.
✓ Branch 3 taken 46859004 times.
|
48644341 | if (vis[cur]) continue; |
| 31 | 46859004 | cu64 d = (u64)cur * p; | |
| 32 |
2/2✓ Branch 1 taken 1439022 times.
✓ Branch 2 taken 45419982 times.
|
46859004 | l[ns] = l[cx] - (d <= lim ? l[s[d / 2] - cnt_p] : s[(div_u64d(n, d) - 1) / 2]) + cnt_p, r[ns++] = cur; |
| 33 | } | ||
| 34 | 1591 | vn = ns; | |
| 35 |
2/2✓ Branch 0 taken 3281224 times.
✓ Branch 1 taken 1591 times.
|
3282815 | for (u32 cx = (lim - 1) / 2, cy = ((lim / p) - 1) | 1; cy >= p; cy -= 2) |
| 36 |
2/2✓ Branch 2 taken 148087094 times.
✓ Branch 3 taken 3281224 times.
|
151368318 | for (cu32 cur = s[cy / 2] - cnt_p, cz = cy * p / 2; cz <= cx; --cx) s[cx] -= cur; |
| 37 | 1591 | ++cnt_p; | |
| 38 | } | ||
| 39 | 29 | l[0] += u64(vn + (cnt_p - 1) * 2) * (vn - 1) / 2; | |
| 40 |
2/2✓ Branch 2 taken 377117 times.
✓ Branch 3 taken 29 times.
|
377146 | flt_ (u32, cx, 1, vn) l[0] -= l[cx]; |
| 41 |
2/2✓ Branch 0 taken 7717 times.
✓ Branch 1 taken 3 times.
|
7720 | flt_ (u32, cx, 1, vn) { |
| 42 | 7717 | cu32 q = r[cx]; | |
| 43 | 7717 | cu64 m = n / q; | |
| 44 | 7717 | cu32 e = s[(m / q - 1) / 2] - cnt_p; | |
| 45 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 7691 times.
|
7717 | if (e < cx + 1) break; |
| 46 | 7691 | u64 t = 0; | |
| 47 |
2/2✓ Branch 3 taken 25914427 times.
✓ Branch 4 taken 7691 times.
|
25922118 | flt_ (u32, cy, cx + 1, e + 1) t += s[(div_u64d(m, r[cy]) - 1) / 2]; |
| 48 | 7691 | l[0] += t - u64(e - cx) * (cnt_p + cx - 1); | |
| 49 | } | ||
| 50 | 29 | return l[0] + 1; | |
| 51 | 29 | } | |
| 52 | |||
| 53 | } // namespace tifa_libs | ||
| 54 |