Skip to main content

Modular Arithmetic, Fast Power and Counting

Editorial

A billion multiplications, or thirty

Count Good Numbers asks how many strings of n digits have an even digit at every even position (counting from 0) and a prime digit (2, 3, 5 or 7) at every odd one. Each even position has 5 choices and each odd one 4, so the count is 5^⌈n/2⌉ × 4^⌊n/2⌋, returned mod 10⁹ + 7, for n up to 10¹⁵. Multiplying it out one factor at a time, timed with a power of 3:

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

const long long M = 1000000007;

int main() {
    long long n = 1000000000;
    auto start = chrono::steady_clock::now();
    long long result = 1;
    for (long long k = 0; k < n; k++) result = result * 3 % M;   // one multiplication per step
    double s = chrono::duration<double>(chrono::steady_clock::now() - start).count();
    printf("3^%lld mod M = %lld in %.2fs\n", n, result, s);
}

// 3^1000000000 mod M = 235939645 in 3.91s
import time

M = 10 ** 9 + 7
n = 10 ** 7                           # a billion would take minutes here
start = time.perf_counter()
result = 1
for _ in range(n):
    result = result * 3 % M           # one multiplication per step
print(f"3^{n} mod M = {result} in {time.perf_counter() - start:.2f}s")

# 3^10000000 mod M = 769346453 in 1.64s
public class Main {
    public static void main(String[] args) {
        final long M = 1000000007;
        long n = 1000000000;
        long start = System.nanoTime();
        long result = 1;
        for (long k = 0; k < n; k++) result = result * 3 % M;   // one multiplication per step
        System.out.printf("3^%d mod M = %d in %.2fs%n", n, result, (System.nanoTime() - start) / 1e9);
    }
}

// 3^1000000000 mod M = 235939645 in 4.45s
const M = 1000000007;
const n = 1000000000;
const start = Date.now();
let result = 1;
for (let k = 0; k < n; k++) result = result * 3 % M;   // one multiplication per step
console.log(`3^${n} mod M = ${result} in ${((Date.now() - start) / 1000).toFixed(2)}s`);

// 3^1000000000 mod M = 235939645 in 13.49s

A billion multiplications took 3.91 s in C++, 4.45 s in Java and 13.49 s in JavaScript; Python managed ten million in 1.64 s. At that C++ rate, 10¹⁵ multiplications would take about 45 days.

Two ideas make it instant. Keep only the remainder: the true count has hundreds of trillions of digits, but its remainder mod 10⁹ + 7 fits in an int, and you can reduce after every operation. And square instead of multiply: x² times x² is x⁴, which times itself is x⁸, so about 50 squarings reach an exponent of 10¹⁵. This lesson builds both, then counts with them: Pascal's rule, what C(n, k) counts, and C(n, k) mod a prime.

A dashed band, the bits of 13, from the top: four tiles, 1, 1, 0 and 1, the ones in orange, with the words start at 3¹, square, then × 3, square, and square, then × 3 under them; under those, four blue tiles joined by arrows, 3^1, 3^3, 3^6 and 3^13, with their values 3, 27, 729 and 1,594,323 and their exponents 1, 2 × 1 + 1 = 3, 2 × 3 = 6 and 2 × 6 + 1 = 13. Below, a ticked card, 5 multiplications: 3 squarings, one per bit after the first, and 2 times × 3, one per 1 after the first; and a crossed card, 12 multiplications: 3 × 3 × 3 × … one factor at a time, n − 1 of them, 12 for 13, a billion for 10⁹. Focal line: every bit: square; every 1 bit: one more × x. Title: square, and multiply on a 1.

Keeping only the remainder

Addition, subtraction and multiplication all survive taking remainders first: (a + b) mod m is ((a mod m) + (b mod m)) mod m, and the same with × and −. So a long computation can take the remainder after every step and never grow:

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

int main() {
    printf("(15 + 23) %% 7 = %d, and from the remainders %d\n", (15 + 23) % 7, (15 % 7 + 23 % 7) % 7);
    printf("(15 * 23) %% 7 = %d, and from the remainders %d\n", (15 * 23) % 7, (15 % 7) * (23 % 7) % 7);
    printf("(5 - 8) %% 7 = %d; add 7 first: %d\n", (5 - 8) % 7, ((5 - 8) % 7 + 7) % 7);
}

// (15 + 23) % 7 = 3, and from the remainders 3
// (15 * 23) % 7 = 2, and from the remainders 2
// (5 - 8) % 7 = -3; add 7 first: 4
print(f"(15 + 23) % 7 = {(15 + 23) % 7}, and from the remainders {(15 % 7 + 23 % 7) % 7}")
print(f"(15 * 23) % 7 = {(15 * 23) % 7}, and from the remainders {(15 % 7) * (23 % 7) % 7}")
print(f"(5 - 8) % 7 = {(5 - 8) % 7}; add 7 first: {((5 - 8) % 7 + 7) % 7}")

# (15 + 23) % 7 = 3, and from the remainders 3
# (15 * 23) % 7 = 2, and from the remainders 2
# (5 - 8) % 7 = 4; add 7 first: 4
public class Main {
    public static void main(String[] args) {
        System.out.println("(15 + 23) % 7 = " + (15 + 23) % 7 + ", and from the remainders " + (15 % 7 + 23 % 7) % 7);
        System.out.println("(15 * 23) % 7 = " + (15 * 23) % 7 + ", and from the remainders " + (15 % 7) * (23 % 7) % 7);
        System.out.println("(5 - 8) % 7 = " + (5 - 8) % 7 + "; add 7 first: " + ((5 - 8) % 7 + 7) % 7);
    }
}

// (15 + 23) % 7 = 3, and from the remainders 3
// (15 * 23) % 7 = 2, and from the remainders 2
// (5 - 8) % 7 = -3; add 7 first: 4
console.log(`(15 + 23) % 7 = ${(15 + 23) % 7}, and from the remainders ${(15 % 7 + 23 % 7) % 7}`);
console.log(`(15 * 23) % 7 = ${(15 * 23) % 7}, and from the remainders ${(15 % 7) * (23 % 7) % 7}`);
console.log(`(5 - 8) % 7 = ${(5 - 8) % 7}; add 7 first: ${((5 - 8) % 7 + 7) % 7}`);

