#1066948 (C++17) No.3106 Simple Math Problem 3

提出ソース
結果

問題	No.3106 Simple Math Problem 3
ユーザー	ecottea
提出日時	2025-04-13 05:04:02
言語	C++17 (gcc 13.3.0 + boost 1.87.0)
結果	AC
実行時間	397 ms / 2,000 ms
コード長	15,991 bytes
コンパイル時間	6,262 ms
コンパイル使用メモリ	300,284 KB
実行使用メモリ	61,036 KB
最終ジャッジ日時	2025-04-13 05:04:23
合計ジャッジ時間	17,695 ms
ジャッジサーバーID （参考情報）	judge5 / judge1
このコードへのチャレンジ
（要ログイン）
ファイルパターン	結果
sample	AC * 4
other	AC * 42
権限があれば一括ダウンロードができます
ソースコード

raw source code
#pragma GCC optimize("O3") // たまにバグる
#pragma GCC optimize("unroll-loops")


#ifndef HIDDEN_IN_VS // 折りたたみ用

// 警告の抑制
#define _CRT_SECURE_NO_WARNINGS

// ライブラリの読み込み
#include <bits/stdc++.h>
using namespace std;

// 型名の短縮
using ll = long long; using ull = unsigned long long; // -2^63 ～ 2^63 = 9e18（int は -2^31 ～ 2^31 = 2e9）
using pii = pair<int, int>;	using pll = pair<ll, ll>;	using pil = pair<int, ll>;	using pli = pair<ll, int>;
using vi = vector<int>;		using vvi = vector<vi>;		using vvvi = vector<vvi>;	using vvvvi = vector<vvvi>;
using vl = vector<ll>;		using vvl = vector<vl>;		using vvvl = vector<vvl>;	using vvvvl = vector<vvvl>;
using vb = vector<bool>;	using vvb = vector<vb>;		using vvvb = vector<vvb>;
using vc = vector<char>;	using vvc = vector<vc>;		using vvvc = vector<vvc>;
using vd = vector<double>;	using vvd = vector<vd>;		using vvvd = vector<vvd>;
template <class T> using priority_queue_rev = priority_queue<T, vector<T>, greater<T>>;
using Graph = vvi;

// 定数の定義
const double PI = acos(-1);
int DX[4] = { 1, 0, -1, 0 }; // 4 近傍（下，右，上，左）
int DY[4] = { 0, 1, 0, -1 };
int INF = 1001001001; ll INFL = 4004004003094073385LL; // (int)INFL = INF, (int)(-INFL) = -INF;

// 入出力高速化
struct fast_io { fast_io() { cin.tie(nullptr); ios::sync_with_stdio(false); cout << fixed << setprecision(18); } } fastIOtmp;

// 汎用マクロの定義
#define all(a) (a).begin(), (a).end()
#define sz(x) ((int)(x).size())
#define lbpos(a, x) (int)distance((a).begin(), std::lower_bound(all(a), (x)))
#define ubpos(a, x) (int)distance((a).begin(), std::upper_bound(all(a), (x)))
#define Yes(b) {cout << ((b) ? "Yes\n" : "No\n");}
#define rep(i, n) for(int i = 0, i##_len = int(n); i < i##_len; ++i) // 0 から n-1 まで昇順
#define repi(i, s, t) for(int i = int(s), i##_end = int(t); i <= i##_end; ++i) // s から t まで昇順
#define repir(i, s, t) for(int i = int(s), i##_end = int(t); i >= i##_end; --i) // s から t まで降順
#define repe(v, a) for(const auto& v : (a)) // a の全要素（変更不可能）
#define repea(v, a) for(auto& v : (a)) // a の全要素（変更可能）
#define repb(set, d) for(int set = 0, set##_ub = 1 << int(d); set < set##_ub; ++set) // d ビット全探索（昇順）
#define repis(i, set) for(int i = lsb(set), bset##i = set; i < 32; bset##i -= 1 << i, i = lsb(bset##i)) // set の全要素（昇順）
#define repp(a) sort(all(a)); for(bool a##_perm = true; a##_perm; a##_perm = next_permutation(all(a))) // a の順列全て（昇順）
#define uniq(a) {sort(all(a)); (a).erase(unique(all(a)), (a).end());} // 重複除去
#define EXIT(a) {cout << (a) << endl; exit(0);} // 強制終了
#define inQ(x, y, u, l, d, r) ((u) <= (x) && (l) <= (y) && (x) < (d) && (y) < (r)) // 半開矩形内判定

