Counting by divisors and harmonic-sum tricks

Why this matters when you build with AI

An AI tool asked for "the sum of floor(n / i) for every i from 1 to n" will often write one loop from 1 to n. That code is correct. For n = 10^6 it finishes in a few milliseconds. For n = 10^12 it runs 10^12 times, which takes far longer than a contest allows, so the judge stops it with a time limit exceeded verdict (TLE). On a small test the output looks right, and nothing in the code warns you. You catch this by timing the program at two sizes and checking how the time grows.

The same kind of gap appears with gcd. A loop that takes every pair of numbers from an array of 200000 values and computes their gcd does about 4 x 10^10 gcd calls. The array looks small in the code, and the cost is hidden in the nested loop.

This lesson gives three tools. The first counts the steps of a loop over multiples, using the harmonic sum. The second groups the values of floor(n / i) into blocks, so a sum over n terms takes about 2 times the square root of n steps. The third turns a sum of gcd values over pairs into a sum over divisor buckets, which is a count of how many numbers in the input each d divides.

The examples are C++17. We compiled them with g++ 15.2.0 using -std=c++17 -O2 -Wall -Wextra, and we tested each routine against a brute force on random small inputs. The timings below come from one machine, so your numbers will differ.

The idea: loops over multiples

How many steps a loop over multiples takes

Here is a loop that many contest solutions use.

for (int d = 1; d <= n; d++)
    for (int m = d; m <= n; m += d)
        // one step per multiple m of d

For each d, the inner loop visits the multiples d, 2d, 3d, and so on, up to n. There are floor(n / d) of them. So the total number of steps is the sum of floor(n / d) for d from 1 to n.

Take n = 10 and count by hand. For d = 1 the loop visits 10 numbers. For d = 2 it visits 5 (2, 4, 6, 8, 10). For d = 3 it visits 3 (3, 6, 9). For d = 4 it visits 2 (4, 8). For d = 5 it visits 2 (5, 10). For each of d = 6 to 10 it visits one number. The total is 10 + 5 + 3 + 2 + 2 + 5 = 27. The program below prints 27 for n = 10.

The sum of 1/d for d from 1 to n is called the harmonic number H_n. Since floor(n / d) is at most n / d, the step count is at most n times H_n. H_n grows like ln(n) plus about 0.577, so the nested loop costs about n ln(n) steps. It does not cost n squared. For n = 10^6 the exact count is 13,970,034 and n times H_n is about 14,392,727, so the bound is close.

This is why a sieve that marks multiples of every d (the method in the sieves lesson of this module) is fast. It is also why a loop over multiples is a good fallback when you need, for each d, the list of numbers it divides.

Here is the routine that counts the steps. We use it to check the formula, not to solve a problem.

// Number of (d, multiple) pairs visited by the nested loop: about n ln n.
ll harmonic_steps(ll n) {
    ll s = 0;
    for (ll d = 1; d <= n; d++)
        for (ll m = d; m <= n; m += d) s++;
    return s;
}

Divisor buckets

A bucket is a counter for one value d. The bucket c_d counts how many numbers in an input array are divisible by d. You fill all buckets with a loop over multiples. For each d, add up the frequency of every multiple m of d.

For the array [4, 6, 8, 9], the buckets are:

Building all buckets costs about M ln M steps when the values are at most M. That is the harmonic bound again.

Floor-division blocks

Why the values repeat

Look at floor(20 / i) for i from 1 to 20.

There are 20 terms but only 8 distinct values. The values repeat in runs. Inside one run every i gives the same quotient, and we call each run a block.

The general fact is this. For a fixed n, floor(n / i) takes at most about 2 times the square root of n distinct values. Here is why. If i is at most the square root of n, each i can give its own value, so there are at most sqrt(n) such values. If i is larger than the square root of n, then the quotient n / i is smaller than the square root of n, so there are at most sqrt(n) possible quotients. Together that is at most 2 sqrt(n) blocks.

Walking the blocks

Start a block at l. Let q = floor(n / l). Let r = floor(n / q). Then floor(n / i) equals q for every i from l to r, and r is the last such i. So the block is l to r, and the next block starts at r + 1.

Check this with n = 20. Start at l = 7. Then q = floor(20 / 7) = 2 and r = floor(20 / 2) = 10. The block is 7 to 10, as the table shows. The next start is 11. Then q = 1 and r = 20, so the block is 11 to 20.

The sum over all i is then the sum over blocks of q times the block length. The code below does that.

// Sum of floor(n / i) for i = 1..n, using blocks of equal quotient.
ll sum_floor_blocks(ll n) {
    ll total = 0;
    for (ll l = 1, r; l <= n; l = r + 1) {
        ll q = n / l;
        r = n / q;              // the block is l..r, and n / i == q for all of them
        total += q * (r - l + 1);
    }
    return total;
}

