MLTDVD - Editorial


Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4

Author: beevo
Tester & Editorialist: iceknight1093




Elementary Number Theory, Prime factorization using a sieve


You have an array A containing N integers upto 10^7.
Define a sequence of arrays A^{(k)}, where:

  • A^{(0)} = A
  • A^{(k)} = \frac{L}{G} \cdot A^{(k-1)} for k \geq 1, where G and L are respectively the gcd and lcm of A^{(k-1)}

Answer Q queries; each giving you an integer k and asking you to find the gcd and lcm of A^{(k)}, modulo 10^9 + 7.


Let G_k and L_k denote the gcd and lcm of A^{(k)}.

Let’s first find the answers for K = 0, i.e, G_0 and L_0.
Finding G_0 is trivial by just iterating across the array; however, L_0 can be extremely large — computing it directly it impossible.

However, notice that A_i \leq 10^7 means that even if L_0 is very large, all its prime factors will be \leq 10^7.
This means we can find L_0 in terms of its prime factorization.

Recall that if we have two integers x and y such that

\begin{aligned} x = p_1^{a_1}p_2^{a_2}\ldots p_k ^ {a_k} \\ y = p_1^{b_1}p_2^{b_2}\ldots p_k ^ {b_k} \\ \end{aligned}

then \text{lcm}(x, y) = p_1^{\max(a_1, b_1)}p_2^{\max(a_2, b_2)}\ldots p_k^{\max(a_k, b_k)}

Of course, this generalizes to the \text{lcm} of more than two integers as well: the power of each prime p_i in the \text{lcm} will be the maximum power of p_i present in one of the numbers.

Using this idea, let’s compute L_0.
Consider some prime p \leq 10^7. To know the power of p in L_0, we want to know the maximum power of p across all A_i.

In particular, we have the following algorithm:

  • Let \text{mx}_p denote the maximum power of p present in some A_i.
    Initially, \text{mx}_p = 0 for all p.
  • Then, for each i from 1 to N:
    • Find all prime factors of A_i; let them be p_1, p_2, \ldots, p_k with respective powers a_1, a_2, \ldots, a_k.
    • For each 1 \leq i \leq k, set \text{mx}_{p_i} \gets \max(\text{mx}_{p_i}, a_i)

At the end of this process, \text{mx}_p will be the value we want, for every p.

For this to run quickly, we need to be able to quickly prime-factorize the A_i values.
This can be done using a slight modification to the sieve of Eratosthenes, storing a prime divisor of each integer instead of just a boolean value.
An overview of this can be seen in this blog.

Any integer x has \leq \log_2 x prime factors (since each prime factor is at least 2), so the entire process above takes \mathcal{O}(N\log 10^7) time.

Once this is done, L_0 is the product of p^{\text{mx}_p} across all primes \leq p.

For k \gt 0, notice that for any integers x, a_1, a_2, \ldots, a_k we have:

  • \gcd(x\cdot a_1, x\cdot a_2, \ldots, x\cdot a_k) = x\cdot \gcd(a_1, a_2, \ldots, a_k)
  • \text{lcm}(x\cdot a_1, x\cdot a_2, \ldots, x\cdot a_k) = x\cdot \text{lcm}(a_1, a_2, \ldots, a_k)

That is, multiplying all the numbers of a sequence by a constant multiplies the gcd/lcm of the sequence by the same constant.
Proofs of these can be found here.

Further, note that if the gcd and lcm are multiplied by the same constant, their ratio remains the same!
In particular, notice that A^{(k)} is obtained from A^{(k-1)} via multiplication by a constant.

This means that \displaystyle \frac{L_k}{G_k} = \frac{L_0}{G_0} for all k \geq 0, since each array is obtained in succession via multiplication with a constant.

However, \displaystyle\frac{L_k}{G_k} is also the constant we multiply A^{(k)} with to get A^{(k+1)}.
So, we simply have

\begin{aligned} G_k = G_0 \cdot \left( \frac{L_0}{G_0}\right) ^k \\ L_k = L_0 \cdot \left( \frac{L_0}{G_0}\right) ^k \end{aligned}

which can both be found in \mathcal{O}(\log k) time using binary exponentiation, since we already know G_0 and L_0.


\mathcal{O}(10^7 \log\log {10^7} + N\log 10^7 + Q\log {MOD}) per testcase.


Author's code (C++)
#include <bits/stdc++.h>

#define el '\n'

typedef long long ll;
typedef long double ld;

#define Beevo ios_base::sync_with_stdio(0); cin.tie(0); cout.tie(0);

using namespace std;

const int N = 1e7 + 5, M = 1e9 + 7;

int spf[N], maxFreq[N];

int mul(int a, int b) {
    return 1LL * a * b % M;

int modPow(int b, int p) {
    if (!p)
        return 1;

    int x = modPow(b, p / 2);

    return p % 2 == 0 ? mul(x, x) : mul(b, mul(x, x));

int modInvFer(int n) {
    return modPow(n, M - 2);

void sieve() {
    for (int i = 1; i < N; i++)
        spf[i] = i;

    for (int i = 2; i * i < N; i++) {
        if (spf[i] == i) {
            for (int j = i * 2; j < N; j += i)
                spf[j] = min(spf[j], i);

void factorize(int n) {
    int pf, freq;
    while (n > 1) {
        pf = spf[n], freq = 0;

        while (n % pf == 0)
            freq++, n /= pf;

        maxFreq[pf] = max(maxFreq[pf], freq);

void testCase() {

    int n;
    cin >> n;

    int g = 0;
    for (int i = 0, a; i < n; i++)
        cin >> a, g = __gcd(g, a), factorize(a);

    int l = 1;
    for (int i = 1; i < N; i++)
        l = mul(l, modPow(i, maxFreq[i]));

    int q;
    cin >> q;

    int k, p, d = mul(l, modInvFer(g));
    while (q--) {
        cin >> k;

        p = modPow(d, k);

        cout << mul(g, p) << ' ' << mul(l, p) << el;

signed main() {

    int t = 1;
//    cin >> t;

    while (t--)
Editorialist's code (Python)
from math import gcd
lim = 10**7 + 10
mod = 10**9 + 7
primefac = [1]*lim
for i in range(2, lim):
    if primefac[i] == 1:
        for j in range(i, lim, i): primefac[j] = i

n = int(input())
a = list(map(int, input().split()))
g, l = a[0], 1
mxpow = [0]*lim
for x in a:
    g = gcd(x, g)
    while x > 1:
        p = primefac[x]
        ct = 0
        while x%p == 0:
            x //= p
            ct += 1
        mxpow[p] = max(mxpow[p], ct)
for x in range(2, lim):
    if primefac[x] != x: continue
    l *= pow(x, mxpow[x], mod)
    l %= mod

for _ in range(int(input())):
    k = int(input())
    mul = pow(l, k, mod) * pow(g, k*(mod-2), mod) % mod
    print(g * mul % mod, l * mul % mod)

Time complexity should be \mathcal{O}(\dots+Q\log(\max k)), not \mathcal{O}(\dots+Q\log MOD).

Depends on your implementation.
If you use modulo inverse to divide for each query (which my code does), it’s indeed Q \log{MOD}.