Showing posts with label Number Theory. Show all posts
Showing posts with label Number Theory. Show all posts

Wednesday, October 28, 2026

GCD and LCM in C++: Euclid's Algorithm Explained

The greatest common divisor (GCD) is one of the oldest algorithms in mathematics, and it still shows up constantly in contest problems: simplifying fractions, synchronising cycles, and in clever tricks that turn a hard-looking problem into one line. This post covers how Euclid's algorithm works, why it is fast, and the LCM overflow trap.

Definitions

gcd(a, b) is the largest number that divides both a and b. lcm(a, b) is the smallest positive number that both a and b divide. For example gcd(12, 18) = 6 and lcm(12, 18) = 36. They are linked by a simple identity:

gcd(a, b) × lcm(a, b) = a × b

Euclid's algorithm

The key fact is gcd(a, b) = gcd(b, a mod b). Any number that divides both a and b also divides a − kb for any k, and a mod b is exactly such a number. The reverse also holds, so both pairs have the same common divisors and therefore the same greatest one.

Keep replacing the pair until the second number is 0. Then gcd(g, 0) = g.

Example for gcd(252, 105):

gcd(252, 105) -> 252 mod 105 = 42
gcd(105, 42)  -> 105 mod 42  = 21
gcd(42, 21)   -> 42 mod 21   = 0
gcd(21, 0)    = 21

Why it is fast

After two steps, the larger number always drops to less than half of what it was. So the number of steps is O(log min(a, b)), at most about 90 steps for numbers up to 1018. The slowest inputs are consecutive Fibonacci numbers, which is a fun fact worth remembering.

Code (C++)

This program reads n numbers and prints the GCD and LCM of all of them.

#include <bits/stdc++.h>
using namespace std;

long long gcd_ll(long long a, long long b) {
    while (b != 0) {
        long long r = a % b;
        a = b;
        b = r;
    }
    return a;
}

int main() {
    int n;
    cin >> n;
    vector<long long> a(n);
    for (auto &x : a) cin >> x;

    long long g = 0, l = 1;
    for (long long x : a) {
        g = gcd_ll(g, x);            // gcd(0, x) = x, so 0 is a safe start
        l = l / gcd_ll(l, x) * x;    // divide first to keep numbers small
    }
    cout << g << ' ' << l << '\n';
    return 0;
}

Since C++17 you can also use std::gcd and std::lcm from <numeric>. Writing it yourself once is still worth doing so you know what is inside.

The LCM overflow trap

The obvious formula a * b / gcd(a, b) multiplies first. For two numbers near 1010 the product is 1020, which overflows long long even though the LCM itself might be small. Dividing first, a / gcd(a, b) * b, avoids the giant intermediate value. The final LCM can still be too large to store if the numbers share few factors, so check the problem's limits.

Handy properties

  • gcd(a, 0) = a, which makes 0 a safe starting value when taking the GCD of a list.
  • The GCD of a list can be computed one element at a time: gcd(a, b, c) = gcd(gcd(a, b), c).
  • gcd(a, b) = gcd(a, b − a). This lets you rewrite expressions like gcd(a1 + x, a2 + x, …) in terms of differences that do not depend on x, which is the whole trick behind Codeforces 1458A – Row GCD.
  • For factorials, gcd(A!, B!) is just min(A, B)!, because the smaller factorial divides the larger one. That is all Codeforces 822A – I'm bored with life asks for.

Practice problems

For more number theory, see the Lucky Division walkthrough.

Binary Exponentiation: Compute a^b mod m in O(log b)

Many problems ask for an answer “modulo 109 + 7”, and many of those need a large power such as 2n or ab with b up to 1018. Multiplying a by itself b times is hopeless. Binary exponentiation does it in about 60 multiplications.

The idea: square instead of multiply

Notice that a8 = ((a2)2)2. That is only three multiplications instead of seven. Repeated squaring gives us a, a2, a4, a8, a16, … very cheaply.

Any exponent can be written as a sum of powers of two, which is exactly its binary representation. For example 13 = 8 + 4 + 1, which is 1101 in binary. So

