algo
This documentation is automatically generated by online-judge-tools/verification-helper
tests/Sum_of_Multiplicative_Function.test.cpp
- View this file on GitHub
- Last update: 2025-11-28 15:36:49+07:00
- Problem: https://judge.yosupo.jp/problem/sum_of_multiplicative_function
Depends on
Code
#define PROBLEM "https://judge.yosupo.jp/problem/sum_of_multiplicative_function"
#include "../misc/macros.h"
#include "../math/Min25.h"
#include "../math/ModInt.h"
using Fp = modint<469762049>;
void solve() {
Min25<Fp> solver;
i64 n;
Fp a, b;
cin >> n >> a >> b;
solver.init(n);
cout << solver.solve(a, b, [&](i64 p, int e) {
return a * e + b * p;
}) << '\n';
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(0);
int tc = 1;
cin >> tc;
while (tc--) solve();
}#line 1 "tests/Sum_of_Multiplicative_Function.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/sum_of_multiplicative_function"
#line 1 "misc/macros.h"
// #pragma GCC optimize("Ofast,unroll-loops") // unroll long, simple loops
// #pragma GCC target("avx2,fma") // vectorizing code
// #pragma GCC target("lzcnt,popcnt,abm,bmi,bmi2") // for fast bitset operation
#include <bits/extc++.h>
#include <tr2/dynamic_bitset>
using namespace std;
using namespace __gnu_pbds; // ordered_set, gp_hash_table
// using namespace __gnu_cxx; // rope
// for templates to work
#define all(x) (x).begin(), (x).end()
#define sz(x) (int) (x).size()
#define pb push_back
#define eb emplace_back
using i32 = int32_t;
using u32 = uint32_t;
using i64 = int64_t;
using u64 = uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;
using ld = long double;
using pii = pair<i32, i32>;
using vi = vector<i32>;
// fast map
const int RANDOM = chrono::high_resolution_clock::now().time_since_epoch().count();
struct chash { // customize hash function for gp_hash_table
int operator()(int x) const { return x ^ RANDOM; }
};
gp_hash_table<int, int, chash> table;
/* ordered set
find_by_order(k): returns an iterator to the k-th element (0-based)
order_of_key(k): returns the number of elements in the set that are strictly less than k
*/
template <class T>
using ordered_set = tree<T, null_type, less<T>, rb_tree_tag, tree_order_statistics_node_update>;
/* rope
rope <int> cur = v.substr(l, r - l + 1);
v.erase(l, r - l + 1);
v.insert(v.mutable_begin(), cur);
*/
#line 1 "math/Min25.h"
template <class T>
struct Min25 {
i64 n;
int sq;
vector<int> primes, id1, id2;
vector<i64> vals;
vector<T> g0, g1; // g0: sum p^0, g1: sum p^1
int id(i64 x) { return x <= sq ? id1[x] : id2[n / x]; }
void init(i64 N) {
n = N, sq = sqrt(n);
primes.clear();
vector<bool> is_p(sq + 1, true);
for (int i = 2; i <= sq; ++i) {
if (is_p[i]) {
primes.pb(i);
for (int j = i * 2; j <= sq; j += i) is_p[j] = false;
}
}
vals.clear(), id1.assign(sq + 1, 0), id2.assign(sq + 1, 0);
for (i64 l = 1, r; l <= n; l = r + 1) {
i64 v = n / l;
r = n / v;
vals.pb(v);
if (v <= sq) id1[v] = sz(vals) - 1;
else id2[n / v] = sz(vals) - 1;
}
g0.resize(sz(vals)), g1.resize(sz(vals));
T inv2 = T(1) / T(2);
for (int i = 0; i < sz(vals); ++i) {
T v = T(vals[i]);
g0[i] = v - 1;
g1[i] = v * (v + 1) * inv2 - 1;
}
for (int p : primes) {
T sp0 = g0[id(p - 1)], sp1 = g1[id(p - 1)];
i64 p2 = (i64) p * p;
T tp = T(p);
for (int i = 0; i < sz(vals); ++i) {
if (vals[i] < p2) break;
int k = id(vals[i] / p);
g0[i] -= g0[k] - sp0;
g1[i] -= tp * (g1[k] - sp1);
}
}
}
// A, B: f(p) = A*1 + B*p
// func: (p, e) -> f(p^e) trả về T
template <class Func>
T solve(T A, T B, Func f_pe) {
vector<T> s_fp(sz(primes) + 1);
for (int i = 0; i < sz(primes); ++i)
s_fp[i + 1] = s_fp[i] + A + B * T(primes[i]);
auto S = [&](auto&& self, i64 x, int j) -> T {
if (x <= 1 || (j < sz(primes) && primes[j] > x)) return 0;
int k = id(x);
T ans = A * g0[k] + B * g1[k];
ans -= s_fp[j];
for (int i = j; i < sz(primes); ++i) {
i64 p = primes[i];
if (p * p > x) break;
i64 pe = p;
for (int e = 1; pe * p <= x; ++e) {
ans += f_pe(p, e) * self(self, x / pe, i + 1);
ans += f_pe(p, e + 1);
pe *= p;
}
}
return ans;
};
return S(S, n, 0) + 1;
}
};
#line 2 "math/ModInt.h"
template <int mod>
struct modint {
using M = modint;
static_assert(mod > 0 && mod <= 2147483647);
static constexpr int modulo = mod;
static constexpr u32 r1 = []() {
u32 r1 = mod;
for (int i = 0; i < 5; ++i) r1 *= 2 - mod * r1;
return -r1;
}();
static constexpr u32 r2 = -u64(mod) % mod;
static u32 reduce(u64 x) {
u32 y = u32(x) * r1, r = (x + u64(y) * mod) >> 32;
return r >= mod ? r - mod : r;
}
u32 x;
modint() : x(0) {}
modint(i64 x) : x(reduce(u64(x % mod + mod) * r2)) {}
M& operator+=(const M& a) {
if ((x += a.x) >= mod) x -= mod;
return *this;
}
M& operator-=(const M& a) {
if ((x += mod - a.x) >= mod) x -= mod;
return *this;
}
M& operator*=(const M& a) {
x = reduce(u64(x) * a.x);
return *this;
}
M& operator/=(const M& a) { return *this *= a.inv(); }
M operator-() const { return M(0) - *this; }
M operator+(const M& a) const { return M(*this) += a; }
M operator-(const M& a) const { return M(*this) -= a; }
M operator*(const M& a) const { return M(*this) *= a; }
M operator/(const M& a) const { return M(*this) /= a; }
bool operator==(const M& a) const { return x == a.x; }
bool operator!=(const M& a) const { return x != a.x; }
M pow(u64 k) const {
M res(1), b = *this;
while (k) {
if (k & 1) res *= b;
b *= b, k >>= 1;
}
return res;
}
M inv() const { return pow(mod - 2); }
friend ostream& operator<<(ostream& os, const M& a) {
return os << reduce(a.x);
}
friend istream& operator>>(istream& is, M& a) {
i64 v;
is >> v;
a = M(v);
return is;
}
};
u64 modmul(u64 x, u64 y, u64 m) { return u128(x) * y % m; }
u64 modpow(u64 x, u64 k, u64 m) {
u64 res = 1;
while (k) {
if (k & 1) res = modmul(res, x, m);
x = modmul(x, x, m);
k >>= 1;
}
return res;
}
#line 6 "tests/Sum_of_Multiplicative_Function.test.cpp"
using Fp = modint<469762049>;
void solve() {
Min25<Fp> solver;
i64 n;
Fp a, b;
cin >> n >> a >> b;
solver.init(n);
cout << solver.solve(a, b, [&](i64 p, int e) {
return a * e + b * p;
}) << '\n';
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(0);
int tc = 1;
cin >> tc;
while (tc--) solve();
}