Embercode
✕
Dynamic ProgrammingLesson 29 of 32

45 min read

Probability and counting DP with modulo

Why this matters

Many questions ask "how many ways?" or "how likely?". How many ways can three dice add up to 10? What is the chance that a password of length 8 has exactly two digits? How many paths of length 50 avoid a wall? The answers grow very fast. A table of counts that works for 5 dice can silently break for 40 dice, because the numbers no longer fit in the variable that stores them.

This lesson uses dice because a die is easy to picture. The same table idea works for any "count the ways to build a total one step at a time" question. You will learn three things: how to build the count table, how to turn counts into probabilities, and how to keep the numbers small enough to store by taking a remainder after each step.

The earlier lessons in this module covered the table idea for coin change and knapsack. Here the table counts ways instead of finding a best value. The modulo step, which keeps only the remainder after a division, is the new part.

Dice, totals and counting

A standard die has faces 1 to 6. When you roll one die, the total is one of six numbers, and each is equally likely. Two dice give a total from 2 to 12. The total 7 can come from six pairs: 1 and 6, 2 and 5, 3 and 4, 4 and 3, 5 and 2, 6 and 1. Two dice have 6 times 6 = 36 equally likely outcomes, so the chance of a total of 7 is 6 out of 36, which is 1 in 6.

Counting by hand works for two dice. For 40 dice it does not. We need a method that builds the answer from smaller parts. That method is dynamic programming: store the answer for a small problem in a table, then use those stored answers to fill in the next row.

Here is the row of totals for two dice. The numbers are the number of ways to get each total from 2 to 12:

total:  2  3  4  5  6  7  8  9 10 11 12
ways:   1  2  3  4  5  6  5  4  3  2  1

The total 2 has one way (1 and 1). The total 3 has two ways (1 and 2, or 2 and 1). The row rises to 6 at total 7 and falls back to 1 at total 12. That shape is the sum of two dice, and the table below shows how each cell gets its value from the row before.

Counting with a table

One die

Let ways[s] mean "the number of ways to get the total s using the dice rolled so far". With zero dice, the only total is 0, and it happens in exactly one way: roll nothing. So before any dice are rolled, ways[0] is 1 and every other entry is 0.

After one die with faces 1 to 6, the totals 1 to 6 each have one way. So ways is 0, 1, 1, 1, 1, 1, 1 for totals 0 to 6.

Adding one more die

Suppose we know the counts for i dice. What are the counts for i + 1 dice? Take a total s with i + 1 dice. The last die shows some face f from 1 to 6. The first i dice must then total s - f. So the number of ways to get s is the sum over all faces f of the number of ways to get s - f with i dice.

In a formula, with m faces per die:

new[s] = old[s - 1] + old[s - 2] + ... + old[s - m]

Only terms where s - f is 0 or more count. A total below 0 is impossible, so those terms are left out.

A worked example: two dice, total 7

We start with the one-die row, old[1] to old[6] equal to 1 and every other entry 0. To find the two-dice count for total 7, add the one-die counts at 6, 5, 4, 3, 2 and 1:

new[7] = old[6] + old[5] + old[4] + old[3] + old[2] + old[1]
       = 1 + 1 + 1 + 1 + 1 + 1
       = 6

The same method gives the total 2: only old[1] is a valid term, because old[0] would need the first die to show 0. So new[2] = old[1] = 1. That matches the row above.

The recurrence works the same way for any number of dice. Each row is built from the row before it, so a table with one row per die is enough.

The program

Here is a complete program. It takes three numbers on the command line: the number of dice, the number of faces, and the total we want. It prints the number of ways, the number of outcomes, the probability as a modular value, and the probability as a decimal number.

#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <utility>
#include <vector>
 
constexpr std::uint64_t MOD = 1'000'000'007;
 