a13 = a8 × a4 × a1

We walk through the bits of b from the lowest. At each step, the current a holds a2k. If the bit is 1, multiply it into the result. Then square a and shift b right.

Tracing a = 3, b = 13

b (binary)lowest bitcurrent aresult
1101131 × 3 = 3
110093
111813 × 81 = 243
116561243 × 6561 = 1 594 323

And 313 is indeed 1 594 323, found with a handful of multiplications.

Code (C++)

#include <bits/stdc++.h>
using namespace std;

const long long MOD = 1'000'000'007;

// Computes (a^b) mod m in O(log b) multiplications.
long long power(long long a, long long b, long long m) {
    long long result = 1 % m;
    a %= m;
    while (b > 0) {
        if (b & 1) result = result * a % m;   // this bit of b is set
        a = a * a % m;                        // a, a^2, a^4, a^8, ...
        b >>= 1;
    }
    return result;
}

int main() {
    int q;
    cin >> q;
    while (q--) {
        long long a, b;
        cin >> a >> b;
        cout << power(a, b, MOD) << '\n';
    }
    return 0;
}

The input format matches CSES 1095 – Exponentiation: a number of queries, then a and b on each line.

Complexity

The loop runs once per bit of b, so O(log b) multiplications. For b up to 1018 that is at most 60 iterations.

Details that cause Wrong Answer

  • Overflow. Both factors are below 109 + 7, so their product is below about 1018, which fits in long long. With int it overflows. If the modulus is close to 1018, even long long is not enough and you need __int128.
  • Reduce a first. If a itself is larger than the modulus, take a %= m before the loop.
  • b = 0. Anything to the power 0 is 1, including 00 by the usual convention in these problems. The code starts result at 1 and handles it without a special case.

Where it shows up

  • Modular inverse. When the modulus p is prime, Fermat's little theorem gives a−1 ≡ ap−2 (mod p). That is how you “divide” under a modulus, for example when computing binomial coefficients.
  • Huge exponents. To compute abc mod p, reduce the exponent bc modulo p − 1 first (again by Fermat), then use binary exponentiation twice.
  • Matrix exponentiation. The same loop works with matrices instead of numbers, which computes the n-th Fibonacci number in O(log n).

Practice problems

Sieve of Eratosthenes: Find All Primes up to N Fast (C++)

Checking whether one number is prime is easy: try dividing it by everything up to its square root. But when a problem needs to know about every number up to 106 or 107, repeating that check for each number is too slow. The Sieve of Eratosthenes, which is more than two thousand years old, finds all of them at once.

The idea

Write down all numbers from 2 to n, and assume they are all prime. Then:

  1. Take the smallest number that is still marked prime. Call it p. It really is prime, because no smaller number divided it.
  2. Cross out every multiple of p: they have p as a factor, so they are composite.
  3. Move on to the next number that is still marked and repeat.

Whatever survives is prime.

Example with n = 30

  • p = 2: cross out 4, 6, 8, …, 30.
  • p = 3: cross out 9, 15, 21, 27 (12, 18, 24 and 30 were already gone).
  • p = 5: cross out 25 (30 was already gone).
  • p = 7: 7 × 7 = 49 > 30, so we can stop.

The survivors are 2, 3, 5, 7, 11, 13, 17, 19, 23, 29: ten primes.

Two small optimisations that matter

Start crossing at p × p. Any smaller multiple k × p with k < p has a prime factor smaller than p, so it was already crossed out earlier. For p = 5 above, 10, 15 and 20 were already gone, so we start at 25.

Stop the outer loop at √n. Once p × p > n, there is nothing left to cross out. Every composite number up to n has a prime factor at most √n, so it has already been handled.

Code (C++)

#include <bits/stdc++.h>
using namespace std;

int main() {
    int n;
    cin >> n;

    // is_prime[x] is true if x is prime, for 0 <= x <= n
    vector<bool> is_prime(n + 1, true);
    is_prime[0] = false;
    if (n >= 1) is_prime[1] = false;

    for (long long p = 2; p * p <= n; p++) {
        if (!is_prime[p]) continue;
        for (long long m = p * p; m <= n; m += p)
            is_prime[m] = false;
    }

    int count = 0;
    for (int x = 2; x <= n; x++)
        if (is_prime[x]) count++;
    cout << count << '\n';
    return 0;
}