// 汎用関数の定義
template <class T> inline ll powi(T n, int k) { ll v = 1; rep(i, k) v *= n; return v; }
template <class T> inline bool chmax(T& M, const T& x) { if (M < x) { M = x; return true; } return false; } // 最大値を更新（更新されたら true を返す）
template <class T> inline bool chmin(T& m, const T& x) { if (m > x) { m = x; return true; } return false; } // 最小値を更新（更新されたら true を返す）
template <class T> inline T getb(T set, int i) { return (set >> i) & T(1); }
template <class T> inline T smod(T n, T m) { n %= m; if (n < 0) n += m; return n; } // 非負mod

// 演算子オーバーロード
template <class T, class U> inline istream& operator>>(istream& is, pair<T, U>& p) { is >> p.first >> p.second; return is; }
template <class T> inline istream& operator>>(istream& is, vector<T>& v) { repea(x, v) is >> x; return is; }
template <class T> inline vector<T>& operator--(vector<T>& v) { repea(x, v) --x; return v; }
template <class T> inline vector<T>& operator++(vector<T>& v) { repea(x, v) ++x; return v; }

#endif // 折りたたみ用


#if __has_include(<atcoder/all>)
#include <atcoder/all>
using namespace atcoder;

#ifdef _MSC_VER
#include "localACL.hpp"
#endif

using mint = modint998244353;
//using mint = static_modint<(int)1e9 + 7>;
//using mint = modint; // mint::set_mod(m);

namespace atcoder {
	inline istream& operator>>(istream& is, mint& x) { ll x_; is >> x_; x = x_; return is; }
	inline ostream& operator<<(ostream& os, const mint& x) { os << x.val(); return os; }
}
using vm = vector<mint>; using vvm = vector<vm>; using vvvm = vector<vvm>; using vvvvm = vector<vvvm>; using pim = pair<int, mint>;
#endif


#ifdef _MSC_VER // 手元環境（Visual Studio）
#include "local.hpp"
#else // 提出用（gcc）
inline int popcount(int n) { return __builtin_popcount(n); }
inline int popcount(ll n) { return __builtin_popcountll(n); }
inline int lsb(int n) { return n != 0 ? __builtin_ctz(n) : 32; }
inline int lsb(ll n) { return n != 0 ? __builtin_ctzll(n) : 64; }
inline int msb(int n) { return n != 0 ? (31 - __builtin_clz(n)) : -1; }
inline int msb(ll n) { return n != 0 ? (63 - __builtin_clzll(n)) : -1; }
#define dump(...)
#define dumpel(v)
#define dump_math(v)
#define input_from_file(f)
#define output_to_file(f)
#define Assert(b) { if (!(b)) { vc MLE(1<<30); EXIT(MLE.back()); } } // RE の代わりに MLE を出す
#endif


ll naive(ll N) {
	ll res = 0;

	repi(A, 1, N) {
		ll sum = 0;

		repi(B, 1, A - 1) {
			if (B * (A % B) < A) {
				sum += A % B;
			}
		}
		//dump(A, sum);

		res += sum;
	}

	return res;
}


ll TLE(ll N) {
	ll res = 0;

	repi(A, 1, N) {
		ll sum = 0;

		repi(B, 1, A - 1) {
			if (A / B >= (A + B) / (B + 1)) {
				sum += A % B;
			}
		}
		//dump(A, sum);

		res += sum;
	}

	return res;
}


ll TLE2(ll N) {
	ll res = 0;

	repi(A, 1, N) {
		ll sum = 0;

		repi(B, 1, A - 1) {
			if (A % B <= A / (B + 1)) {
				sum += A % B;
			}
		}
		//dump(A, sum);

		res += sum;
	}

	return res;
}


ll TLE3(ll N) {
	ll res = 0;

	repi(m, 1, N / 2) {
		ll sum = 0;

		repi(B, m + 1, N / m - 1) {
			ll add = (N - m) / B - m + 1;

			sum += add;
		}
		//dump(m, sum);

		res += m * sum;
	}

	return res;
}


