library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub Harui-i/library

:heavy_check_mark: test/unittest/unittest-multiple-divisor-moebius-transform.test.cpp

Depends on

Code

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

#include "template/template.hpp"
#include "math/modint.hpp"

#include "convolution/multiple-zeta-moebius-transform.hpp"
#include "convolution/divisor-zeta-moebius-transform.hpp"

using mint = modint998244353;

mint op(mint a, mint b) {
  return a + b;
}

mint invop(mint a, mint b) {
  return a - b;
}

// mapで倍数のゼータ変換とメビウス変換やるやつのテスト
void test() {
  random_device seed_gen;
  mt19937 rng(seed_gen());


  auto map_transform = [&rng](int N) {
    map<ll, mint> multiple_map;
    map<ll, mint> divisor_map;

    vector<mint> multiple_vec(N + 1);
    vector<mint> divisor_vec(N + 1);

    for (ll i = 1; i * i <= N; i++) {
      if (N % i == 0) {
        int r1 = rng();
        int r2 = rng();
        multiple_map[i] = r1;
        multiple_vec[i] = r1;

        multiple_map[N / i] = r2;
        multiple_vec[N / i] = r2;
      }
    }

    divisor_map = multiple_map;
    divisor_vec = multiple_vec;

    // multiple
    {
      map<ll, mint> multiple_map2 = multiple::zeta_transform<ll,mint,op>(multiple_map);
      vector<mint> multiple_vec2 = multiple::zeta_transform_naive<mint,op>(multiple_vec);

      for (pair<ll, mint> plmi : multiple_map2) {
        assert(multiple_vec2[plmi.first] == plmi.second);
      }

      map<ll, mint> multiple_map3 = multiple::moebius_transform<ll,mint,invop>(multiple_map2);
      vector<mint> multiple_vec3 = multiple::moebius_transform_naive<mint,invop>(multiple_vec2);

      assert(multiple_map == multiple_map3 && "multiple transform for map");
      assert(multiple_vec == multiple_vec3 && "multiple transform for vector");

      for (pair<ll, mint> plmi : multiple_map3) {
        assert(multiple_vec3[plmi.first] == plmi.second);
      }
    }

    // divisor
    {
      map<ll, mint> divisor_map2 = divisor::zeta_transform<ll,mint,op>(divisor_map);
      vector<mint> divisor_vec2 = divisor::zeta_transform_naive<mint,op>(divisor_vec);

      for (pair<ll, mint> plmi : divisor_map2) {
        assert(divisor_vec2[plmi.first] == plmi.second && "divisor zeta transform");
      }

      map<ll, mint> divisor_map3 = divisor::moebius_transform<ll,mint,invop>(divisor_map2);
      vector<mint> divisor_vec3 = divisor::moebius_transform_naive<mint,invop>(divisor_vec2);

      assert(divisor_map == divisor_map3 && "multiple transform for map");
      assert(divisor_vec == divisor_vec3 && "multiple transform for vector");

      for (pair<ll, mint> plmi : divisor_map3) {
        assert(divisor_vec3[plmi.first] == plmi.second && "divisor moebius transform");
      }
    }

    };

  for (int i = 0; i < 1000; i++) {
    map_transform(rng() % 10000 + 1);
  }
  return;
}

int main() {
  ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
  test();
  int A, B; cin >> A >> B;
  cout << A + B << endl;
}
#line 1 "test/unittest/unittest-multiple-divisor-moebius-transform.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"

#line 1 "template/template.hpp"
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned int uint;
template<class T> inline bool chmax(T& a, const T& b) {if (a<b) {a=b; return true;} return false;}
template<class T> inline bool chmin(T& a, const T& b) {if (b<a) {a=b; return true;} return false;}
const int INTINF = 1000001000;
const int INTMAX = 2147483647;
const ll LLMAX = 9223372036854775807;
const ll LLINF = 1000000000000000000;
#line 1 "math/modint.hpp"



#line 1 "math/external_gcd.hpp"



#line 5 "math/external_gcd.hpp"

// g,x,y
template<typename T>
constexpr std::tuple<T, T, T> extendedGCD(T a, T b) {
    T x0 = 1, y0 = 0, x1 = 0, y1 = 1;
    while (b != 0) {
        T q = a / b;
        T r = a % b;
        a = b;
        b = r;
        
        T xTemp = x0 - q * x1;
        x0 = x1;
        x1 = xTemp;
        
        T yTemp = y0 - q * y1;
        y0 = y1;
        y1 = yTemp;
    }
    return {a, x0, y0};
}

#line 5 "math/modint.hpp"

template<int MOD>
struct static_modint {
    int value;

    constexpr static_modint() : value(0) {}

    constexpr static_modint(long long v) {
        value = int(((v % MOD) + MOD) % MOD);
    }