// Number of distinct values of floor(n / i) for i = 1..n.
int count_blocks(ll n) {
    int c = 0;
    for (ll l = 1, r; l <= n; l = r + 1) {
        ll q = n / l;
        r = n / q;
        c++;
    }
    return c;
}

The loop variable r must become r + 1 for the next block. If you write l = r by mistake, the same block repeats forever, and the program never ends. The symptom is a process that uses 100 percent of one core and never prints.

For the sums we used the 64-bit type long long (ll). The sum for n = 10^12 is 27785452449086, which is about 2.8 times 10^13. That does not fit in a 32-bit int.

We ran count_blocks on several values of n. The numbers below come from our run.

The block count stays close to 2 sqrt(n), as the argument says it must.

A sum identity that the blocks check

The sum of floor(n / i) over i from 1 to n equals the sum of d(k) over k from 1 to n, where d(k) is the number of divisors of k. Both sides count the pairs (a, b) of positive integers with a times b at most n. The left side counts them by a, and the right side counts them by k = a times b. We checked the identity by program for n = 1, 2, 10, 100, 1000 and 100000, and both sides matched at each value. For n = 10 both sides are 27.

This identity is useful when a problem gives you the divisor sum and you want a block-based count, or the other way round.

Gcd sums with divisor buckets

The identity

Let phi(d) be Euler's totient of d, the count of numbers from 1 to d that share no factor with d. For example phi(6) = 2, because only 1 and 5 share no factor with 6. (The totient lesson in this module proves the facts used here.)

The key identity is this. For any two positive integers a and b,

gcd(a, b) = sum of phi(d) over all d that divide both a and b

The reason is short. The numbers that divide both a and b are exactly the divisors of g = gcd(a, b). The sum of phi(d) over the divisors of g is g itself.

Check it on three pairs. For gcd(4, 6) = 2, the common divisors are 1 and 2, and phi(1) + phi(2) = 1 + 1 = 2. For gcd(4, 9) = 1, the only common divisor is 1, and phi(1) = 1. For gcd(6, 9) = 3, the common divisors are 1 and 3, and phi(1) + phi(3) = 1 + 2 = 3.

From pairs to buckets

Now sum the identity over all ordered pairs (i, j) of an array. Swap the order of the sums. For each d, the pairs that both have d as a common divisor are the pairs where both numbers are in bucket d. There are c_d times c_d such ordered pairs, including the pairs where i = j. So

sum over all ordered pairs (i, j) of gcd(a_i, a_j) = sum over d of phi(d) times c_d squared

Here is the calculation for [4, 6, 8, 9]. The bucket counts from the table above are c_1 = 4, c_2 = 3, c_3 = 2, c_4 = 2, c_5 = 0, c_6 = 1, c_7 = 0, c_8 = 1, c_9 = 1. The totients are phi(1) = 1, phi(2) = 1, phi(3) = 2, phi(4) = 2, phi(5) = 4, phi(6) = 2, phi(7) = 6, phi(8) = 4, phi(9) = 6. Terms with c_d = 0 vanish. The total is

1 x 16 + 1 x 9 + 2 x 4 + 2 x 4 + 2 x 1 + 4 x 1 + 6 x 1 = 53

The program in the worked example prints 53 for this array.

The code builds the buckets with the multiples loop, then sums phi(d) times c_d squared. The totients come from a sieve.

// Euler's totient for all values up to N (simple sieve).
vector<int> phi_table(int N) {
    vector<int> phi(N + 1);
    for (int i = 0; i <= N; i++) phi[i] = i;
    for (int p = 2; p <= N; p++)
        if (phi[p] == p)  // p is prime
            for (int m = p; m <= N; m += p) phi[m] -= phi[m] / p;
    return phi;
}

// Sum of gcd(a_i, a_j) over all ordered pairs i, j (including i == j).
// gcd(a, b) = sum of phi(d) over d dividing both a and b.
// So the total is sum over d of phi(d) * c_d^2, where c_d = how many a_i are divisible by d.
ll gcd_pair_sum(const vector<int>& a, int M) {
    vector<ll> freq(M + 1, 0);
    for (int x : a) freq[x]++;
    vector<int> phi = phi_table(M);
    ll total = 0;
    for (int d = 1; d <= M; d++) {
        ll cd = 0;                              // divisor bucket for d
        for (int m = d; m <= M; m += d) cd += freq[m];
        total += (ll)phi[d] * cd * cd;
    }
    return total;
}