mint TLE4(ll N) {
	mint res = 0;

	int sqrtN = (int)sqrt(N);

	repi(m, 1, sqrtN) {
		mint sum = 0;

		repi(B, m + 1, N / m - 1) {
			sum += (N - m) / B;
		}
		//dump(m, sum);

		sum += mint(-m + 1) * max(N / m - m - 1, 0LL);

		res += m * sum;
	}

	return res;
}


//【商列挙】O(√N)
/*
* 区間 [1..N] を N/i = q（切り捨て）となる半開区間 i∈(il..ir] に分割し，
* i について昇順にそれぞれに対して f(il, ir, q) を呼び出す．
* なお各範囲においては N mod i は公差 -q の等差数列を成す．
*/
template <class T, class FUNC>
void quotient_range(T N, const FUNC& f) {
	// 参考 : https://ei1333.github.io/luzhiled/snippets/math/quotient-range.html
	// verify : https://judge.yosupo.jp/problem/enumerate_quotients

	//【方法】
	// N/i の商が q となるような i の範囲を考える．条件を i について整理すると
	//		q = floor(N/i)
	//		⇔ q ≦ N/i < q+1
	//		⇔ i q ≦ N < i(q+1)
	//		⇔ N/(q+1) < i ≦ N/q  (⇔ floor(N/(q+1)) < i ≦ floor(N/q))
	// となる．
	//
	// この幅が 1 以下であれば，q に対応する i は高々 1 個である．その条件は
	//		N/q - N/(q+1) ≦ 1
	//		⇔ (q+1)N - q N ≦ q(q+1)
	//		⇔ N ≦ q(q+1)
	// である．条件をやや弱めて
	//		N ≦ q^2 ⇔ √N ≦ q
	// としてもオーダーに影響はない．

	//（例）
	// 例えば N = 15 のときは (0..15] を以下のように分割できる：
	//		i の範囲		q=N/i	N mod i
	//		(0..1]		15		[0]
	//		(1..2]		7		[1]
	//		(2..3]		5		[0]
	//		(3..5]		3		[3, 0]
	//		(5..7]		2		[3, 1]
	//		(7..15]		1		[7, 6, 5, 4, 3, 2, 1, 0]

	T sqrt_n = (T)(sqrt(N) - 1e-9);

	// q に対応する i が高々 1 個の部分は i ごとに愚直に考える．
	T i_max = N / (sqrt_n + 1);
	for (T i = 1; i <= i_max; ++i) f(i - 1, i, N / i);

	// そうでない部分は q ごとにまとめて考える．
	T il, ir = i_max;
	for (T q = sqrt_n; q >= 1; --q) {
		il = ir;
		ir = N / q;
		f(il, ir, q);
	}

	/* f の定義の雛形
	using T = ll;
	auto f = [&](T il, T ir, T q) {

	};
	quotient_range(N, f);
	*/
}