    constexpr static_modint& operator+=(const static_modint& other) {
        if ((value += other.value) >= MOD) value -= MOD;
        return *this;
    }

    constexpr static_modint& operator-=(const static_modint& other) {
        if ((value -= other.value) < 0) value += MOD;
        return *this;
    }

    constexpr static_modint& operator*=(const static_modint& other) {
        value = int((long long)value * other.value % MOD);
        return *this;
    }

    constexpr static_modint operator+(const static_modint& other) const {
        return static_modint(*this) += other;
    }

    constexpr static_modint operator-(const static_modint& other) const {
        return static_modint(*this) -= other;
    }

    constexpr static_modint operator*(const static_modint& other) const {
        return static_modint(*this) *= other;
    }

    constexpr static_modint pow(long long exp) const {
        static_modint base = *this, res = 1;
        while (exp > 0) {
            if (exp & 1) res *= base;
            base *= base;
            exp >>= 1;
        }
        return res;
    }

    constexpr static_modint inv() const {
        //return pow(MOD - 2);
        int g,x,y;
        tie(g,x,y) = extendedGCD(value, MOD);
        assert(g==1);
        if (x < 0) {
            x += MOD;
        }
        //cerr << g << " " << x << " " << y << " " << value << endl;
        //assert((((long)x*value)%MOD + MOD)%MOD == 1);
        return x;
    }

    constexpr static_modint& operator/=(const static_modint& other) {
        return *this *= other.inv();
    }

    constexpr static_modint operator/(const static_modint& other) const {
        return static_modint(*this) /= other;
    }

    constexpr bool operator!=(const static_modint& other) const {
        return val() != other.val();
    }

    constexpr bool operator==(const static_modint& other) const {
        return val() == other.val();
    }

    int val() const {
      return this->value;
    }

    friend std::ostream& operator<<(std::ostream& os, const static_modint& mi) {
        return os << mi.value;
    }

    friend std::istream& operator>>(std::istream& is, static_modint& mi) {
        long long x;
        is >> x;
        mi = static_modint(x);
        return is;
    }
};

template <int mod>
using modint = static_modint<mod>;
using modint998244353  = modint<998244353>;
using modint1000000007 = modint<1000000007>;


#line 5 "test/unittest/unittest-multiple-divisor-moebius-transform.test.cpp"

#line 1 "convolution/multiple-zeta-moebius-transform.hpp"



#line 5 "convolution/multiple-zeta-moebius-transform.hpp"

namespace multiple {

  // 倍数についてのゼータ変換。 g_n = \Sigma_{n|m} f_m なる g を求める。
  // n|mというのは、m%n==0という意味。
  // O(N log N) (調和級数)
  // うまくやるとO(Nlog(log(N)))にできることがよく知られているが、難しいしlogは定数なので妥協する。
  template <typename T, T (*op)(T, T)>
  std::vector<T> zeta_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = 1; i <= N; i++) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = op(g[i], f[j]);
      }
    }

    return g;
  }

  // 倍数についてのメビウス変換
  // f_n = \Sigma_{n|m} g_m なる g を求める。
  // O(N log N) (調和級数)
  template <typename T, T(*invop)(T, T)>
  std::vector<T> moebius_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;
    for (int i = N; i >= 1; i--) {
      for (int j = 2 * i; j <= N; j += i) {
        g[i] = invop(g[i], g[j]);
      }
    }
    return g;
  }


  template <typename I, typename T, T(*op)(T, T)>
  std::map<I, T> zeta_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (std::pair<I, T> pit : ret) {
      for (auto p2itr = ret.rbegin(); (*p2itr).first != pit.first; p2itr++) {
        if ((*p2itr).first % pit.first == 0) {
          ret[pit.first] = op(ret[pit.first], (*p2itr).second);
        }
      }
    }

    return ret;
  }


  template <typename I, typename T, T (*invop)(T, T)>
  std::map<I, T> moebius_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (auto p1itr = ret.rbegin(); p1itr != ret.rend(); p1itr++) {
      for (auto p2itr = ret.rbegin(); p2itr != p1itr; p2itr++) {
        if ((*p2itr).first % (*p1itr).first == 0) {
          (*p1itr).second = invop((*p1itr).second, (*p2itr).second);
        }
      }
    }

    return ret;
  }

} // namespace multiple


#line 1 "convolution/divisor-zeta-moebius-transform.hpp"



#line 5 "convolution/divisor-zeta-moebius-transform.hpp"

namespace divisor {
  // 約数についてのゼータ変換。 g_n = \Sigma_{m|n} f_m なる g を求める。
  template <typename T, T(*op)(T, T) >
  std::vector<T> zeta_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;

    for (int i = 1; i <= N; i++) {
      for (int j = 2 * i; j <= N; j += i) {
        g[j] = op(g[j], f[i]);
      }
    }

