NumberTheory/Yosupo/Unit_test_math.test.cpp

View this file on GitHub
Last update: 2026-07-26 19:50:18+07:00
Problem: https://judge.yosupo.jp/problem/aplusb

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"

#include "../../template.h"

#include "../Math/Divisors.h"

#include "../Math/Factorial.h"

#include "../Math/Misc.h"


void solve() {
    const int N = 200000;

    // count_divisor vs sieve_count_divisors

    auto dv = sieve_count_divisors(N);
    for (int n = 1; n <= N; n++) assert(count_divisor(n) == dv[n]);

    // sum_of_divisor vs prefix sum of sieve_sum_divisors

    auto sd = sieve_sum_divisors(N);
    ll acc = 0;
    for (int n = 1; n <= N; n++) {
        acc = (acc + sd[n]) % MOD;
        if (n % 4096 == 0 || n == N) assert(sum_of_divisor(n) == acc);
    }

    // factmod vs brute force

    for (ll p : {2LL, 3LL, 97LL, 9973LL}) {
        ll r = 1;
        for (ll n = 1; n <= 5000; n++) {
            ll x = n;
            while (x % p == 0) x /= p;
            r = r * (x % p) % p;
            assert(factmod(n, p) == r);
        }
    }

    // GCD, LCM

    for (ll a = -60; a <= 60; a++) for (ll b = -60; b <= 60; b++) {
        ll g = GCD(a, b), l = LCM(a, b);
        if (a && b) {
            assert(llabs(a) % g == 0 && llabs(b) % g == 0);
            assert(l == llabs(a * b) / g);
        } else assert(l == 0);
    }

    // base conversion round trip

    for (int base = 2; base <= 36; base++)
        for (ll v = 0; v <= 3000; v++)
            assert(convert_decimal(decimal_to_base(v, base), base) == v);

    ll a, b; cin >> a >> b;
    cout << a + b << '\n';
}

#line 1 "NumberTheory/Yosupo/Unit_test_math.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"

#line 2 "template.h"

#include <bits/stdc++.h>

using namespace std;
 
#define ll long long
#define MOD (ll)(1e9+7)
#define all(x) (x).begin(),(x).end()
#define unique(x) x.erase(unique(all(x)), x.end())
#define INF32 ((1ull<<31)-1)
#define INF64 ((1ull<<63)-1)
#define inf (ll)1e18

#define vi vector<int>
#define pii pair<int, int>
#define pll pair<ll, ll>
#define fi first
#define se second

mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
ll get_rand(ll r) { return uniform_int_distribution<ll>(0, r - 1)(rng); }

const int mod = 998244353;

void solve();

int main(){
    ios_base::sync_with_stdio(false);cin.tie(NULL);
    // cin.exceptions(cin.failbit);

    // int t; cin >> t;

    // while(t--)

        solve();
    cerr << "\nTime run: " << 1000 * clock() / CLOCKS_PER_SEC << "ms" << '\n';
    return 0;
}
#line 2 "NumberTheory/Math/Binary_exponentiation.h"

using u128 = __uint128_t;
using i128 = __int128;

ll binMul(ll a, ll b, ll M) { return (i128)a * b % M; }

// long double trick

// require: mantissa 64 bit, x86 gcc/clang

ll binMul2(ll a, ll b, ll M) {
    ll q = (ll)((long double)a * b / M);
    ll r = (ll)((unsigned ll)a * b - (unsigned ll)q * M);
    return r < 0 ? r + M : (r >= M ? r - M : r);
}

ll binMul3(ll a, ll b, ll M) {
    unsigned long long ua = a % M, um = M, res = 0;
    while (b) {
        if (b & 1) { res += ua; if (res >= um) res -= um; }
        ua <<= 1; if (ua >= um) ua -= um;
        b >>= 1;
    }
    return res;
}

ll binPow(ll a, ll b, ll M) {
    a %= M;
    ll res = 1 % M;
    while (b) {
        if (b & 1) res = (i128)res * a % M;
        a = (i128)a * a % M;
        b /= 2;
    }
    return res;
}
#line 3 "NumberTheory/Math/MillerRabin.h"


bool millerTest(ll a, ll n, ll k, ll m) {
    ll mod = binPow(a, m, n);
    if (mod == 1 || mod == n - 1) return true;
    for (int l = 1; l < k; l++) {
        mod = (u128)mod * mod % n;
        if (mod == n - 1) return true;
    }
    return false;
}

// Miller rabin

bool MillerRabin0(ll n) {
    if (n < 4) return n == 2 || n == 3;
    if (n % 2 == 0) return false;
    ll k = 0, m = n - 1;
    while(m % 2 == 0) {
        m /= 2;
        k++;
    }
    for (int i = 0; i < 5; i++) {
        ll a = get_rand(n-3) + 2;
        if (!millerTest(a, n, k, m)) return false;
    }
    return true;
}