mint TLE5(ll N) {
	mint res = 0;

	int sqrtN = (int)sqrt(N);

	int cnt = 0;

	repi(m, 1, sqrtN) {
		mint sum = 0;

		//repi(B, m + 1, N / m - 1) {
		//	sum += (N - m) / B;
		//}
				
		using T = ll;
		auto f = [&](T il, T ir, T q) {
			cerr << ir << " ";
			chmax<ll>(il, m);
			chmin<ll>(ir, N / m - 1);
			if (il < ir) {
				sum += mint(q) * (ir - il);
				cnt++;
			}
		};
		quotient_range(N - m, f);
		cerr << endl;

		dump(m, sum);
		
		sum += mint(-m + 1) * max(N / m - m - 1, 0LL);

		res += m * sum;
	}

	//dump("cnt:", cnt); // O(N)

	return res;
}
/*
N=100:
m=1: 2 3 4 5 6 7 8 9 11 12 14 16 19 24 33 49 99		99
m=2:   3 4 5 6 7 8 9 10 12 14 16 19 24 32 49		98
m=3:     4 5 6 7 8 9 10 12 13 16 19 24 32			97
m=4:       5 6 7 8 9 10 12 13 16 19 24				96
m=5:         6 7 8 9 10 11 13 15 19					95
m=6:           7 8 9 10 11 13 15					94
m=7:             8 9 10 11 13						93
m=8:               9 10 11							92
m=9:                 10								91
約数のとこだけしか変わっていない．差分更新できそう？

N=500
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 26 27 29 31 33 35 38 41 45 49 55 62 71 83 99 124 166 249 499
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 26 27 29 31 33 35 38 41 45 49 55 62 71 83 99 124 166 249 498
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 26 27 29 31 33 35 38 41 45 49 55 62 71 82 99 124 165 248 497
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 26 27 29 31 33 35 38 41 45 49 55 62 70 82 99 124 165 248 496
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 26 27 29 30 33 35 38 41 45 49 55 61 70 82 99 123 165 247 495
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 26 27 29 30 32 35 38 41 44 49 54 61 70 82 98 123 164 247 494
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 29 30 32 35 37 41 44 49 54 61 70 82 98 123 164 246 493
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 30 32 35 37 41 44 49 54 61 70 82 98 123 164 246 492
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 30 32 35 37 40 44 49 54 61 70 81 98 122 163 245 491
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 30 32 35 37 40 44 49 54 61 70 81 98 122 163 245 490
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 30 32 34 37 40 44 48 54 61 69 81 97 122 163 244 489
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 30 32 34 37 40 44 48 54 61 69 81 97 122 162 244 488
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 30 32 34 37 40 44 48 54 60 69 81 97 121 162 243 487
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 30 32 34 37 40 44 48 54 60 69 81 97 121 162 243 486
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 28 30 32 34 37 40 44 48 53 60 69 80 97 121 161 242 485
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 28 30 32 34 37 40 44 48 53 60 69 80 96 121 161 242 484
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 23 24 25 26 28 30 32 34 37 40 43 48 53 60 69 80 96 120 161 241 483
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 24 25 26 28 30 32 34 37 40 43 48 53 60 68 80 96 120 160 241 482
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 24 25 26 28 30 32 34 37 40 43 48 53 60 68 80 96 120 160 240 481
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 24 25 26 28 30 32 34 36 40 43 48 53 60 68 80 96 120 160 240 480
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 25 26 28 29 31 34 36 39 43 47 53 59 68 79 95 119 159 239 479
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 25 26 28 29 31 34 36 39 43 47 53 59 68 79 95 119 159 239 478
総更新回数は調和級数で抑えられそう．
*/


mint TLE6(ll N) {
	int sqrtN = (int)sqrt(N);

	mint sum = 0;

	vl x{ 0 };
	vm w;

	using T = ll;
	auto f = [&](T il, T ir, T q) {
		x.push_back(ir);
		w.push_back(q);

		chmax<ll>(il, 1);
		chmin<ll>(ir, N / 1 - 1);
		if (il < ir) {
			sum += mint(q) * (ir - il);
		}
	};
	quotient_range(N - 1, f);
	w.push_back(0);
	//dump("x:", x);
	//dump("w:", w);
	//dump("sum:", sum);

	mint res = sum;

	vi r;
	repi(i, 1, sqrtN) r.push_back(N % i);
	//dump("r:", r);

	repi(m, 2, sqrtN) {
		if (N / m - m - 1 <= 0) break;
		//dump("-------------- m:", m, "-----------------");

		//repi(B, m + 1, N / m - 1) {
		//	sum += (N - m) / B;
		//}
		//dump("B_min:", m + 1, "Bmax:", N / m - 1);

		repi(i, 0, sqrtN - 1) {
			r[i] = smod(r[i] - 1, i + 1);
			if (r[i] == 0) {
				x[sz(x) - 1 - i]--;

				int pos = sz(x) - 1 - i - 1;
				if (m - 1 <= pos && pos <= sz(x) - m) {
					sum -= w[sz(x) - 1 - i - 1];
				}

				pos = sz(x) - 1 - i;
				if (m - 1 <= pos && pos <= sz(x) - m) {
					sum += w[pos];
				}

				w[i]--;

				if (m - 1 <= i && i <= sz(x) - m) {
					sum -= x[i + 1] - x[i];
				}
			}
		}
		//dump("sum:", sum);

		sum -= w[m - 1] * (x[m] - x[m - 1]);
		//dump("sum:", sum);

		int pos = sz(w) - 1 - (m - 1);
		sum -= w[pos] * (x[pos + 1] - x[pos]);
		//dump("sum:", sum);

		res += m * sum;

		//dump("x:", x);
		//dump("w:", w);
		//dump("sum:", sum);
	}

	repi(m, 1, sqrtN) {
		if (N / m - m - 1 <= 0) break;
		res += mint(m) * (-m + 1) * (N / m - m - 1);
	}

	return res;
}