// Number of ways to roll `target` in total with `dice` dice,
// each die showing 1 to `faces`, counted modulo MOD.
std::uint64_t count_ways(int dice, int faces, int target) {
    std::vector<std::uint64_t> prev(target + 1, 0);
    prev[0] = 1;  // zero dice: total 0 in exactly one way
    for (int i = 1; i <= dice; ++i) {
        std::vector<std::uint64_t> curr(target + 1, 0);
        for (int s = 0; s <= target; ++s) {
            for (int f = 1; f <= faces && f <= s; ++f) {
                curr[s] = (curr[s] + prev[s - f]) % MOD;
            }
        }
        prev = std::move(curr);
    }
    return prev[target];
}
 
// The same table with doubles: each face has chance 1 / faces.
double probability(int dice, int faces, int target) {
    std::vector<double> prev(target + 1, 0.0);
    prev[0] = 1.0;
    for (int i = 1; i <= dice; ++i) {
        std::vector<double> curr(target + 1, 0.0);
        for (int s = 0; s <= target; ++s) {
            for (int f = 1; f <= faces && f <= s; ++f) {
                curr[s] += prev[s - f] / faces;
            }
        }
        prev = std::move(curr);
    }
    return prev[target];
}
 
// (base ^ exp) mod MOD by repeated squaring.
std::uint64_t pow_mod(std::uint64_t base, std::uint64_t exp) {
    std::uint64_t result = 1;
    base %= MOD;
    while (exp > 0) {
        if (exp & 1) result = result * base % MOD;
        base = base * base % MOD;
        exp >>= 1;
    }
    return result;
}
 
int main(int argc, char* argv[]) {
    if (argc != 4) {
        std::cerr << "usage: dice_mod DICE FACES TOTAL\n";
        return EXIT_FAILURE;
    }
    int dice = std::atoi(argv[1]);
    int faces = std::atoi(argv[2]);
    int target = std::atoi(argv[3]);
 
    std::uint64_t ways = count_ways(dice, faces, target);
    std::uint64_t outcomes = pow_mod(faces, dice);  // faces^dice, mod MOD
    // Division mod a prime p: x / y equals x * y^(p-2) mod p (Fermat's little theorem).
    std::uint64_t prob_mod = ways * pow_mod(outcomes, MOD - 2) % MOD;
 
    std::cout << "ways mod p:      " << ways << '\n';
    std::cout << "outcomes mod p:  " << outcomes << '\n';
    std::cout << "prob mod p:      " << prob_mod << '\n';
    std::cout << "prob as double:  " << probability(dice, faces, target) << '\n';
}

Compile it and run it on the two-dice example:

$ g++ -std=c++20 -Wall -o dice_mod dice_mod.cpp
$ ./dice_mod 2 6 7
ways mod p:      6
outcomes mod p:  36
prob mod p:      166666668
prob as double:  0.166667

The count is 6, the outcomes are 36, and the probability is 6 / 36 = 1 / 6. The decimal form is 0.166667 because std::cout prints six significant digits by default. The modular form is a different kind of answer, and the next section explains it.

The table of probabilities

The probability can also be built as a table, without counting at all. Each die contributes a chance of 1 / faces to each face. So a total's probability is the sum of the probabilities of the totals before it, each divided by faces. The probability function in the program does exactly that. It uses double, which stores a decimal number with about 15 to 16 significant digits. For small problems the answer is exact enough. For very large problems the answer is a tiny number that the double cannot hold precisely, and that is where the modular form helps.

Probability is ways divided by outcomes

The probability of a total is the number of ways divided by the number of equally likely outcomes. With n dice of m faces, the number of outcomes is m to the power n. For 2 dice with 6 faces, that is 6 to the power 2, which is 36.

Integer division does not give a probability. In C++, 6 / 36 with integers is 0, because the result is rounded down. The program never does that. It takes the modular inverse instead.

What a modular inverse is

Working with a remainder is called modular arithmetic. We keep only the remainder after dividing by a number called the modulus. The modulus here is MOD = 1,000,000,007, a prime number that is often used in contests and code.

Division in modular arithmetic works through an inverse. The inverse of a number y is a number z such that y * z leaves remainder 1 when divided by MOD. For a prime MOD, Fermat's little theorem says that y to the power MOD - 2 gives that inverse. Then x / y becomes x * y^(MOD-2) with a remainder taken.