// Miller Rabin deterministic version

bool MillerRabin(ll n) {
    if (n <= 1) return false;
    ll k = 0, m = n-1;
    while (m % 2 == 0) {
        m /= 2;
        k++;
    }
    for (int a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        if (n == a) return true;
        if (!millerTest(a, n, k, m)) return false;
    }
    return true;
}
#line 2 "NumberTheory/Math/Sieve.h"

vector<int> sieve(int n) {
    vector<int> nt(n+1, 1);
    nt[0] = nt[1] = 0;
    for (int i = 2; i * i <= n; i++) {
        if (!nt[i]) continue;
        for (int j = i * i; j <= n; j += i)
            nt[j] = 0;
    }
    return nt;
}

vector<int> segmentSieve(ll l, ll r){
    ll sq = sqrtl(r);
    while (sq * sq > r) sq--;
    while (sq+1 <= r/(sq+1)) sq++;
    vector<int> small(sq+1, 1);
    for (ll i = 2; i*i <= sq; i++) {
        if (!small[i]) continue;
        for (ll j = i*i; j <= sq; j += i) {
            small[j] = 0;
        }
    }
    vector<int> prime(r-l+1, 1);
    for(ll p = 2; p <= sq; p++){
        if (!small[p]) continue;
        ll lim = max(p*p, (l+p-1)/p*p);
        for(ll j = lim; j <= r; j += p) {
            prime[j-l] = 0;
        }
    }
    if (l == 0) prime[0] = 0;
    if (l == 0 && r > l) prime[1] = 0;
    if (l == 1) prime[1-l] = 0;
    return prime;
}

vector<int> listPrime(int l, int r) {
    vector<int> prime = segmentSieve(l, r);
    vector<int> listPi;
    for (int i = l; i <= r; i++) {
        if (prime[i-l]) listPi.push_back(i);
    }
    return listPi;
}

vector<int> sieve_count_divisors(int n) {
    vector<int> divisors(n+1, 0);
    for (ll i = 1; i*i <= n; i++) {
        for (ll j = i*i; j <= n; j += i) {
            divisors[j] += 2;
        }
        divisors[i*i]--;
    }
    return divisors;
}

vector<int> sieve_sum_divisors(int n) {
    vector<int> sumDiv(n+1, 0);
    for (int i = 1; i*i <= n; i++) {
        for (int j = i*i; j <= n; j += i) {
            sumDiv[j] += i;
            if (i*i != j) sumDiv[j] += j/i;
        }
    }
    return sumDiv;
}

vector<int> segmentSieveDivisors(int l, int r){
    vector<int> div(r-l+1, 0);
    for(ll i = 1; i*i <= r; i++){
        ll lim = max(i*i, (l+i-1)/i*i);
        for(ll j = lim; j <= r; j += i) {
            div[j-l] += 2;
        }
        if (i * i >= l) div[i*i-l]--;
    }
    if (l == 0) div[0] = 0;
    return div;
}

vector<vector<int>> sieve_prime_divisors(int n) {
    vector<int> prime = sieve(n);
    vector<vector<int>> div(n+1);
    for (int i = 2; i <= n; i++) {
        if (!prime[i]) continue;
        for (int j = i; j <= n; j += i) {
            div[j].push_back(i);
        }
    }
    return div;
}
#line 4 "NumberTheory/Math/Divisors.h"

using u128 = __uint128_t;

// sum of all divisor [1, n]

// https://usaco.guide/problems/cses-1082-sum-of-divisors/solution

ll sum_of_divisor(ll n) {
    ll res = 0, i = 1;
    while (i <= n) {
        ll l = i;
        ll r = n / (n / i);
        ll val = n / l;
        ll csc = (u128)(r - l + 1) * (l + r) / 2 % MOD;
        res = (res + csc * val) % MOD;
        i = r + 1;
    }
    return res;
}

// count number of divisor up to n <= 1e18

// AC: https://codeforces.com/gym/100753 (Probblem F)

int count_divisor(ll n) {
    if (n == 1) return 1;
    static vector<int> prime = listPrime(1, 1e6+5);
    int ans = 1;
    for (int p : prime) {
        if (1ll * p * p * p > n) break;
        int cnt = 0;
        while (n % p == 0) {
            n /= p;
            cnt++;
        }
        ans *= cnt + 1;
    }
    auto isSqrt = [&](ll n)->bool {
        ll c = sqrtl(n);
        return c * c == n;
    };
    if (n == 1) return ans;
    if (MillerRabin(n)) ans *= 2;
    else if (isSqrt(n)) ans *= 3;
    else ans *= 4;
    return ans;
}
#line 2 "NumberTheory/Math/Factorial.h"


// https://wiki.vnoi.info/translate/he/Wilsons-theorem

// Wilson theorem