// (15 + 23) % 7 = 3, and from the remainders 3
// (15 * 23) % 7 = 2, and from the remainders 2
// (5 - 8) % 7 = -3; add 7 first: 4
Language(5 − 8) % 7((5 − 8) % 7 + 7) % 7
C++−34
Python44
Java−34
JavaScript−34

Subtraction is the one to watch. In C++, Java and JavaScript the remainder of a negative number is negative (Digits and Number Manipulation showed why), so a difference of two remainders needs m added back before the last mod. Python's % is never negative for a positive m. Division doesn't survive at all: it's done by multiplying with an inverse, below.

Why 10⁹ + 7. It's prime, so every remainder except 0 has an inverse. It's below 2³¹, so a remainder fits a 32-bit int. And two remainders multiply to less than (10⁹ + 7)², about 1.0 × 10¹⁸, below 2⁶³ ≈ 9.2 × 10¹⁸, so their product fits a 64-bit integer before it's reduced.

Multiplying two remainders

That last point decides the types. In C++ and Java, multiply remainders as 64-bit values (long long, long); two ints multiplied as ints overflow first. In JavaScript every number is a double, exact only up to 2⁵³ ≈ 9.0 × 10¹⁵, and a product of two remainders can reach 10¹⁸:

A dashed panel, the largest whole number each type holds exactly, on a log scale from 10⁹: three bars, a 32-bit int ending at 2³¹ in purple, a double ending at 2⁵³ in orange and a 64-bit integer ending at 2⁶³ in green; a red dashed line at (10⁹ + 6)² ≈ 1.0 × 10¹⁸ crosses them, with a cross on the int and double rows and a tick on the 64-bit row, and 10⁹: one remainder at the left. Below, a ticked card, 64-bit integers: 64, C++ long long, Java long, BigInt, the product, then the remainder; and a crossed card, JavaScript numbers: 63, the product was rounded to a double before the % ever ran. Focal line: multiply remainders in 64 bits, never in a double. Title: where two remainders fit.
#include <bits/stdc++.h>
using namespace std;

const long long M = 1000000007;

int main() {
    long long a = 999999999, b = 999999999;   // two remainders, multiplied as 64-bit values
    printf("%lld\n", a * b % M);
}

// 64
M = 10 ** 9 + 7
a = b = 999999999
print(a * b % M)                      # Python's ints are exact at any size

# 64
public class Main {
    public static void main(String[] args) {
        final long M = 1000000007;
        long a = 999999999, b = 999999999;   // two remainders, multiplied as 64-bit values
        System.out.println(a * b % M);
    }
}

// 64
const M = 1000000007;
const a = 999999999, b = 999999999;
console.log("as numbers:", a * b % M);                         // a double holds integers exactly only up to 2^53
console.log("as BigInt:", Number(BigInt(a) * BigInt(b) % BigInt(M)));

// as numbers: 63
// as BigInt: 64

999,999,999 × 999,999,999 mod (10⁹ + 7) is 64. JavaScript's plain numbers gave 63: the product, 999,999,998,000,000,001, was rounded to the nearest double before the % ever ran. BigInt is exact. Multiplied as 32-bit ints in C++, the same product wrapped and gave 808,348,673. Python's ints are exact at any size, so Python only has to remember to reduce.

Fast power

x^n is (x^(n/2))² when n is even, and x × x^(n−1) when it's odd, so every squaring halves what's left of the exponent. The loop below does it from the low end of n's binary digits: x, x², x⁴, x⁸ and so on are squared up one per bit, and each power whose bit of n is 1 is multiplied into the result. 13 is 1101 in binary, so x¹³ = x⁸ × x⁴ × x¹:

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

const long long M = 1000000007;

long long power(long long x, long long n) {   // x^n mod M by squaring: O(log n)
    long long result = 1;
    x %= M;
    while (n > 0) {
        if (n & 1) result = result * x % M;   // this bit of n is 1: multiply this power of x in
        x = x * x % M;                        // x, x², x⁴, x⁸, ...: square for the next bit
        n >>= 1;
    }
    return result;
}

int main() {
    printf("%lld %lld %lld\n", power(2, 10), power(3, 1000000000), power(2, 1000000000000000000LL));
}

// 1024 235939645 719476260
M = 10 ** 9 + 7


def power(x, n):                      # x^n mod M by squaring: O(log n)
    result = 1
    x %= M
    while n > 0:
        if n & 1:                     # this bit of n is 1: multiply this power of x in
            result = result * x % M
        x = x * x % M                 # x, x², x⁴, x⁸, ...: square for the next bit
        n >>= 1
    return result


print(power(2, 10), power(3, 10 ** 9), power(2, 10 ** 18))   # Python's own pow(x, n, M) does the same

# 1024 235939645 719476260
public class Main {
    static final long M = 1000000007;

    static long power(long x, long n) {   // x^n mod M by squaring: O(log n)
        long result = 1;
        x %= M;
        while (n > 0) {
            if ((n & 1) == 1) result = result * x % M;   // this bit of n is 1: multiply this power of x in
            x = x * x % M;                             // x, x², x⁴, x⁸, ...: square for the next bit
            n >>= 1;
        }
        return result;
    }

    public static void main(String[] args) {
        System.out.println(power(2, 10) + " " + power(3, 1000000000L) + " " + power(2, 1000000000000000000L));
    }
}

// 1024 235939645 719476260
const M = 1000000007n;

function power(x, n) {                 // x^n mod M by squaring: O(log n); BigInt, since x * x passes 2^53
    let result = 1n;
    x %= M;
    while (n > 0n) {
        if (n & 1n) result = result * x % M;   // this bit of n is 1: multiply this power of x in
        x = x * x % M;                         // x, x², x⁴, x⁸, ...: square for the next bit
        n >>= 1n;
    }
    return result;
}