    return g;
  }

  // 約数についてのメビウス変換。 f_n = \Sigma_{m|n} g_m なる g を求める。
  template <typename T, T(*invop)(T, T)>
  std::vector<T> moebius_transform_naive(const std::vector<T>& f) {
    int N = f.size() - 1;
    std::vector<T> g = f;

    for (int i = 1; i <= N; i++) {
      for (int j = i * 2; j <= N; j += i) {
        g[j] = invop(g[j], g[i]);
      }
    }

    return g;
  }

  template <typename I, typename T, T(*op)(T, T)>
  std::map<I, T> zeta_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;
    for (auto p2itr = mp.rbegin(); p2itr != mp.rend(); p2itr++) {
      for (auto p1itr = next(p2itr); p1itr != mp.rend(); p1itr++) {
        if ((*p2itr).first % (*p1itr).first == 0) {
          ret[(*p2itr).first] = op(ret[(*p2itr).first], ret[(*p1itr).first]);
        }
      }
    }

    return ret;
  }


  template <typename I, typename T, T(*op)(T, T)>
  std::map<I, T> moebius_transform(const std::map<I, T>& mp) {
    std::map<I, T> ret = mp;

    for (auto p1itr = ret.begin(); p1itr != ret.end(); p1itr++) {
      for (auto p2itr = next(p1itr); p2itr != ret.end(); p2itr++) {
        if ((*p2itr).first % (*p1itr).first == 0) {
          ret[(*p2itr).first] = invop(ret[(*p2itr).first], ret[(*p1itr).first]);
        }
      }
    }

    return ret;
  }
} // namespace divisor


#line 8 "test/unittest/unittest-multiple-divisor-moebius-transform.test.cpp"

using mint = modint998244353;

mint op(mint a, mint b) {
  return a + b;
}

mint invop(mint a, mint b) {
  return a - b;
}

// mapで倍数のゼータ変換とメビウス変換やるやつのテスト
void test() {
  random_device seed_gen;
  mt19937 rng(seed_gen());


  auto map_transform = [&rng](int N) {
    map<ll, mint> multiple_map;
    map<ll, mint> divisor_map;

    vector<mint> multiple_vec(N + 1);
    vector<mint> divisor_vec(N + 1);

    for (ll i = 1; i * i <= N; i++) {
      if (N % i == 0) {
        int r1 = rng();
        int r2 = rng();
        multiple_map[i] = r1;
        multiple_vec[i] = r1;

        multiple_map[N / i] = r2;
        multiple_vec[N / i] = r2;
      }
    }

    divisor_map = multiple_map;
    divisor_vec = multiple_vec;

    // multiple
    {
      map<ll, mint> multiple_map2 = multiple::zeta_transform<ll,mint,op>(multiple_map);
      vector<mint> multiple_vec2 = multiple::zeta_transform_naive<mint,op>(multiple_vec);

      for (pair<ll, mint> plmi : multiple_map2) {
        assert(multiple_vec2[plmi.first] == plmi.second);
      }

      map<ll, mint> multiple_map3 = multiple::moebius_transform<ll,mint,invop>(multiple_map2);
      vector<mint> multiple_vec3 = multiple::moebius_transform_naive<mint,invop>(multiple_vec2);

      assert(multiple_map == multiple_map3 && "multiple transform for map");
      assert(multiple_vec == multiple_vec3 && "multiple transform for vector");

      for (pair<ll, mint> plmi : multiple_map3) {
        assert(multiple_vec3[plmi.first] == plmi.second);
      }
    }

    // divisor
    {
      map<ll, mint> divisor_map2 = divisor::zeta_transform<ll,mint,op>(divisor_map);
      vector<mint> divisor_vec2 = divisor::zeta_transform_naive<mint,op>(divisor_vec);

      for (pair<ll, mint> plmi : divisor_map2) {
        assert(divisor_vec2[plmi.first] == plmi.second && "divisor zeta transform");
      }

      map<ll, mint> divisor_map3 = divisor::moebius_transform<ll,mint,invop>(divisor_map2);
      vector<mint> divisor_vec3 = divisor::moebius_transform_naive<mint,invop>(divisor_vec2);

      assert(divisor_map == divisor_map3 && "multiple transform for map");
      assert(divisor_vec == divisor_vec3 && "multiple transform for vector");

      for (pair<ll, mint> plmi : divisor_map3) {
        assert(divisor_vec3[plmi.first] == plmi.second && "divisor moebius transform");
      }
    }

    };

  for (int i = 0; i < 1000; i++) {
    map_transform(rng() % 10000 + 1);
  }
  return;
}

int main() {
  ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
  test();
  int A, B; cin >> A >> B;
  cout << A + B << endl;
}
Back to top page