The function pow_mod computes a power and its remainder in a loop of squarings, so it takes about 30 steps for an exponent near one billion. Using the program's numbers, the probability 6 / 36 becomes 6 * pow_mod(36, MOD - 2) % MOD. That gives 166666668. This is the same as 1 / 6 in the modular world. You can check it by hand: 6 times 166666668 is 1000000008, which is MOD + 1, so its remainder is 1.

The modular answer is exact. The decimal answer is an approximation. Use the modular answer when the question asks for an integer that a judge or a checker will compare, and use the decimal answer when a person wants to read the probability.

Why the table overflows

An overflow happens when a number is too large for the type that stores it. A 32-bit int holds values up to 2,147,483,647. A 64-bit std::uint64_t holds values up to about 18 quintillion, which is 18,446,744,073,709,551,615.

The number of ways grows fast. With 6 faces, the number of outcomes for 20 dice is 6 to the power 20, which is 3,656,158,440,062,976. That fits in a 64-bit number but not in an int. With 40 dice it is 6 to the power 40, which is about 1.3 times 10 to the power 31. That does not fit in any standard integer type.

The next program is the same table without the modulo. It counts with an int for 20 dice and with a 64-bit integer for 40 dice:

#include <cstdint>
#include <iostream>
#include <vector>
 
// Same table as count_ways, but without the modulo step.
// T is the number type used for the counts.
template <typename T>
T count_no_mod(int dice, int faces, int target) {
    std::vector<T> prev(target + 1, 0);
    prev[0] = 1;
    for (int i = 1; i <= dice; ++i) {
        std::vector<T> curr(target + 1, 0);
        for (int s = 0; s <= target; ++s) {
            for (int f = 1; f <= faces && f <= s; ++f) {
                curr[s] = curr[s] + prev[s - f];
            }
        }
        prev = std::move(curr);
    }
    return prev[target];
}
 
int main() {
    // 6^20 is far above the largest int, so the int version overflows.
    std::cout << "int, 20 dice:      " << count_no_mod<int>(20, 6, 60) << '\n';
    std::cout << "uint64, 40 dice:   " << count_no_mod<std::uint64_t>(40, 6, 100) << '\n';
    std::cout << "uint64, 20 dice:   " << count_no_mod<std::uint64_t>(20, 6, 60) << '\n';
}

The output is:

int, 20 dice:      379280431
uint64, 40 dice:   6574346175042375576
uint64, 20 dice:   81987009993775

The last line is correct. The exact count for 20 dice and total 60 is 81,987,009,993,775, and the program printed that. The first line is wrong. The exact count is 81,987,009,993,775, which is far more than the int can hold, so the stored value wrapped around and became 379,280,431. The second line is also wrong. The 64-bit number wrapped around, and the printed value is the exact count taken modulo 2 to the power 64. Nothing in the output tells you that the value is wrong. The program did not stop and did not print a warning.

Some compilers can catch the int case when you ask for it. With -fsanitize=undefined, the program stops to report the first signed overflow it sees:

$ g++ -std=c++20 -Wall -fsanitize=undefined -o dice_overflow dice_overflow.cpp
$ ./dice_overflow
dice_overflow.cpp:15:35: runtime error: signed integer overflow: 2065397387 + 202789093 cannot be represented in type 'int'
int, 20 dice:      379280431
uint64, 40 dice:   6574346175042375576
uint64, 20 dice:   81987009993775

The sanitizer checks the program while it runs. The message shows the file, the line, and the two numbers that overflowed. The unsigned 64-bit lines do not trigger it, because unsigned overflow is defined to wrap around in C++. That is why the unsigned lines look plausible and are still wrong.

Keeping numbers small with a modulo

The fix is to take the remainder after each addition. If every stored value is below MOD, then the table never holds a number above MOD. Each addition is one stored value plus one more value below MOD, so the sum before the remainder is below twice MOD, which fits easily in a 64-bit integer.

The line that does this is:

curr[s] = (curr[s] + prev[s - f]) % MOD;