// n > 1 is prime <=> (n-1)! ≡ -1 (mod n)

//

// Proof:

//    a^(n-2) ≡ a^(-1) (mod n) (Fermat's little theorem)

// => a^(n-2) * a ≡ 1 (mod n)

// Set b = a^(n-2)

// => ab ≡ 1 (mod n)

//

// Have a = b <=> a^2 ≡ 1 (mod n) <=> a = 1 or a = n-1

//

// So if a = 2,3,...,n-2 then a != b

// => we have (n-3)/2 distinct pairs (because with each a, b is unique)

// so we multiple all paris

// => 2.3...(n-2) ≡ 1 (mod n)

// => (n-1)! ≡ n-1 (mod n)


// Proof by contradiction:

// if n is not prime => n have divisors in range [2, n)

//                   => gcd((n-1)!, n) > 1

//                   => gcd(n-1, n) > 1 (contradiction)



// use for small prime p <= 46341

// slow version

int factmod0(int n, int p) {
    vector<int> f(p);
    f[0] = 1;
    for (int i = 1; i < p; i++) {
        f[i] = f[i-1] * i % p;
    }
    int res = 1;
    while (n > 1) {
        if ((n/p) % 2) res = p - res;
        res = res * f[n % p] % p;
        n /= p;
    }
    return res;
}


// Optimize version using static to cache, can work with multiple calling

ll factmod(ll n, ll p) {
    static ll lastp = -1;
    static vector<ll> f;
    if (p != lastp) {
        lastp = p;
        f.assign(p, 1);
        for (ll i = 1; i < p; i++) f[i] = f[i-1] * i % p;
    }
    ll res = 1;
    while (n > 1) {
        if ((n/p) % 2) res = p - res;
        res = res * f[n % p] % p;
        n /= p;
    }
    return res;
}
#line 2 "NumberTheory/Math/Misc.h"


// for both negative + positive val

ll GCD(ll a, ll b) {
    a = llabs(a), b = llabs(b);
    return (!b ? a : GCD(b, a % b));
}
ll LCM(ll a, ll b) { return (!a || !b) ? 0 : llabs(a / GCD(a, b) * b); }

// Only positive

// ll GCD(ll a, ll b) { return (!b ? a : GCD(b, a % b)); }

// ll LCM(ll a, ll b) { return a / GCD(a, b) * b; }


// logb(a)

double log_base(ll a, ll b) { return log(a) / log(b); }

// use for 2 <= BASE <= 36

string decimal_to_base(ll n, int BASE) {
    if (!n) return "0";
    bool neg = (n < 0 ? 1 : 0);
    n = llabs(n);
    string num = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
    string res;
    while(n) {
        res.push_back(num[n % BASE]);
        n /= BASE;
    }
    if (neg) res += '-';
    reverse(all(res));
    return res;
}

// any BASE to decimal

// use ctype lib

ll convert_decimal(string s, int BASE) {
    auto val = [&](char c) { c = toupper(c); return (isdigit(c) ? c-'0' : c-'A'+10); };
    ll n = s.size(), p = 1, res = 0;
    for (int i = n-1; i >= 0; i--) {
        res += val(s[i]) * p;
        p = p * BASE;
    }
    return res;
}
#line 7 "NumberTheory/Yosupo/Unit_test_math.test.cpp"

void solve() {
    const int N = 200000;

    // count_divisor vs sieve_count_divisors

    auto dv = sieve_count_divisors(N);
    for (int n = 1; n <= N; n++) assert(count_divisor(n) == dv[n]);

    // sum_of_divisor vs prefix sum of sieve_sum_divisors

    auto sd = sieve_sum_divisors(N);
    ll acc = 0;
    for (int n = 1; n <= N; n++) {
        acc = (acc + sd[n]) % MOD;
        if (n % 4096 == 0 || n == N) assert(sum_of_divisor(n) == acc);
    }

    // factmod vs brute force

    for (ll p : {2LL, 3LL, 97LL, 9973LL}) {
        ll r = 1;
        for (ll n = 1; n <= 5000; n++) {
            ll x = n;
            while (x % p == 0) x /= p;
            r = r * (x % p) % p;
            assert(factmod(n, p) == r);
        }
    }

    // GCD, LCM

    for (ll a = -60; a <= 60; a++) for (ll b = -60; b <= 60; b++) {
        ll g = GCD(a, b), l = LCM(a, b);
        if (a && b) {
            assert(llabs(a) % g == 0 && llabs(b) % g == 0);
            assert(l == llabs(a * b) / g);
        } else assert(l == 0);
    }

    // base conversion round trip

    for (int base = 2; base <= 36; base++)
        for (ll v = 0; v <= 3000; v++)
            assert(convert_decimal(decimal_to_base(v, base), base) == v);

    ll a, b; cin >> a >> b;
    cout << a + b << '\n';
}