fork download
copy 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define ll long long int
  4. #define endl "\n"
  5.  
  6. // Use mult64() if the number fits in 64 bits, otherwise change all to mult128() (slower)
  7. inline ll mult64(const ll &a, const ll &b, const ll &mod)
  8. {
  9.     return [&]() -> ll
  10.     {
  11.         return mod > INT32_MAX ? (ll)((__int128_t)(a)*b % mod) : (long long)(a)*b % mod;
  12.     }();
  13. }
  14. namespace bigNumber
  15. {
  16.     using u128 = __uint128_t;
  17.  
  18.     // ---------- 64-bit limb helpers (fast) ----------
  19.     static inline void mult64(uint64_t a, uint64_t b, uint64_t &lo, uint64_t &hi)
  20.     {
  21.         __uint128_t p = (__uint128_t)a * b;
  22.         lo = (uint64_t)p;
  23.         hi = (uint64_t)(p >> 64);
  24.     }
  25.     static inline uint64_t add64(uint64_t a, uint64_t b, uint64_t &carry)
  26.     {
  27.         __uint128_t s = (__uint128_t)a + b + carry;
  28.         carry = (uint64_t)(s >> 64);
  29.         return (uint64_t)s;
  30.     }
  31.     // Fallback generic (slow): kept for setup (e.g., computing R^2 mod N once)
  32.     u128 mult128(u128 a, u128 b, u128 mod)
  33.     {
  34.         u128 result = 0;
  35.         for (a %= mod; b > 0; a <<= 1, b >>= 1)
  36.         {
  37.             a >= mod ? a -= mod : 0;
  38.             if (b & 1)
  39.                 result += a, result >= mod ? result -= mod : 0;
  40.         }
  41.         return result;
  42.     }
  43.     // ---------- Montgomery (CIOS, 2x64-bit limbs) ----------
  44.     // Compute n0' = -N^{-1} mod 2^64 (low limb only)
  45.     static inline uint64_t inv64_2k(uint64_t n0)
  46.     {
  47.         uint64_t x = 1; // initial approx
  48.         for (int i = 6; i > 0; i--)
  49.             x = (__uint128_t)x * (2 - (__uint128_t)n0 * x);
  50.         return x; // x ≡ n0^{-1}
  51.     }
  52.  
  53.     inline pair<u128, u128> montModInv(u128 N)
  54.     {
  55.         return {0, (u128)(0 - inv64_2k(N))};
  56.     }
  57.  
  58.     // Fast Montgomery multiplication: returns a*b*R^{-1} mod N, where R=2^128
  59.     inline u128 montMult(u128 a, u128 b, u128 N, u128 N_n0prime)
  60.     {
  61.         uint64_t n0 = (uint64_t)N, n1 = (uint64_t)(N >> 64);
  62.         uint64_t a0 = (uint64_t)a, a1 = (uint64_t)(a >> 64);
  63.         uint64_t b0 = (uint64_t)b, b1 = (uint64_t)(b >> 64);
  64.         uint64_t n0p = (uint64_t)N_n0prime;
  65.         uint64_t t0 = 0, t1 = 0, t2 = 0;
  66.  
  67.         auto round_step = [&](uint64_t ai)
  68.         {
  69.             uint64_t lo, hi, carry;
  70.             // t += ai * b
  71.             mult64(ai, b0, lo, hi);
  72.             carry = 0;
  73.             t0 = add64(t0, lo, carry);
  74.             t1 = add64(t1, hi, carry);
  75.             t2 = add64(t2, 0, carry);
  76.             mult64(ai, b1, lo, hi);
  77.             carry = 0;
  78.             t1 = add64(t1, lo, carry);
  79.             t2 = add64(t2, hi, carry);
  80.             // m = t0 * n0' (mod 2^64)
  81.             uint64_t m = (uint64_t)((__uint128_t)t0 * n0p);
  82.             // t += m * N
  83.             mult64(m, n0, lo, hi);
  84.             carry = 0;
  85.             t0 = add64(t0, lo, carry);
  86.             t1 = add64(t1, hi, carry);
  87.             t2 = add64(t2, 0, carry);
  88.             mult64(m, n1, lo, hi);
  89.             carry = 0;
  90.             t1 = add64(t1, lo, carry);
  91.             t2 = add64(t2, hi, carry);
  92.             // shift by one limb
  93.             t0 = t1;
  94.             t1 = t2;
  95.             t2 = 0;
  96.         };
  97.  
  98.         round_step(a0);
  99.         round_step(a1);
  100.  
  101.         __uint128_t res = (((__uint128_t)t1 << 64) | t0);
  102.         if (res >= N)
  103.             res -= N;
  104.         return (u128)res;
  105.     }
  106. }
  107. using namespace bigNumber;
  108.  
  109. template <typename T>
  110. inline T F(T x, T c, T mod, T inv) // Pollard-rho function
  111. {
  112.     x = montMult(x, x, mod, inv);
  113.     x = x >= mod - c ? x - mod + c : x + c;
  114.     return x;
  115. }
  116.  
  117. template <typename T>
  118. inline T __abs(T N)
  119. {
  120.     if (N < 0)
  121.         return -N;
  122.  
  123.     return N;
  124. }
  125.  
  126. template <typename T>
  127. T Pollard_Brent(T N)
  128. {
  129.     if (!(N & 1))
  130.         return 2;
  131.  
  132.     // Random Number Linear Congruential Generator MMIX from D.E. Knuth
  133.     static u128 rng = 0xdeafbeefff;
  134.     uint64_t a = rng * 6364136223846793005ull + 1442695040888963407ull;
  135.     uint64_t b = a * 6364136223846793005ull + 1442695040888963407ull;
  136.     rng = (a + b) ^ (a * b);
  137.  
  138.     T X0 = 1 + a % (N - 1);
  139.     T C = 1 + b % (N - 1);
  140.     T X = X0; // X1
  141.     T gcdVal = 1;
  142.     T q = 1;
  143.     T Xs, Xt;
  144.     T m = 128;
  145.     u128 inv = montModInv(N).second;
  146.     T L = 1;
  147.     while (gcdVal == 1)
  148.     {
  149.         Xt = X;
  150.         for (size_t i = 1; i < L; i++)
  151.             X = F(X, C, N, inv);
  152.  
  153.         int k = 0;
  154.         while (k < L && gcdVal == 1)
  155.         {
  156.             Xs = X;
  157.             for (size_t i = 0; i < m && i < L - k; i++)
  158.             {
  159.                 X = F(X, C, N, inv);
  160.                 q = montMult(q, Xt > X ? Xt - X : X - Xt, N, inv);
  161.             }
  162.             gcdVal = __gcd(q, N);
  163.             k += m;
  164.         }
  165.         L += L;
  166.     }
  167.     if (gcdVal == N) // Failure
  168.     {
  169.         do
  170.         {
  171.             Xs = F(Xs, C, N, inv);
  172.             gcdVal = __gcd(Xs > Xt ? Xs - Xt : Xt - Xs, N);
  173.         } while (gcdVal == 1);
  174.     }
  175.     return gcdVal;
  176. }
  177.  
  178. template <typename T>
  179. T modPow(T N, T power, T mod)
  180. {
  181.     if (N % mod == 0 || N == 0)
  182.         return 0;
  183.     if (N == 1 || power == 0)
  184.         return 1;
  185.     T res{1};
  186.     while (power)
  187.     {
  188.         if (power & 1)
  189.             res = mult64(res, N, mod);
  190.         N = mult64(N, N, mod);
  191.         power >>= 1;
  192.     }
  193.     return res;
  194. }
  195.  
  196. template <typename T>
  197. bool isPrime(T N)
  198. {
  199.     if (N < 2 || N % 6 % 4 != 1)
  200.         return (N | 1) == 3;
  201.     T d = N - 1;
  202.     int s{};
  203.     while (!(d & 1))
  204.         d >>= 1, ++s;
  205.     for (const T &a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022})
  206.     {
  207.         T p = modPow(a % N, d, N), i = s;
  208.         while (p != 1 && p != N - 1 && a % N && i--)
  209.             p = mult64(p, p, N);
  210.         if (p != N - 1 && i != s)
  211.             return false;
  212.     }
  213.     return true;
  214. }
  215.  
  216. template <typename T>
  217. void primeFactorize(T N, map<T, T> &primeFactors)
  218. {
  219.     if (N == 1)
  220.         return;
  221.  
  222.     if (!isPrime(N))
  223.     {
  224.         T Y = Pollard_Brent(N);
  225.         primeFactorize(Y, primeFactors);
  226.         primeFactorize(N / Y, primeFactors);
  227.     }
  228.     else
  229.     {
  230.         primeFactors[N]++;
  231.         return;
  232.     }
  233. }
  234.  
  235. template <typename T>
  236. void getAllFactors(T N, vector<T> &factors)
  237. {
  238.     factors = {1};
  239.     map<T, T> freq;
  240.     primeFactorize(N, freq);
  241.  
  242.     for (auto &[p, cnt] : freq)
  243.     {
  244.         vector<T> temp;
  245.         T pw = 1;
  246.         for (int i = 0; i <= cnt; i++, pw *= p)
  247.         {
  248.             for (const T &f : factors)
  249.                 temp.push_back(f * pw);
  250.         }
  251.         factors.swap(temp);
  252.     }
  253.     sort(factors.begin(), factors.end());
  254. }
  255.  
  256. template <typename T>
  257. T countDivisors(T N)
  258. {
  259.     map<T, T> primeFactors;
  260.     primeFactorize(N, primeFactors);
  261.     T ans{1};
  262.     for (const auto &[p, exp] : primeFactors)
  263.         ans *= (exp + 1);
  264.     return ans;
  265. }
  266.  
  267. // GCC's implementation for I/O of 128-bit integers
  268. using int128 = signed __int128;
  269. using uint128 = unsigned __int128;
  270.  
  271. namespace int128_io
  272. {
  273.  
  274.     inline auto char_to_digit(int chr)
  275.     {
  276.         return static_cast<int>(isalpha(chr) ? 10 + tolower(chr) - 'a' : chr - '0');
  277.     }
  278.  
  279.     inline auto digit_to_char(int digit)
  280.     {
  281.         return static_cast<char>(digit > 9 ? 'a' + digit - 10 : '0' + digit);
  282.     }
  283.  
  284.     template <class integer>
  285.     inline auto to_int(const std::string &str, size_t *idx = nullptr, int base = 10)
  286.     {
  287.         size_t i = idx != nullptr ? *idx : 0;
  288.         const auto n = str.size();
  289.         const auto neg = str[i] == '-';
  290.         integer num = 0;
  291.         if (neg)
  292.             ++i;
  293.         while (i < n)
  294.             num *= base, num += char_to_digit(str[i++]);
  295.         if (idx != nullptr)
  296.             *idx = i;
  297.         return neg ? -num : num;
  298.     }
  299.  
  300.     template <class integer>
  301.     inline auto to_string(integer num, int base = 10)
  302.     {
  303.         const auto neg = num < 0;
  304.         std::string str;
  305.         if (neg)
  306.             num = -num;
  307.         do
  308.             str += digit_to_char(num % base), num /= base;
  309.         while (num > 0);
  310.         if (neg)
  311.             str += '-';
  312.         std::reverse(str.begin(), str.end());
  313.         return str;
  314.     }
  315.  
  316.     inline auto next_str(std::istream &stream)
  317.     {
  318.         std::string str;
  319.         stream >> str;
  320.         return str;
  321.     }
  322.  
  323.     template <class integer>
  324.     inline auto &read(std::istream &stream, integer &num)
  325.     {
  326.         num = to_int<integer>(next_str(stream));
  327.         return stream;
  328.     }
  329.  
  330.     template <class integer>
  331.     inline auto &write(std::ostream &stream, integer num) { return stream << to_string(num); }
  332. }
  333.  
  334. inline auto &operator>>(istream &stream, int128 &num) { return int128_io::read(stream, num); }
  335. inline auto &operator<<(ostream &stream, int128 num) { return int128_io::write(stream, num); }
  336. inline auto &operator>>(istream &stream, uint128 &num) { return int128_io::read(stream, num); }
  337. inline auto &operator<<(ostream &stream, uint128 num) { return int128_io::write(stream, num); }
  338.  
  339. int main()
  340. {
  341.     ios_base::sync_with_stdio(false);
  342.     cin.tie(nullptr);
  343. #ifndef ONLINE_JUDGE
  344.     freopen("input.txt", "r", stdin);
  345.     freopen("Output.txt", "w", stdout);
  346. #endif //! ONLINE_JUDGE
  347.     int t = 1;
  348.     // cin >> t;
  349.     while (t--)
  350.     {
  351.         int n, m, q;
  352.         cin >> n >> m >> q;
  353.         vector<ll> a(n), b(m);
  354.         ll sumA{}, sumB{};
  355.         bitset<400001> inA, inB;
  356.  
  357.         auto adjust = [&](ll N)
  358.         {
  359.             return N + 200000; // An offset to make all numbers positive
  360.         };
  361.         for (int i{}; i < n; i++)
  362.             cin >> a[i], sumA += a[i], inA[adjust(a[i])] = 1;
  363.  
  364.         for (int j{}; j < m; j++)
  365.             cin >> b[j], sumB += b[j], inB[adjust(b[j])] = 1;
  366.  
  367.         /*
  368.         (a1 * b1 + a1 * b2 ... a1 * bm)
  369.         (a2 * b1 + a2 * b2 ... a2 * bm)
  370.         ...
  371.         ...
  372.         Take ai as a common factor of each row, you get:
  373.         ai * sumB
  374.         Accumulate columns, you get:
  375.         a1 * sumB + a2 * sumB + a3 * sumB ... + an * sumB
  376.         Take sumB as a common factor, you get:
  377.         sum of matrix = B = sumA * sumB
  378.  
  379.         After any operation:
  380.         x = (sumA - ai) * (sumB - bj)
  381.         Factorize x, x = f * d
  382.         sumA - ai = f
  383.         ai = sumA - f
  384.         We need to check that (sumA - f) and (sumB - d) both exists
  385.         */
  386.         auto findInA = [&](ll x)
  387.         {
  388.             x = adjust(x);
  389.             return (x >= 0 && x <= 400000 && inA[x]);
  390.         };
  391.         auto findInB = [&](ll x)
  392.         {
  393.             x = adjust(x);
  394.             return (x >= 0 && x <= 400000 && inB[x]);
  395.         };
  396.         while (q--)
  397.         {
  398.             ll x;
  399.             cin >> x;
  400.             bool ok = false;
  401.             vector<u128> vc;
  402.             getAllFactors(u128(__abs(x)), vc);
  403.             for (const ll &d : vc)
  404.             {
  405.                 // x = (SumB − b[j]) ⋅ (SumA − a[i])
  406.                 ll f1 = d;
  407.                 ll f2 = x / d;
  408.                 ok |= (findInA(sumA - f1) && findInB(sumB - f2));
  409.                 ok |= (findInA(sumA - f2) && findInB(sumB - f1));
  410.                 f1 *= -1, f2 *= -1;
  411.                 ok |= (findInA(sumA - f1) && findInB(sumB - f2));
  412.                 ok |= (findInA(sumA - f2) && findInB(sumB - f1));
  413.             }
  414.             if (ok)
  415.                 cout << "YES\n";
  416.             else
  417.                 cout << "NO\n";
  418.         }
  419.     }
  420.     return 0;
  421. }
Success #stdin #stdout 0.01s 5324KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/UlimF8
language:
C++ (gcc 8.3)
created:
2 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!