console.log(`${power(2n, 10n)} ${power(3n, 1000000000n)} ${power(2n, 1000000000000000000n)}`);

// 1024 235939645 719476260

3^1,000,000,000 took one squaring per bit of n, which has 30 bits, and one multiplication per 1 bit, 13 of them: under 60 multiplications in all, where the loop above did a billion. O(log n). 2^(10¹⁸) took one squaring for each of its exponent's 60 bits. Every multiplication is reduced mod 10⁹ + 7, and JavaScript uses BigInt for it.

Count Good Numbers is two fast powers multiplied together:

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

const long long M = 1000000007;

long long power(long long x, long long n) {   // x^n mod M by squaring: O(log n)
    long long result = 1;
    x %= M;
    while (n > 0) {
        if (n & 1) result = result * x % M;   // this bit of n is 1: multiply this power of x in
        x = x * x % M;                        // x, x², x⁴, x⁸, ...: square for the next bit
        n >>= 1;
    }
    return result;
}

long long countGoodNumbers(long long n) {   // 5 even digits at even positions, 4 primes at odd ones
    return power(5, (n + 1) / 2) * power(4, n / 2) % M;
}

int main() {
    for (long long n : {1LL, 4LL, 50LL, 1000000000000000LL}) printf("%lld ", countGoodNumbers(n));
    printf("\n");
}

// 5 400 564908303 711414395
M = 10 ** 9 + 7


def power(x, n):                      # x^n mod M by squaring: O(log n)
    result = 1
    x %= M
    while n > 0:
        if n & 1:                     # this bit of n is 1: multiply this power of x in
            result = result * x % M
        x = x * x % M                 # x, x², x⁴, x⁸, ...: square for the next bit
        n >>= 1
    return result