The bucket counts are at most N, so c_d times c_d is at most N squared, which is 4 x 10^10 for N = 200000. Multiplied by phi(d), which is at most 10^6, the sum can reach about 4 x 10^16. That still fits in a 64-bit long long.

Worked example: the gcd sum of an array

The problem. You get N numbers a_1 to a_N. Print the sum of gcd(a_i, a_j) over all ordered pairs (i, j), including pairs where i equals j. Limits: 1 <= N <= 200000 and 1 <= a_i <= 10^6. The answer fits in a 64-bit integer.

The input format. The first line holds N. The second line holds the N numbers, separated by spaces.

The solution. Find the largest value M in the array. Fill the frequency table, build the buckets c_d with multiples, and sum phi(d) times c_d squared. The full program is below. It is the same code as the functions above, plus a main that reads the input.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
 
// Euler's totient for all values up to N (simple sieve).
vector<int> phi_table(int N) {
    vector<int> phi(N + 1);
    for (int i = 0; i <= N; i++) phi[i] = i;
    for (int p = 2; p <= N; p++)
        if (phi[p] == p)  // p is prime
            for (int m = p; m <= N; m += p) phi[m] -= phi[m] / p;
    return phi;
}
 
// Sum of gcd(a_i, a_j) over all ordered pairs i, j (including i == j).
// gcd(a, b) = sum of phi(d) over d dividing both a and b.
// So the total is sum over d of phi(d) * c_d^2, where c_d = how many a_i are divisible by d.
ll gcd_pair_sum(const vector<int>& a, int M) {
    vector<ll> freq(M + 1, 0);
    for (int x : a) freq[x]++;
    vector<int> phi = phi_table(M);
    ll total = 0;
    for (int d = 1; d <= M; d++) {
        ll cd = 0;                              // divisor bucket for d
        for (int m = d; m <= M; m += d) cd += freq[m];
        total += (ll)phi[d] * cd * cd;
    }
    return total;
}
 
int main() {
    int n;
    cin >> n;
    vector<int> a(n);
    int M = 1;
    for (int& x : a) {
        cin >> x;
        M = max(M, x);
    }
    cout << gcd_pair_sum(a, M) << "\n";
    return 0;
}

Real output on small inputs. We compiled the program and ran it on two inputs.

$ printf "4\n4 6 8 9\n" | ./ex11
53
$ printf "3\n2 3 6\n" | ./ex11
23

Check the first answer by hand. The diagonal pairs (i = j) give gcd(a, a) = a, so they add 4 + 6 + 8 + 9 = 27. The unordered pairs with i < j have gcds 2, 4, 1, 2, 3, 1 for the pairs (4,6), (4,8), (4,9), (6,8), (6,9), (8,9), which add to 13. Each unordered pair counts twice in the ordered sum, so the total is 27 + 2 x 13 = 53. The second answer is 2 + 3 + 6 = 11 on the diagonal, plus 2 x (1 + 2 + 3) = 12, for 23.

We also compared the program with a brute force on 2000 random arrays of up to 60 values, with values up to 300. All of them agreed. The brute force loops over every ordered pair and calls gcd directly.

Timing on one machine. For N = 200000 random values up to 10^6, the routine took about 180 milliseconds on our machine and printed 444631347449. The brute force would make 4 x 10^10 gcd calls, which is far beyond a contest limit.

Complexity

The brute force for the gcd sum takes N squared gcd calls. The bucket method replaces that with about M log M work, where M is the largest value and not the array length. So it is fast when the values are small, and it is slow when M is large, such as 10^9. Then the frequency table does not fit in memory. Use the prime-factor methods from the factorisation lesson instead.

Common mistakes

1. The sum overflows a 32-bit int. The sum of floor(n / i) for n = 10^12 is about 2.8 x 10^13. The symptom is a wrong value, sometimes negative. Signed overflow is undefined behaviour in C++, so the output can be anything. Use long long.
2. The block loop does not advance. Writing l = r instead of l = r + 1 repeats one block forever. The symptom is a program that never ends and never prints.
3. The bucket array is sized for the wrong bound. If freq has size N, but the values go up to M, then freq[x] writes past the end of the vector. The symptom is a crash or a result that changes from run to run. Size the array by the largest value.
4. The gcd loop is quadratic. A double loop over an array of 200000 values does 2 x 10^10 gcd calls for pairs with i < j. The symptom is a time limit exceeded verdict on large tests, even though the sample passes.
5. Only the sample was tested. The sample has N = 4. Overflow, a wrong bucket size and a slow loop show up only at large sizes. Test at the largest limit before you submit.
6. Trusting an AI summary of the complexity. An AI tool may say a loop is O(n log n) without checking the inner loop. Count the inner loop's steps yourself for one value of d, then sum over d.