The program prints how many primes there are up to n. For n = 106 it prints 78 498, and for n = 107 it prints 664 579. Once is_prime is built, each “is x prime?” question is a single array lookup.

Complexity

For each prime p, the inner loop does about n / p steps. Summing n/2 + n/3 + n/5 + … over primes gives O(n log log n), which is almost linear in practice. Memory is O(n). vector<bool> stores one bit per entry, so a sieve up to 107 needs only about 1.25 MB.

Common mistakes

  • Forgetting that 0 and 1 are not prime.
  • Writing p * p <= n with int when n is close to the int limit. The product can overflow, so use long long for p, as in the code above.
  • Re-running the sieve for every test case. Build it once, up to the largest value you will need.

A useful variation: smallest prime factor

Instead of a true/false flag, store for every number its smallest prime factor. Then you can factorise any number up to n in O(log n) by repeatedly dividing by its smallest prime factor. This is the standard way to factorise many numbers quickly.

Practice problems

For a gentler number theory warm-up, see Lucky Division (122A).

Monday, August 6, 2018

Codeforces 122A Lucky Division Solution in C (Generating Lucky Numbers)

Lucky Division (122A) is a small number theory problem, and it is a good place to practice one habit: when the input is small, generate everything and check it. It also has a trap. If you list the lucky numbers by hand, it is easy to miss some, and an earlier version of this very post did exactly that.

The problem in short

A positive integer is lucky if its decimal digits are only 4s and 7s (for example 4, 47, 744). A number is almost lucky if it is divisible by some lucky number. Given n (1 ≤ n ≤ 1000), print YES if it is almost lucky, otherwise NO.

Full statement: Codeforces 122A – Lucky Division.

Which lucky numbers do we need?

A divisor of n can never be larger than n, and n ≤ 1000. So we need every lucky number up to 1000. With 1 digit there are 2 (4, 7), with 2 digits there are 4 (44, 47, 74, 77) and with 3 digits there are 8, giving 14 in total:

4 7 44 47 74 77 444 447 474 477 744 747 774 777

The old version of this post skipped 77, 744, 747, 774 and 777. Some of those are harmless to skip, since 77 and 777 are multiples of 7 and 744 is a multiple of 4. But 747 and 774 are not divisible by any smaller lucky number, so the old code wrongly printed NO for them. That is why it is safer to list the full set than to reason about which ones you can drop.

Approach

Loop over the 14 lucky numbers. If any of them divides n exactly, print YES. If none do, print NO.

Examples

  • n = 47: 47 is itself lucky, so it divides itself → YES.
  • n = 16: 16 is divisible by 4 → YES.
  • n = 78: 78 = 2 × 3 × 13. None of the 14 lucky numbers divide it → NO.

Solution in C

#include <stdio.h>

/* All lucky numbers that are <= 1000 (the problem's limit). */
int lucky[] = {4, 7, 44, 47, 74, 77, 444, 447, 474, 477, 744, 747, 774, 777};

int main(void) {
    int n;
    scanf("%d", &n);

    int count = sizeof(lucky) / sizeof(lucky[0]);
    for (int i = 0; i < count; i++) {
        if (n % lucky[i] == 0) {
            printf("YES\n");
            return 0;
        }
    }
    printf("NO\n");
    return 0;
}

Going further: generating the list

If the limit were 109 instead of 1000, you would not want to type the list. You can generate lucky numbers recursively: start from 0, and from any number x produce 10x + 4 and 10x + 7, stopping once the value passes the limit. There are only 2 + 4 + … + 29 = 1022 lucky numbers below 109, so even that is tiny.

Complexity

14 divisions per input: O(1) time and memory.

Practice next

For another “check a small fixed set” idea, see 282A Bit++. For a math shortcut problem, try 4A Watermelon.