def count_good_numbers(n):            # 5 even digits at even positions, 4 primes at odd ones
    return power(5, (n + 1) // 2) * power(4, n // 2) % M


print(*(count_good_numbers(n) for n in (1, 4, 50, 10 ** 15)))

# 5 400 564908303 711414395
public class Main {
    static final long M = 1000000007;

    static long power(long x, long n) {   // x^n mod M by squaring: O(log n)
        long result = 1;
        x %= M;
        while (n > 0) {
            if ((n & 1) == 1) result = result * x % M;   // this bit of n is 1: multiply this power of x in
            x = x * x % M;                             // x, x², x⁴, x⁸, ...: square for the next bit
            n >>= 1;
        }
        return result;
    }

    static long countGoodNumbers(long n) {   // 5 even digits at even positions, 4 primes at odd ones
        return power(5, (n + 1) / 2) * power(4, n / 2) % M;
    }

    public static void main(String[] args) {
        StringBuilder line = new StringBuilder();
        for (long n : new long[]{1, 4, 50, 1000000000000000L}) line.append(countGoodNumbers(n)).append(" ");
        System.out.println(line.toString().trim());
    }
}

// 5 400 564908303 711414395
const M = 1000000007n;

function power(x, n) {                 // x^n mod M by squaring: O(log n); BigInt, since x * x passes 2^53
    let result = 1n;
    x %= M;
    while (n > 0n) {
        if (n & 1n) result = result * x % M;   // this bit of n is 1: multiply this power of x in
        x = x * x % M;                         // x, x², x⁴, x⁸, ...: square for the next bit
        n >>= 1n;
    }
    return result;
}

function countGoodNumbers(n) {         // 5 even digits at even positions, 4 primes at odd ones
    return power(5n, (n + 1n) / 2n) * power(4n, n / 2n) % M;
}

console.log([1n, 4n, 50n, 1000000000000000n].map(countGoodNumbers).join(" "));

// 5 400 564908303 711414395

n = 4 gives 5² × 4² = 400, n = 50 gives 564,908,303 (the problem's own example), and n = 10¹⁵ gives 711,414,395, at once.

Square and Multiply. The loop above reads n's bits from the low end; you can also read them from the top. Start at x, the top bit; for each bit after it, square what you have, and multiply by x once more if the bit is 1. Decide each step yourself.

Square and Multiply: read n's bits from the top, square, and multiply on each 1

Pow(x, n): negative exponents

Pow(x, n) is fast power on a real number x, with n a 32-bit int that may be negative: x^−n is (1/x)^n. The trap is n = −2,147,483,648: its negation, 2,147,483,648, doesn't fit an int, so C++ and Java negate it as a 64-bit value. JavaScript's bit operators work on 32-bit ints, so it halves with arithmetic instead:

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

double myPow(double x, int n) {
    long long e = n;                       // -(-2^31) doesn't fit an int: negate as 64 bits
    if (e < 0) {
        x = 1 / x;                         // x^-e = (1/x)^e
        e = -e;
    }
    double result = 1;
    while (e > 0) {
        if (e & 1) result *= x;
        x *= x;
        e >>= 1;
    }
    return result;
}

int main() {
    printf("%.5f %.5f %.5f %.5f %.5f\n", myPow(2.0, 10), myPow(2.1, 3), myPow(2.0, -2), myPow(2.0, INT_MIN),
           myPow(1.0, INT_MIN));
}

// 1024.00000 9.26100 0.25000 0.00000 1.00000
def my_pow(x, n):
    if n < 0:
        x, n = 1 / x, -n              # x^-n = (1/x)^n; Python's ints can't overflow here
    result = 1.0
    while n > 0:
        if n & 1:
            result *= x
        x *= x
        n >>= 1
    return result


print(" ".join(f"{v:.5f}" for v in (my_pow(2.0, 10), my_pow(2.1, 3), my_pow(2.0, -2), my_pow(2.0, -2 ** 31),
                                      my_pow(1.0, -2 ** 31))))

# 1024.00000 9.26100 0.25000 0.00000 1.00000
public class Main {
    static double myPow(double x, int n) {
        long e = n;                        // -(-2^31) doesn't fit an int: negate as 64 bits
        if (e < 0) {
            x = 1 / x;                     // x^-e = (1/x)^e
            e = -e;
        }
        double result = 1;
        while (e > 0) {
            if ((e & 1) == 1) result *= x;
            x *= x;
            e >>= 1;
        }
        return result;
    }

    public static void main(String[] args) {
        System.out.println(String.format("%.5f %.5f %.5f %.5f %.5f", myPow(2.0, 10), myPow(2.1, 3), myPow(2.0, -2),
            myPow(2.0, Integer.MIN_VALUE), myPow(1.0, Integer.MIN_VALUE)));
    }
}

// 1024.00000 9.26100 0.25000 0.00000 1.00000
function myPow(x, n) {
    if (n < 0) {
        x = 1 / x;                         // x^-n = (1/x)^n
        n = -n;
    }
    let result = 1;
    while (n > 0) {                        // halve with arithmetic: JavaScript's bit operators are 32-bit
        if (n % 2 === 1) result *= x;
        x *= x;
        n = Math.floor(n / 2);
    }
    return result;
}

console.log([myPow(2.0, 10), myPow(2.1, 3), myPow(2.0, -2), myPow(2.0, -(2 ** 31)), myPow(1.0, -(2 ** 31))]
    .map((v) => v.toFixed(5)).join(" "));

// 1024.00000 9.26100 0.25000 0.00000 1.00000

2 to the power −2³¹ is so small it rounds to 0, and 1 to any power is 1. 2.1³ printed as 9.26100 at five decimals; a double can't hold 2.1 exactly, so products like this drift in the last few digits.

Fast power on matrices

Squaring works for any multiplication where (A × B) × C = A × (B × C), matrices included. The Fibonacci numbers step forward by one matrix, so F(n) is a power of it: [[1, 1], [1, 0]]^n holds F(n + 1), F(n) and F(n − 1). The same loop, with a 2 × 2 product in place of a number's:

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

const long long M = 1000000007;
typedef array<array<long long, 2>, 2> Matrix;

Matrix times(const Matrix& A, const Matrix& B) {   // 2 × 2 matrix product, mod M
    Matrix C{};
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            C[i][j] = (A[i][0] * B[0][j] + A[i][1] * B[1][j]) % M;
    return C;
}

long long fib(long long n) {           // [[1, 1], [1, 0]]^n holds F(n + 1), F(n) and F(n - 1)
    Matrix result{{{1, 0}, {0, 1}}}, step{{{1, 1}, {1, 0}}};
    while (n > 0) {                    // the same squaring as power(x, n)
        if (n & 1) result = times(result, step);
        step = times(step, step);
        n >>= 1;
    }
    return result[0][1];
}

int main() {
    printf("%lld %lld %lld\n", fib(10), fib(90), fib(1000000000000000000LL));
}

// 55 210345902 209783453
M = 10 ** 9 + 7


def times(A, B):                      # 2 × 2 matrix product, mod M
    return [[(A[i][0] * B[0][j] + A[i][1] * B[1][j]) % M for j in range(2)] for i in range(2)]


def fib(n):                           # [[1, 1], [1, 0]]^n holds F(n + 1), F(n) and F(n - 1)
    result, step = [[1, 0], [0, 1]], [[1, 1], [1, 0]]
    while n > 0:                      # the same squaring as power(x, n)
        if n & 1:
            result = times(result, step)
        step = times(step, step)
        n >>= 1
    return result[0][1]


print(fib(10), fib(90), fib(10 ** 18))

# 55 210345902 209783453
public class Main {
    static final long M = 1000000007;

    static long[][] times(long[][] A, long[][] B) {   // 2 × 2 matrix product, mod M
        long[][] C = new long[2][2];
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                C[i][j] = (A[i][0] * B[0][j] + A[i][1] * B[1][j]) % M;
        return C;
    }

    static long fib(long n) {          // [[1, 1], [1, 0]]^n holds F(n + 1), F(n) and F(n - 1)
        long[][] result = {{1, 0}, {0, 1}}, step = {{1, 1}, {1, 0}};
        while (n > 0) {                // the same squaring as power(x, n)
            if ((n & 1) == 1) result = times(result, step);
            step = times(step, step);
            n >>= 1;
        }
        return result[0][1];
    }

    public static void main(String[] args) {
        System.out.println(fib(10) + " " + fib(90) + " " + fib(1000000000000000000L));
    }
}

// 55 210345902 209783453
const M = 1000000007n;

function times(A, B) {                 // 2 × 2 matrix product, mod M (BigInt: the products pass 2^53)
    return [0, 1].map((i) => [0, 1].map((j) => (A[i][0] * B[0][j] + A[i][1] * B[1][j]) % M));
}

function fib(n) {                      // [[1, 1], [1, 0]]^n holds F(n + 1), F(n) and F(n - 1)
    let result = [[1n, 0n], [0n, 1n]], step = [[1n, 1n], [1n, 0n]];
    while (n > 0n) {                   // the same squaring as power(x, n)
        if (n & 1n) result = times(result, step);
        step = times(step, step);
        n >>= 1n;
    }
    return result[0][1];
}

console.log(`${fib(10n)} ${fib(90n)} ${fib(1000000000000000000n)}`);

// 55 210345902 209783453

F(10) = 55, F(90) mod 10⁹ + 7 = 210,345,902, and F(10¹⁸) mod 10⁹ + 7 = 209,783,453, in about 60 squarings of the matrix where the usual loop would take 10¹⁸ additions. Any recurrence that's a fixed linear step, like F(n) = F(n − 1) + F(n − 2), can be raised to a power this way; Dynamic Programming meets such recurrences.

Counting: Pascal's rule

C(n, k), "n choose k", is the number of ways to pick k of n things. Pick the k from n by looking at the last thing: either it's in (and k − 1 more come from the other n − 1) or it's out (and all k come from them). So C(n, k) = C(n − 1, k − 1) + C(n − 1, k), each number the sum of the two above it, with 1s down both edges. Pascal's Triangle builds the rows that way:

Two dashed panels. Left, rows 0 to 5 of Pascal's triangle as tiles: 1; 1 1; 1 2 1; 1 3 3 1; 1 4 6 4 1; 1 5 10 10 5 1, the edges green, the 4 and 6 of row 4 blue with arrows down to the first 10 of row 5, which is orange and ringed. Right, C(5, 2) = 10: pick 2 of 5 things; the 5th is in, so 1 more of the other 4, C(4, 1) = 4; the 5th is out, so both from the 4, C(4, 2) = 6; 4 + 6 = 10; the edges are 1: one way to take all, or none. Focal line: inner numbers: the two above, added; edges: 1. Title: each number is the two above it.
#include <bits/stdc++.h>
using namespace std;

vector<vector<int>> generate(int numRows) {   // each inner number is the two above it added
    vector<vector<int>> rows = {{1}};
    for (int r = 1; r < numRows; r++) {
        vector<int> row(r + 1, 1);
        for (int k = 1; k < r; k++) row[k] = rows[r - 1][k - 1] + rows[r - 1][k];
        rows.push_back(row);
    }
    return rows;
}

int main() {
    for (auto& row : generate(6)) {
        for (size_t k = 0; k < row.size(); k++) printf(k ? " %d" : "%d", row[k]);
        printf("\n");
    }
}

// 1
// 1 1
// 1 2 1
// 1 3 3 1
// 1 4 6 4 1
// 1 5 10 10 5 1
def generate(num_rows):               # each inner number is the two above it added
    rows = [[1]]
    for r in range(1, num_rows):
        above = rows[-1]
        rows.append([1] + [above[k - 1] + above[k] for k in range(1, r)] + [1])
    return rows


for row in generate(6):
    print(*row)

# 1
# 1 1
# 1 2 1
# 1 3 3 1
# 1 4 6 4 1
# 1 5 10 10 5 1
import java.util.*;

public class Main {
    static List<List<Integer>> generate(int numRows) {   // each inner number is the two above it added
        List<List<Integer>> rows = new ArrayList<>();
        rows.add(List.of(1));
        for (int r = 1; r < numRows; r++) {
            List<Integer> above = rows.get(r - 1), row = new ArrayList<>();
            row.add(1);
            for (int k = 1; k < r; k++) row.add(above.get(k - 1) + above.get(k));
            row.add(1);
            rows.add(row);
        }
        return rows;
    }

    public static void main(String[] args) {
        for (List<Integer> row : generate(6)) {
            StringBuilder line = new StringBuilder();
            for (int v : row) line.append(line.length() > 0 ? " " : "").append(v);
            System.out.println(line);
        }
    }
}

// 1
// 1 1
// 1 2 1
// 1 3 3 1
// 1 4 6 4 1
// 1 5 10 10 5 1
function generate(numRows) {           // each inner number is the two above it added
    const rows = [[1]];
    for (let r = 1; r < numRows; r++) {
        const above = rows[r - 1], row = [1];
        for (let k = 1; k < r; k++) row.push(above[k - 1] + above[k]);
        row.push(1);
        rows.push(row);
    }
    return rows;
}

for (const row of generate(6)) console.log(row.join(" "));

// 1
// 1 1
// 1 2 1
// 1 3 3 1
// 1 4 6 4 1
// 1 5 10 10 5 1

Row 5 is 1 5 10 10 5 1: C(5, 2) = 10 = C(4, 1) + C(4, 2) = 4 + 6. Row r takes r − 1 additions, so n rows take O(n²), which is fine for the problem's 30 rows and for any C(n, k) with n in the thousands.

What C(n, k) counts

Many counts turn into C(n, k) once one answer is written as a choice of places:

  • Ordered picks. Filling k slots in order from n things, none used twice, has n choices for the first, n − 1 for the second, and so on: P(n, k) = n × (n − 1) × … × (n − k + 1). Each set of k things comes out in k! orders, so C(n, k) = P(n, k) / k!, which is n! / (k! (n − k)!).
  • Grid paths. A walk from the top-left cell of an r × c grid to the bottom-right one, moving right or down, makes r − 1 moves down and c − 1 right. Choosing which of its r + c − 2 moves go down picks the walk: C(r + c − 2, r − 1).
  • Stars and bars. Splitting n identical things into k groups, empty ones allowed, is a row of n dots and k − 1 bars, each bar marking where a group ends. Choosing the bars' places picks the split: C(n + k − 1, k − 1).
Two dashed panels. Left, paths across a 3 × 7 grid: 21 tiles, a path of 9 blue tiles from the top-left to the bottom-right with a line along it, dark across and orange down, the rest pale; its 8 moves, in order, as tiles R, R, D, R, R, R, D, R, the two Ds in orange; choose where the 2 downs go: C(8, 2) = 28; the code walked all 28. Right, 5 things into 3 groups: three boxes holding 2, 1 and 2 blue dots; as 7 places, 5 dots and 2 bars: tiles dot, dot, bar, dot, bar, dot, dot, the bars in orange; choose where the 2 bars go: C(7, 2) = 21; the code listed all 21. Focal line: one answer, one choice of places: count with C(n, k). Title: a choice of places.

Each one listed answer by answer, then counted by its formula:

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

long long arrange(int n, int k) {          // P(n, k) = n (n - 1) ... (n - k + 1): ordered picks
    long long result = 1;
    for (int i = 0; i < k; i++) result *= n - i;
    return result;
}

long long choose(int n, int k) {           // C(n, k), exact at every step: the partial result is C(n - k + i, i)
    long long result = 1;
    for (int i = 1; i <= k; i++) result = result * (n - k + i) / i;
    return result;
}

long long paths(int r, int c) {            // walks from the top-left cell to the bottom-right, right or down
    if (r == 1 || c == 1) return 1;
    return paths(r - 1, c) + paths(r, c - 1);
}

int main() {
    int ordered = 0, splits = 0;
    for (int a = 0; a < 5; a++)            // 3 of 5 in order, none twice
        for (int b = 0; b < 5; b++)
            for (int c = 0; c < 5; c++)
                if (a != b && b != c && a != c) ordered++;
    for (int a = 0; a <= 5; a++)           // 5 into 3 groups, empty ones allowed
        for (int b = 0; b <= 5; b++)
            for (int c = 0; c <= 5; c++)
                if (a + b + c == 5) splits++;
    printf("3 of 5 in order: %d listed, P(5, 3) = %lld\n", ordered, arrange(5, 3));
    printf("paths across 3 x 7: %lld walked, C(8, 2) = %lld\n", paths(3, 7), choose(8, 2));
    printf("5 into 3 groups: %d listed, C(7, 2) = %lld\n", splits, choose(7, 2));
}

// 3 of 5 in order: 60 listed, P(5, 3) = 60
// paths across 3 x 7: 28 walked, C(8, 2) = 28
// 5 into 3 groups: 21 listed, C(7, 2) = 21
def arrange(n, k):                    # P(n, k) = n (n - 1) ... (n - k + 1): ordered picks
    result = 1
    for i in range(k):
        result *= n - i
    return result


def choose(n, k):                     # C(n, k), exact at every step: the partial result is C(n - k + i, i)
    result = 1
    for i in range(1, k + 1):
        result = result * (n - k + i) // i
    return result


def paths(r, c):                      # walks from the top-left cell to the bottom-right, right or down
    if r == 1 or c == 1:
        return 1
    return paths(r - 1, c) + paths(r, c - 1)


ordered = sum(1 for a in range(5) for b in range(5) for c in range(5)
              if a != b and b != c and a != c)           # 3 of 5 in order, none twice
splits = sum(1 for a in range(6) for b in range(6) for c in range(6)
             if a + b + c == 5)                          # 5 into 3 groups, empty ones allowed
print(f"3 of 5 in order: {ordered} listed, P(5, 3) = {arrange(5, 3)}")
print(f"paths across 3 x 7: {paths(3, 7)} walked, C(8, 2) = {choose(8, 2)}")
print(f"5 into 3 groups: {splits} listed, C(7, 2) = {choose(7, 2)}")

# 3 of 5 in order: 60 listed, P(5, 3) = 60
# paths across 3 x 7: 28 walked, C(8, 2) = 28
# 5 into 3 groups: 21 listed, C(7, 2) = 21
public class Main {
    static long arrange(int n, int k) {    // P(n, k) = n (n - 1) ... (n - k + 1): ordered picks
        long result = 1;
        for (int i = 0; i < k; i++) result *= n - i;
        return result;
    }

    static long choose(int n, int k) {     // C(n, k), exact at every step: the partial result is C(n - k + i, i)
        long result = 1;
        for (int i = 1; i <= k; i++) result = result * (n - k + i) / i;
        return result;
    }

    static long paths(int r, int c) {      // walks from the top-left cell to the bottom-right, right or down
        if (r == 1 || c == 1) return 1;
        return paths(r - 1, c) + paths(r, c - 1);
    }

    public static void main(String[] args) {
        int ordered = 0, splits = 0;
        for (int a = 0; a < 5; a++)        // 3 of 5 in order, none twice
            for (int b = 0; b < 5; b++)
                for (int c = 0; c < 5; c++)
                    if (a != b && b != c && a != c) ordered++;
        for (int a = 0; a <= 5; a++)       // 5 into 3 groups, empty ones allowed
            for (int b = 0; b <= 5; b++)
                for (int c = 0; c <= 5; c++)
                    if (a + b + c == 5) splits++;
        System.out.println("3 of 5 in order: " + ordered + " listed, P(5, 3) = " + arrange(5, 3));
        System.out.println("paths across 3 x 7: " + paths(3, 7) + " walked, C(8, 2) = " + choose(8, 2));
        System.out.println("5 into 3 groups: " + splits + " listed, C(7, 2) = " + choose(7, 2));
    }
}

// 3 of 5 in order: 60 listed, P(5, 3) = 60
// paths across 3 x 7: 28 walked, C(8, 2) = 28
// 5 into 3 groups: 21 listed, C(7, 2) = 21
function arrange(n, k) {                   // P(n, k) = n (n - 1) ... (n - k + 1): ordered picks
    let result = 1;
    for (let i = 0; i < k; i++) result *= n - i;
    return result;
}

function choose(n, k) {                    // C(n, k), exact at every step: the partial result is C(n - k + i, i)
    let result = 1;
    for (let i = 1; i <= k; i++) result = result * (n - k + i) / i;
    return result;
}

function paths(r, c) {                     // walks from the top-left cell to the bottom-right, right or down
    if (r === 1 || c === 1) return 1;
    return paths(r - 1, c) + paths(r, c - 1);
}

let ordered = 0, splits = 0;
for (let a = 0; a < 5; a++)                // 3 of 5 in order, none twice
    for (let b = 0; b < 5; b++)
        for (let c = 0; c < 5; c++)
            if (a !== b && b !== c && a !== c) ordered++;
for (let a = 0; a <= 5; a++)               // 5 into 3 groups, empty ones allowed
    for (let b = 0; b <= 5; b++)
        for (let c = 0; c <= 5; c++)
            if (a + b + c === 5) splits++;
console.log(`3 of 5 in order: ${ordered} listed, P(5, 3) = ${arrange(5, 3)}`);
console.log(`paths across 3 x 7: ${paths(3, 7)} walked, C(8, 2) = ${choose(8, 2)}`);
console.log(`5 into 3 groups: ${splits} listed, C(7, 2) = ${choose(7, 2)}`);

// 3 of 5 in order: 60 listed, P(5, 3) = 60
// paths across 3 x 7: 28 walked, C(8, 2) = 28
// 5 into 3 groups: 21 listed, C(7, 2) = 21

All three agree: 60 ordered picks of 3 from 5, 28 walks across a 3 × 7 grid, 21 splits of 5 into 3 groups. Listing takes a step per answer, so it only works on small inputs; each formula takes k steps. Unique Paths (Dynamic Programming) counts the same walks with a table, one cell at a time; for its m × n grid, C(m + n − 2, m − 1) gets there at once.

Counting mod a prime: factorials and inverses

For C(100,000, 50,000), with 30,101 digits, the formula C(n, k) = n! / (k! (n − k)!) is the way, mod a prime p. The division needs an inverse: a number b⁻¹ with b × b⁻¹ = 1 mod p. Fermat's little theorem gives it when p is prime: b^(p−1) = 1 mod p, so b^(p−2) is b's inverse, one fast power. Build the factorials mod p once, then their inverses from the top down, and every C(n, k) is three multiplications:

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

const long long M = 1000000007;
const int N = 100000;
long long fact[N + 1], invFact[N + 1];

long long power(long long x, long long n) {
    long long result = 1;
    for (x %= M; n > 0; n >>= 1, x = x * x % M)
        if (n & 1) result = result * x % M;
    return result;
}

long long comb(int n, int k) {         // C(n, k) mod M from the tables
    if (k < 0 || k > n) return 0;
    return fact[n] * invFact[k] % M * invFact[n - k] % M;
}

int main() {
    fact[0] = 1;
    for (int k = 1; k <= N; k++) fact[k] = fact[k - 1] * k % M;
    invFact[N] = power(fact[N], M - 2);      // Fermat: a^(p - 2) is a's inverse mod a prime p
    for (int k = N; k > 0; k--) invFact[k - 1] = invFact[k] * k % M;   // 1/(k-1)! = k * 1/k!
    printf("%lld %lld %lld\n", comb(5, 2), comb(1000, 500), comb(100000, 50000));
}

// 10 159835829 149033233
M = 10 ** 9 + 7
N = 100000
fact = [1] * (N + 1)
for k in range(1, N + 1):
    fact[k] = fact[k - 1] * k % M
inv_fact = [1] * (N + 1)
inv_fact[N] = pow(fact[N], M - 2, M)  # Fermat: a^(p - 2) is a's inverse mod a prime p
for k in range(N, 0, -1):
    inv_fact[k - 1] = inv_fact[k] * k % M   # 1/(k-1)! = k * 1/k!


def comb(n, k):                       # C(n, k) mod M from the tables
    if k < 0 or k > n:
        return 0
    return fact[n] * inv_fact[k] % M * inv_fact[n - k] % M


print(comb(5, 2), comb(1000, 500), comb(100000, 50000))

# 10 159835829 149033233
public class Main {
    static final long M = 1000000007;
    static final int N = 100000;
    static long[] fact = new long[N + 1], invFact = new long[N + 1];

    static long power(long x, long n) {
        long result = 1;
        for (x %= M; n > 0; n >>= 1, x = x * x % M)
            if ((n & 1) == 1) result = result * x % M;
        return result;
    }

    static long comb(int n, int k) {   // C(n, k) mod M from the tables
        if (k < 0 || k > n) return 0;
        return fact[n] * invFact[k] % M * invFact[n - k] % M;
    }

    public static void main(String[] args) {
        fact[0] = 1;
        for (int k = 1; k <= N; k++) fact[k] = fact[k - 1] * k % M;
        invFact[N] = power(fact[N], M - 2);   // Fermat: a^(p - 2) is a's inverse mod a prime p
        for (int k = N; k > 0; k--) invFact[k - 1] = invFact[k] * k % M;   // 1/(k-1)! = k * 1/k!
        System.out.println(comb(5, 2) + " " + comb(1000, 500) + " " + comb(100000, 50000));
    }
}

// 10 159835829 149033233
const M = 1000000007n;
const N = 100000;

function power(x, n) {
    let result = 1n;
    for (x %= M; n > 0n; n >>= 1n, x = x * x % M)
        if (n & 1n) result = result * x % M;
    return result;
}

const fact = new Array(N + 1).fill(1n), invFact = new Array(N + 1).fill(1n);
for (let k = 1; k <= N; k++) fact[k] = fact[k - 1] * BigInt(k) % M;
invFact[N] = power(fact[N], M - 2n);   // Fermat: a^(p - 2) is a's inverse mod a prime p
for (let k = N; k > 0; k--) invFact[k - 1] = invFact[k] * BigInt(k) % M;   // 1/(k-1)! = k * 1/k!

function comb(n, k) {                  // C(n, k) mod M from the tables
    if (k < 0 || k > n) return 0n;
    return fact[n] * invFact[k] % M * invFact[n - k] % M;
}

console.log(`${comb(5, 2)} ${comb(1000, 500)} ${comb(100000, 50000)}`);

// 10 159835829 149033233

C(5, 2) = 10 as before, C(1,000, 500) mod 10⁹ + 7 = 159,835,829, and C(100,000, 50,000) came out the same as Python's exact math.comb reduced. Only one inverse is computed, for 100,000!; the rest follow from 1/(k − 1)! = k × 1/k!. O(n) to build, O(1) per C(n, k).

When the modulus isn't prime, b has an inverse only if gcd(b, m) = 1, and the extended Euclidean algorithm finds it. It runs Euclid's algorithm (Divisors, GCD and Primes) and, on the way back, keeps x and y with a × x + b × y = gcd(a, b). When the gcd is 1, a × x = 1 mod b, so x is the inverse:

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

long long extendedGcd(long long a, long long b, long long& x, long long& y) {   // a·x + b·y = gcd(a, b)
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    long long x1, y1;
    long long g = extendedGcd(b, a % b, x1, y1);
    x = y1;                            // from b·x1 + (a mod b)·y1 = g
    y = x1 - (a / b) * y1;
    return g;
}

string inverse(long long a, long long m) {   // a's inverse mod m, when gcd(a, m) = 1
    long long x, y;
    if (extendedGcd(a, m, x, y) != 1) return "none";
    return to_string(((x % m) + m) % m);   // x can be negative
}

int main() {
    printf("%s %s %s\n", inverse(3, 11).c_str(), inverse(7, 26).c_str(), inverse(4, 26).c_str());
}

// 4 15 none
def extended_gcd(a, b):               # (g, x, y) with a·x + b·y = g = gcd(a, b)
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extended_gcd(b, a % b)
    return g, y1, x1 - (a // b) * y1  # from b·x1 + (a mod b)·y1 = g


def inverse(a, m):                    # a's inverse mod m, when gcd(a, m) = 1
    g, x, _ = extended_gcd(a, m)
    return x % m if g == 1 else "none"


print(inverse(3, 11), inverse(7, 26), inverse(4, 26))

# 4 15 none
public class Main {
    static long[] extendedGcd(long a, long b) {   // {g, x, y} with a·x + b·y = g = gcd(a, b)
        if (b == 0) return new long[]{a, 1, 0};
        long[] r = extendedGcd(b, a % b);
        return new long[]{r[0], r[2], r[1] - (a / b) * r[2]};   // from b·x1 + (a mod b)·y1 = g
    }

    static String inverse(long a, long m) {   // a's inverse mod m, when gcd(a, m) = 1
        long[] r = extendedGcd(a, m);
        if (r[0] != 1) return "none";
        return String.valueOf(((r[1] % m) + m) % m);   // x can be negative
    }

    public static void main(String[] args) {
        System.out.println(inverse(3, 11) + " " + inverse(7, 26) + " " + inverse(4, 26));
    }
}

// 4 15 none
function extendedGcd(a, b) {           // [g, x, y] with a·x + b·y = g = gcd(a, b)
    if (b === 0) return [a, 1, 0];
    const [g, x1, y1] = extendedGcd(b, a % b);
    return [g, y1, x1 - Math.floor(a / b) * y1];   // from b·x1 + (a mod b)·y1 = g
}

function inverse(a, m) {               // a's inverse mod m, when gcd(a, m) = 1
    const [g, x] = extendedGcd(a, m);
    return g === 1 ? ((x % m) + m) % m : "none";   // x can be negative
}

console.log(inverse(3, 11), inverse(7, 26), inverse(4, 26));

// 4 15 none

3 × 4 = 12 = 1 mod 11, and 7 × 15 = 105 = 4 × 26 + 1. 4 has no inverse mod 26: they share the factor 2, so 4 × x is always even and never 1 more than a multiple of 26.

Spotting it

The clues: "return it modulo 10⁹ + 7", an exponent or a count that grows with n ("how many ways", "how many strings"), n far past what a loop can reach (10¹⁵, 10¹⁸), a power with a negative exponent, or a count that's a choice of places (grid paths, splits into groups, C(n, k) for large n). Near misses:

  • Power of Two asks whether n is a power, not what a power is: one bit operation (Bit Manipulation).
  • A count with a rule between neighbours, like binary strings with no two 1s in a row, isn't a single power: that's dynamic programming, though its answer is often mod 10⁹ + 7 too.

Mistakes that cost an interview

Multiplying remainders in 32 bits, or as JavaScript numbers. 999,999,999 × 999,999,999 mod 10⁹ + 7 came out 808,348,673 as 32-bit ints and 63 as JavaScript numbers; it's 64. Use 64-bit integers, or BigInt.

Leaving a difference negative. (5 − 8) % 7 is −3 in C++, Java and JavaScript. Add m before the last mod.

Negating −2³¹ in an int. Pow(x, n) with n = -n in Java:

public class Main {
    static double myPow(double x, int n) {
        if (n < 0) {
            x = 1 / x;
            n = -n;                        // -(-2^31) is -2^31 again in an int
        }
        double result = 1;
        while (n > 0) {
            if ((n & 1) == 1) result *= x;
            x *= x;
            n >>= 1;
        }
        return result;
    }

    public static void main(String[] args) {
        System.out.println(myPow(2.0, Integer.MIN_VALUE) + ", where Math.pow gives " + Math.pow(2.0, Integer.MIN_VALUE));
    }
}

It printed 1.0: −(−2,147,483,648) is −2,147,483,648 again in an int, the loop never ran, and the result stayed 1. Negate into a 64-bit integer.

Counting Good Numbers with a loop. About 45 days at the measured rate, for n = 10¹⁵. Square.

Factorials past 64 bits. 20! is the last that fits:

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

int main() {
    long long f = 1;
    for (int k = 2; k <= 21; k++) {
        if (k == 21) printf("20! = %lld\n", f);
        f *= k;                        // 21! is past 2^63 (in C++ the overflow is undefined behaviour)
    }
    printf("21! = %lld\n", f);
}

// 20! = 2432902008176640000
// 21! = -4249290049419214848
public class Main {
    public static void main(String[] args) {
        long f = 1;
        for (int k = 2; k <= 21; k++) {
            if (k == 21) System.out.println("20! = " + f);
            f *= k;                    // 21! is past 2^63: the long wraps around
        }
        System.out.println("21! = " + f);
    }
}

// 20! = 2432902008176640000
// 21! = -4249290049419214848

21! wrapped to −4,249,290,049,419,214,848. Reduce mod p as you go, and never divide by a factorial.

Dividing remainders. C(100,000, 2) from factorials reduced mod 10⁹ + 7, divided with //, gave 0; multiplying by the inverse gave 999,949,972, which is 100,000 × 99,999 / 2 reduced. A remainder can't be divided, only multiplied by an inverse.

Check yourself

Quiz

1 / 6

Multiple choice

In C++, Java and JavaScript, what does (5 − 8) % 7 give?

Practice

In this section, easy to hard:

  • Pascal's Triangle: each inner number is the two above it added.
  • Pow(x, n): fast power on a double, a negative exponent, and −2³¹.
  • Count Good Numbers: two fast powers mod 10⁹ + 7, for n up to 10¹⁵.

Math Quiz, at the end of this section, mixes all three lessons.