The remainder is not the exact count. It is the count modulo MOD. For a question that asks for the count modulo MOD, this is exactly the answer. For a probability, the program multiplies the count by the inverse of the outcomes, both taken modulo MOD, and the result is the probability modulo MOD.

Here is why the multiplication in pow_mod and in the probability line does not overflow. Each value is below MOD, which is about 1 billion. Their product is below about 1 times 10 to the power 18. A 64-bit unsigned integer holds up to about 1.8 times 10 to the power 19, so the product fits.

The 40-dice example from before, with the modulo, gives these results:

$ ./dice_mod 40 6 100
ways mod p:      675121920
outcomes mod p:  935320003
prob mod p:      200324221
prob as double:  3.32582e-05

The true count for 40 dice and total 100 is 444,579,119,430,204,638,902,177,688. Taken modulo MOD, it is 675,121,920. That matches the program's first line. The decimal probability is 3.32582e-05, which is the same as 0.0000332582 and matches the exact value to the digits shown.

Worked examples of increasing difficulty

Three dice

Three dice with 6 faces give 216 outcomes. The number of ways to total 10 is 27, so the probability is 27 / 216, which is 1 / 8 and 0.125.

$ ./dice_mod 3 6 10
ways mod p:      27
outcomes mod p:  216
prob mod p:      125000001
prob as double:  0.125

The modular probability 125000001 is the value that multiplies to 1 / 8 in the modular world. You can check it by computing 8 * 125000001 modulo MOD, which gives 1.

Chance of a total at or above a target

A common question asks for the chance of getting at least a given total. The table already holds the counts for every total up to the target. For this question we need counts above the target too, so the table must be large enough for the largest possible total, which is dice * faces.

The change is small. Make the table size dice * faces + 1, count every total, and then add up the counts from the target up to the largest total. Here is the idea in the probability table:

// Chance that the total is at least `target`.
double chance_at_least(int dice, int faces, int target) {
    int max_total = dice * faces;
    std::vector<double> prev(max_total + 1, 0.0);
    prev[0] = 1.0;
    for (int i = 1; i <= dice; ++i) {
        std::vector<double> curr(max_total + 1, 0.0);
        for (int s = 0; s <= max_total; ++s) {
            for (int f = 1; f <= faces && f <= s; ++f) {
                curr[s] += prev[s - f] / faces;
            }
        }
        prev = std::move(curr);
    }
    double chance = 0.0;
    for (int s = target; s <= max_total; ++s) chance += prev[s];
    return chance;
}

For 3 dice with 6 faces and target 10, the total chance is 0.625. Out of 216 outcomes, 135 have a total of 10 or more. You can confirm 135 by listing all 216 rolls, which is what the Python itertools module does when checking this lesson's answers.

A table that only keeps the totals it needs

The table above keeps totals up to target, and that is enough for "exact total" questions. Totals above the target can never come back down, because every face adds a positive number. So for the exact count we can stop the table at target, as the first program does. This saves memory when the target is small compared with dice * faces.

Common mistakes

Forgetting the lower bound on the faces loop

The inner loop must not look at totals below zero. If the loop checks every face without the condition f <= s, the index s - f becomes negative. A vector cannot take a negative index, so the program reads memory it does not own. This small program has that mistake:

for (int f = 1; f <= faces; ++f) {  // missing: && f <= s
    curr[s] = (curr[s] + prev[s - f]) % MOD;
}

On the Ubuntu g++ 15.2 install used to test this lesson, the program stops with this message:

/usr/include/c++/15/bits/stl_vector.h:1263: constexpr std::vector<_Tp, _Alloc>::reference std::vector<_Tp, _Alloc>::operator[](size_type) [with _Tp = long unsigned int; _Alloc = std::allocator<long unsigned int>; reference = long unsigned int&; size_type = long unsigned int]: Assertion '__n < this->size()' failed.

The important part is Assertion '__n < this->size()' failed. It means the index was outside the vector. The fix is the condition f <= s in the loop.

Dividing with integers

The probability 6 / 36 is 0 when both numbers are integers. The next program prints that value:

std::uint64_t ways = 6, outcomes = 36;
std::cout << "ways / outcomes = " << ways / outcomes << '\n';

Its output is:

ways / outcomes = 0

No error appears. The answer is simply wrong. Use the modular inverse for modular answers, or double for decimal answers. Do not divide two integer counts directly.

Multiplying two int values before the modulo

The product of two numbers near MOD is about 1 times 10 to the power 18. An int holds only about 2 times 10 to the power 9, so the product overflows long before the modulo can help. This program multiplies two int values:

int a = 1'000'000'006, b = 1'000'000'006;
std::cout << (a * b) % 1'000'000'007 << '\n';

With -fsanitize=undefined, the program reports:

mistakes.cpp:33:25: runtime error: signed integer overflow: 1000000006 * 1000000006 cannot be represented in type 'int'

The fix is to store values in std::uint64_t or long long before the multiplication, as the main program does.

Taking the modulo only at the end

If the table holds exact counts and the modulo is applied only to the final answer, the table overflows in the middle. The 40-dice example above shows the result: the printed number is wrong, and the program does not complain. Take the modulo at every addition, not at the end.

Using a modulus that is not prime

The inverse trick needs MOD to be prime. If MOD is 1,000,000,000, which is not prime, pow_mod(outcomes, MOD - 2) does not give the inverse, and the probability is wrong without any error message. Use a prime modulus, such as 1,000,000,007 or 998,244,353, for this method.

Practice

Exercise 1. How many ways can two dice with 4 faces (faces 1 to 4) total 5? Give the probability as a fraction.

Show the answer

There are 4 ways: (1, 4), (2, 3), (3, 2) and (4, 1). The number of outcomes is 4 times 4 = 16. The probability is 4 / 16, which is 1 / 4. The program ./dice_mod 2 4 5 prints ways mod p: 4 and prob as double: 0.25.

Exercise 2. What is 1 / 6 modulo 1,000,000,007? Check your answer without the program.

Show the answer

The answer is 166,666,668. The check is 6 * 166666668 = 1000000008, and 1000000008 is 1000000007 + 1, so the remainder after dividing by MOD is 1. That means 166,666,668 is the inverse of 6.

Exercise 3. How many ways can three dice with 6 faces total 3? Use the table idea to explain your count before checking it.

Show the answer

Only one way: each die shows 1, which gives 1 + 1 + 1 = 3. In the table, the total 2 after two dice has one way (1 and 1), and the only face that reaches 3 from there is 1. The program ./dice_mod 3 6 3 prints ways mod p: 1.

Exercise 4. Why is the product a * b % MOD safe in a 64-bit integer when a and b are both below MOD?

Show the answer

Both values are below about 1 billion, so their product is below about 1 times 10 to the power 18. A 64-bit unsigned integer holds up to about 1.8 times 10 to the power 19, so the product fits before the remainder is taken.

Exercise 5. Change count_ways so it returns the number of ways for every total from 0 to dice * faces, and print the list for 2 dice with 6 faces.

Show the answer

Make the table size dice * faces + 1 instead of target + 1, and run the loop for every s in that range. The printed list for 2 dice with 6 faces is 0, 0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1 (totals 0 to 12). The values match the row shown in the lesson, with 0 for totals 0 and 1, because the first real total is 2.

Summary

  • Count the ways to reach each total with a table. Each new die adds the counts of the faces totals before it.
  • The probability is the number of ways divided by faces to the power dice. Use a modular inverse for an exact modular answer, and double for a decimal answer.
  • Integer division rounds down, so 6 / 36 is 0. Do not use it for probabilities.
  • Overflow happens silently for unsigned types and is undefined for int. Take the modulo at every addition, and keep every stored value below the modulus.
  • The modulus used with a modular inverse must be prime.
Check yourself: Dice probability (7 questions) →Practice this topic: coding problems for Dynamic Programming →

Stuck? Ask the tutor

Ask about anything in this lesson. Answers are short and use simple words.

Sign in to ask the tutor.

← Reconstructing the answer by storing choices