mint solve(ll N) {
	int sqrtN = (int)sqrt(N);

	mint sum = 0;

	vl x{ 0 };
	vm w;
	x.reserve(2 * sqrtN + 10);
	w.reserve(2 * sqrtN + 10);

	using T = ll;
	auto f = [&](T il, T ir, T q) {
		x.push_back(ir);
		w.push_back(q);

		if (il == 0) return;
		sum += mint(q) * (ir - il);
	};
	quotient_range(N - 1, f);
	w.push_back(0);
//	dump("x:", x);
//	dump("w:", w);
//	dump("sum:", sum);

	int L = sz(x);

	mint res = sum;

	vvi is(sqrtN + 1);
	repi(i, 1, sqrtN) {
		int m0 = N % i;
		int m_max = min(i + 10, sqrtN - 1);
		for (int m = m0 + 1; m <= m_max; m += i) {
			is[m].push_back(i - 1);
		}
	}
//	dumpel(is);

	repi(m, 2, sqrtN) {
		if (N / m - m - 1 <= 0) break;
//		dump("-------------- m:", m, "-----------------");

		//repi(B, m + 1, N / m - 1) {
		//	sum += (N - m) / B;
		//}
//		dump("B_min:", m + 1, "Bmax:", N / m - 1);

		repe(i, is[m]) {
			x[L - 1 - i]--;

			int pos = L - 1 - i - 1;
			if (m - 1 <= pos && pos <= L - m) {
				sum -= w[pos];
			}

			pos = L - 1 - i;
			if (m - 1 <= pos && pos <= L - m) {
				sum += w[pos];
			}

			w[i]--;

			if (m - 1 <= i && i <= L - m) {
				sum -= x[i + 1] - x[i];
			}
		}
//		dump("sum:", sum);

		sum -= w[m - 1] * (x[m] - x[m - 1]);
//		dump("sum:", sum);

		int pos = L - 1 - (m - 1);
		sum -= w[pos] * (x[pos + 1] - x[pos]);
//		dump("sum:", sum);

		res += m * sum;

//		dump("x:", x);
//		dump("w:", w);
//		dump("sum:", sum);
	}

	repi(m, 1, sqrtN) {
		if (N / m - m - 1 <= 0) break;
		res += mint(m) * (-m + 1) * (N / m - m - 1);
	}

	return res;
}


void bug_find() {
#ifdef _MSC_VER
	// 合わない入力例を見つける．

	mute_dump = true;

	mt19937_64 mt((int)time(NULL));
	uniform_int_distribution<ll> rnd(0LL, 1LL << 60);

	rep(hoge, 1000) {
		ll n = rnd(mt) % 1000 + 1;

		auto res_naive = naive(n);
		auto res_solve = solve(n);

		if (res_naive != res_solve) {
			cout << "----------error!----------" << endl;
			cout << "input:" << endl;
			cout << n << endl;
			cout << "results:" << endl;
			cout << res_naive << endl;
			cout << res_solve << endl;
			cout << "--------------------------" << endl;
		}
	}

	mute_dump = false;
	exit(0);
#endif
}


int main() {
	input_from_file("input.txt");
//	output_to_file("output.txt");

//	bug_find();

	ll n;
	cin >> n;

//	dump(naive(n)); dump("=====");
//	dump(TLE(n)); dump("=====");
//	dump(TLE2(n)); dump("=====");
//	dump(TLE3(n)); dump("=====");
//	dump(TLE4(n)); dump("=====");
//	dump(TLE5(n)); dump("=====");
//	dump(TLE6(n)); dump("=====");

	EXIT(solve(n));
}
yukicoder

結果

ソースコード