How to catch these with an AI tool

Ask the tool for a brute force that is obviously correct, and run both on random inputs up to a few thousand. Then time the fast version at sizes 10^4, 10^5 and 10^6, and check that the time grows about as n log n, not as n squared. If the time multiplies by about 100 when n multiplies by 10, the code has a quadratic loop.

For the sums of floor(n / i), check the block count directly. For n = 10^6 the program should report 1999 blocks. A block count near n means the loop is not grouping values.

Hands-on exercise

Write the program below on your own machine, then test it.

Write a C++17 program. Input: a first line with Q (1 <= Q <= 100), then Q lines,
each with one integer n (1 <= n <= 10^12). For each n, print the sum of
floor(n / i) for i from 1 to n. Use the block method with 64-bit integers.
Also write a brute force loop over i. Compare the two on every n from 1 to 3000,
and on 200 random n up to 10^6. Then time your program on 100 copies of n = 10^12.
A good answer: the two methods agree everywhere, the value for n = 10 is 27,
and the 100 queries finish in about one or two seconds on your machine.
Check: the answer uses long long, the block loop ends at n, and the brute force
is a separate function, not a call to the fast one.

Your timings will differ from ours. The sums must match. On our machine one call with n = 10^12 took about 13 milliseconds, so 100 calls take about 1.3 seconds.

Checklist

• Before coding a loop over multiples, write the step count as a sum of floor(n / d).
• Use long long for any sum over n terms.
• In a block loop, set the block end as r = n / (n / l) and the next start as l = r + 1.
• Size each bucket array by the largest value that can occur.
• Compare every fast routine with a brute force on small inputs.
• Time the program at the largest allowed size before you trust it.

Summary

A loop over multiples costs about n ln n in total, because the number of multiples of d is floor(n / d) and the harmonic sum adds those up. The values of floor(n / i) form at most about 2 times the square root of n blocks, so a sum over n terms can be computed in that many steps. For a gcd sum over pairs, the identity gcd(a, b) = sum of phi(d) over common divisors turns the sum into a sum over buckets of counts, which costs about M log M when the values are at most M.

Practice

Practice 1: floor sum by hand

Find the sum of floor(20 / i) for i from 1 to 20 without running code.

Hint: group the i values into the blocks of the table in the floor-division section, and multiply each quotient by the block length.

Show the answer

The blocks have quotients 20, 10, 6, 5, 4, 3, 2 and 1, with lengths 1, 1, 1, 1, 1, 1, 4 and 10. The sum is 20 + 10 + 6 + 5 + 4 + 3 + 2 x 4 + 1 x 10 = 66.

Practice 2: count the blocks

Run count_blocks for n = 1000000 and report the number of blocks. Check it against 2 times sqrt(n).

Hint: the function prints nothing by itself, so call it from a small main that prints its result.

Show the answer

Our run printed 1999 blocks. The bound 2 times sqrt(10^6) is 2000, and the count is below it, as the argument predicts.

Practice 3: buckets for an array

For the array [2, 4, 6, 9], find c_d for every d from 1 to 9, where c_d is the number of array values divisible by d.

Hint: for each d, test each value with the modulo operator, or add up the frequencies of the multiples of d.

Show the answer

The counts are c_1 = 4, c_2 = 3, c_3 = 2, c_4 = 1, c_5 = 0, c_6 = 1, c_7 = 0, c_8 = 0 and c_9 = 1. A short Python check printed the same list.

Practice 4: gcd sum by buckets

Find the sum of gcd(a_i, a_j) over all ordered pairs of the array [2, 4, 6, 9], including the pairs where i equals j.

Hint: the diagonal adds the values themselves. Then find the unordered gcd sum for i less than j and double it.

Show the answer

The diagonal is 2 + 4 + 6 + 9 = 21. The unordered pairs give gcd 2 for (2, 4), 2 for (2, 6), 1 for (2, 9), 2 for (4, 6), 1 for (4, 9) and 3 for (6, 9), which add to 11. The total is 21 + 2 x 11 = 43. The program gcd_pair_sum with the same array gives the same total.

Practice 5: why the multiples loop is n log n

Explain, in two or three sentences, why the nested loop over d and its multiples costs about n ln n steps and not n squared.

Hint: for a fixed d, count how many multiples of d are at most n. Then sum that count over d.

Show the answer

For a fixed d there are floor(n / d) multiples of d up to n. Summing floor(n / d) over d from 1 to n gives at most n times H_n, where H_n is the sum of 1/d. Since H_n is about ln(n), the total is about n ln(n). The steps shrink as d grows, because large d has few multiples.