Digit & Probability DP

Most dynamic programs in this section index their subproblems by something you can point at: a prefix of an array, a remaining capacity, a subset encoded as a mask. Two classic patterns are harder to see because their state lives in the structure of a single object — the decimal digits of a number, or the probabilities on the edges of a process — rather than in an input array. Both still rest on the same two requirements from the principles lesson: overlapping subproblems and a substructure that lets us recombine them. What changes is what the state must remember and what value it accumulates.

Digit DP: counting numbers, not enumerating them

Many problems ask: how many integers in $[L,R]$ satisfy some property of their digits? — no 4 anywhere, digit sum divisible by $7$, no two equal adjacent digits, at most three distinct digits. The ranges are astronomical ($R$ up to $10^{18}$), so we cannot loop $x=L\dots R$. But the property depends only on the digits, and there are at most $19$ of them, so a DP over digit positions runs in time proportional to the number of digits, not their value.

The first move is a standard reduction. Let $f(N)$ count the valid integers in $[0,N]$. Then the count in $[L,R]$ is $f(R)-f(L-1)$, so it suffices to solve the prefix problem $f(N)$, counting valid numbers from $0$ up to a single upper bound $N$.

Why a plain digit-by-digit count is not enough

Fix $N$ and write it as a digit string $N=d_{m-1}{d}_{m-2}\cdots d_0$. We build a candidate number one digit at a time, most-significant first, and at each position choose a digit. The complication is the upper bound. As long as every digit we have placed so far equals the corresponding digit of $N$, the number is still pinned to the boundary, so the next digit may range only up to $d_i$ — go higher and we exceed $N$. But the moment we place a digit strictly below$d_i$, the prefix drops under $N$, and every remaining position is free to use any digit $0$ through $9$ without risk of overflow.

That single bit of history — is the prefix still equal to $N$'s prefix? — is exactly the extra state digit DP needs. Call it the tight flag.

The figure shows the crux: there is exactly one tight path at any depth — the one that has matched $N$ digit-for-digit — and it is the only place where the next digit is capped. Everything hanging off a below-cap choice is free, and free subtrees with the same remaining length and the same carried information are identical, which is precisely the overlapping-subproblem structure that makes this a DP and not a $10^m$ enumeration.¹

The state and recurrence

Beyond position and the tight flag, the state must carry whatever the property needs to be checkable incrementally. For digit sum $\equiv 0(\mathrm{mod}7)$ that is the running sum modulo $7$; for no 4 it is nothing extra (we just forbid the digit 4); for no equal adjacent digits it is the previous digit. Write that problem-specific carry as $s$.

At position $i$ the choosable digits are $d\in \{0,\dots ,\text{cap}\}$, where $\text{cap}=d_i$ when tight and $\text{cap}=9$ when free. Placing $d$ advances to the next position, updates the carry to $s^{\prime }=\text{step}(s,d)$, and keeps the prefix tight only if it was tight and $d$ hit the cap:

The base case is reaching past the last position, $i=-1$: return $1$ if the accumulated carry $s$ marks a valid number (e.g. $s\equiv 0$ for the divisibility property) and $0$ otherwise. The answer is $\text{cnt}(m-1,s_0,\text{true})$, starting tight at the most-significant digit.

This is why the carry $s$ must be chosen to be exactly the information the validity test consumes — no more (or the state space blows up) and no less (or the substructure breaks). Memoizing on $(i,s,\text{tight})$ collapses the exponential tree to a table.²

To ground the recurrence, count the integers in $[0,325]$ that contain no digit 4. Here the carry $s$ is empty — validity depends only on forbidding the digit 4 — so a free state is keyed by position alone. Let $g(i)$ be the number of valid ways to fill $i$ remaining free positions: each place picks any of the $9$ non-4 digits, so $g(i)=9^i$, giving $g(0)=1$, $g(1)=9$, $g(2)=81$. Now walk the tight path of $N=325$, and at each position sum the free completions of the below-cap digits:

• Hundreds (cap $3$): a leading digit in $\{0,1,2\}$ is below the cap and none is 4, so each opens $g(2)=81$ free completions: $3\times 81=243$. Digit 3 keeps the prefix tight and continues down.
• Tens (cap $2$, prefix 3): below-cap digits $\{0,1\}$ are valid, each with $g(1)=9$ completions: $2\times 9=18$. Digit 2 stays tight.
• Units (cap $5$, prefix 32): below-cap digits $\{0,1,2,3\}$ are valid (excluding 4, which is below the cap but forbidden), contributing $4\times g(0)=4$. Digit 5 completes the tight path at $N=325$ itself, which contains no 4, so it counts as $1$.

Summing the tight path: $243+18+4+1=266$, matching a brute-force count over all $326$ integers. The DP touched three positions instead of enumerating them.

Note the memo caches only the free states: tight states lie on the single boundary path, are visited once, and depend on $N$'s specific digits, so caching them would be both useless and unsound across different bounds. With $m\le 19$ positions, a carry $s$ from a small set $S$, and $10$ digit choices, the running time is $O(m\cdot |S|\cdot 10)$ — for no 4 that is a few hundred operations to count over $[0,10^{18}]$. Count Numbers with Unique Digits and Numbers With Repeated Digits are this template with $s$ tracking a $10$-bit mask of used digits; the latter is cleanest as $R-(\text{count with all-distinct digits})$.³

from typing import Callable, Generic, Hashable, TypeVar

# the problem-specific summary carried alongside (position, tight).
Carry = TypeVar("Carry", bound=Hashable)

class DigitProperty(Generic[Carry]):
  """
    A digit property expressed as the three pieces a digit DP needs:\n
    the initial carry, how placing a digit updates it, and whether a fully\n
    built number with a given carry is accepted. Choose `Carry` to be exactly\n
    the information the validity test consumes — no more (the table blows up)\n
    and no less (the substructure breaks).\n
  """

  def __init__(
    self,
    initial: Carry,
    step: Callable[[Carry, int, int], Carry],
    accept: Callable[[Carry], bool],
  ) -> None:
    """
      `initial` is the carry before any digit is placed. `step(carry, digit,\n
      position)` returns the carry after appending `digit` at `position` (the\n
      number of less-significant places still to fill). `accept(carry)` decides\n
      whether a completed number with that carry counts.\n
    """
    self.initial: Carry = initial
    self.step: Callable[[Carry, int, int], Carry] = step
    self.accept: Callable[[Carry], bool] = accept

def count_up_to(bound: int, property: DigitProperty[Carry]) -> int:
  """
    Count integers in [0, bound] whose digits satisfy `property`.\n
    This is the prefix count f(N) of the lesson's DigitCount(N).\n
  """
  if bound < 0:
    return 0

  # decompose the bound into its decimal digits, most-significant first.
  digits: list[int] = [int(character) for character in str(bound)]
  length: int = len(digits)

  # memo holds only free states, keyed by (position, carry); tight states sit
  # on the single boundary path, are visited once, and depend on the bound's
  # own digits, so caching them would be unsound across different bounds.
  memo: dict[tuple[int, Carry], int] = {}

  def recurse(position: int, carry: Carry, tight: bool) -> int:
    # a completed number counts iff its final carry is accepted.
    if position == length:
      return 1 if property.accept(carry) else 0

    # free states repeat across the tree, so serve them from the memo.
    if not tight and (position, carry) in memo:
      return memo[(position, carry)]

    # tight rows are capped at the bound's digit; free rows span 0..9.
    cap: int = digits[position] if tight else 9
    remaining: int = length - position - 1

    # try every digit, staying tight only while we match the cap.
    total: int = 0
    for digit in range(cap + 1):
      next_carry: Carry = property.step(carry, digit, remaining)
      total += recurse(position + 1, next_carry, tight and digit == cap)

    # cache free states only; tight states live on the unique boundary path.
    if not tight:
      memo[(position, carry)] = total
    return total

  return recurse(0, property.initial, True)

def count_in_range(low: int, high: int, property: DigitProperty[Carry]) -> int:
  """
    Count integers in [low, high] satisfying `property`, via f(high) -\n
    f(low - 1). Both bounds are inclusive.\n
  """
  if low > high:
    return 0
  return count_up_to(high, property) - count_up_to(low - 1, property)

def without_forbidden_digit(forbidden: int) -> DigitProperty[bool]:
  """
    Property: the digit `forbidden` (e.g. 4) appears nowhere. The carry is a\n
    single bool — True while the number is still clean.\n
  """
  return DigitProperty(
    initial=True,
    step=lambda clean, digit, _position: clean and digit != forbidden,
    accept=lambda clean: clean,
  )

def digit_sum_divisible_by(modulus: int) -> DigitProperty[int]:
  """
    Property: the sum of digits is divisible by `modulus`. The carry is the\n
    running digit sum reduced modulo `modulus`.\n
  """
  return DigitProperty(
    initial=0,
    step=lambda residue, digit, _position: (residue + digit) % modulus,
    accept=lambda residue: residue == 0,
  )

Probability & Expectation DP: averaging instead of optimizing

The second pattern keeps the DP skeleton but changes the operator. An optimization DP combines child values with $\min$ or $\max$; a probability DP combines them with a probability-weighted sum. The value stored in a state is no longer the best you can do from here but the expected value of the random process started from here. The key fact is linearity of expectation: the expected value of a state is the average of its successors' expected values, each weighted by the transition probability — and this holds whether or not the transitions are independent.

The classic instance is the expected number of steps to absorption in a random walk. Consider a process that moves between states by chance and eventually reaches a terminal (absorbing) state. We want the expected number of steps from each starting state.

Each step costs $1$, and after that step we are in a successor $u$ chosen with probability $p(v,u)$, from which the remaining expected cost is $E[u]$. Linearity of expectation turns this into one equation per state:

The $1$ is the step just taken; the sum is the expected remaining cost, averaged over where that step lands.⁴ For the figure, $E[A]=1+\frac{1}{2}E[B]+\frac{1}{2}\cdot 0$ and $E[B]=1+\frac{1}{2}E[A]+\frac{1}{2}\cdot 0$; solving the two-by-two system gives $E[A]=8/3$ and $E[B]=7/3$.

When is this a DP, and when a linear system?

The recurrence above is a genuine DP — solvable by memoization or a bottom-up sweep — exactly when the dependency graph on states is acyclic, so that each $E[v]$ depends only on states already computed. That is the common, easy case: games that strictly advance (a counter that only grows, a round number that only increases, soup that only gets consumed), where you evaluate states in reverse topological order and read off the answer.

The same $E[v]=1+\sum p(v,u)E[u]$ shape computes expected values other than step counts — replace the constant $1$ by the per-step reward, or drop it and let an accepting state contribute its payoff — and hitting probabilities, where $P[v]={\sum }_{u}p(v,u)P[u]$ with boundary states pinned to $0$ or $1$. The acyclic, DP-friendly versions are the staple of competitive problems.

A worked acyclic example: expected die rolls to a target

Roll a fair $6$-sided die repeatedly, summing the pips, and stop the instant the running total is $\ge T$. What is the expected number of rolls? Let $E[t]$ be the expected additional rolls when the current total is $t<T$ (and $E[t]=0$ for $t\ge T$). One roll lands on $1,\dots ,6$ each with probability $\frac{1}{6}$:

Because every transition strictly increases $t$, the dependency graph is acyclic and we fill $E$ from $t=T-1$ downward — a textbook one-dimensional DP. Soup Servings is the same shape in two dimensions (the two soup volumes only decrease, so the state space is acyclic and finite), and New 21 Game is a probability-DP over a sliding window of point totals, made $O(n)$ by maintaining a running sum of the window of $E$-values rather than re-summing each transition.

The reverse-topological sweep is what makes this $O(T)$ and not a linear-system solve: each $E[t]$ reads only larger totals, already finalized. The window sum trick (subtract the leaving term, add the entering one) drops the inner loop and yields $O(T)$ even when the die has many faces — the same optimization that turns New 21 Game from $O(nk)$ into $O(n)$.

def expected_die_rolls(target: int, faces: int = 6) -> float:
  """
    Expected number of fair `faces`-sided rolls until the running total\n
    reaches at least `target`, starting from total 0. A non-positive target\n
    needs no rolls.\n
  """
  if faces < 1:
    raise ValueError("a die needs at least one face")
  if target <= 0:
    return 0.0

  # expected[total] = expected additional rolls from this total; totals >=
  # target are absorbing with value 0 and are read as such past the array end.
  expected: list[float] = [0.0 for _ in range(target + 1)]

  # window_sum tracks E[total+1] + ... + E[total+faces] so the inner loop over
  # outcomes collapses to one subtract-and-add per state.
  window_sum: float = 0.0
  for total in range(target - 1, -1, -1):
    # slide the window to cover total+1..total+faces: add the new low end,
    # drop the outcome that just fell off the high end.
    window_sum += expected[total + 1] if total + 1 <= target else 0.0
    leaving: int = total + 1 + faces
    if leaving <= target:
      window_sum -= expected[leaving]

    # one roll plus the average of where it lands.
    expected[total] = 1.0 + window_sum / faces

  return expected[0]

from typing import Generic, Hashable, TypeVar

# the label identifying a state of the chain.
Label = TypeVar("Label", bound=Hashable)

class State(Generic[Label]):
  """
    One state of the chain: its label and its outgoing transitions as a map\n
    from successor label to probability. An empty transition map marks an\n
    absorbing state — the process stops there.\n
  """

  def __init__(self, label: Label) -> None:
    self.label: Label = label
    self.transitions: dict[Label, float] = {}

  def add_transition(self, target: Label, probability: float) -> None:
    """
      Record a step to `target` taken with `probability`.\n
    """
    self.transitions[target] = self.transitions.get(target, 0.0) + probability

  @property
  def is_absorbing(self) -> bool:
    """
      Whether this state has no outgoing transitions.\n
    """
    return len(self.transitions) == 0

  def __repr__(self) -> str:
    return f"State({self.label!r})"

class MarkovChain(Generic[Label]):
  """
    A discrete-time Markov chain over labelled states.\n
    States with no outgoing transitions are absorbing.\n
  """

  def __init__(self) -> None:
    self._states: dict[Label, State[Label]] = {}

  def state(self, label: Label) -> State[Label]:
    """
      The state for `label`, creating it on first reference.\n
    """
    if label not in self._states:
      self._states[label] = State(label)
    return self._states[label]

  def add_transition(
    self, source: Label, target: Label, probability: float
  ) -> None:
    """
      Add a `source -> target` step taken with `probability`. Both states are\n
      created as needed.\n
    """
    self.state(target)
    self.state(source).add_transition(target, probability)

  @property
  def labels(self) -> list[Label]:
    """
      Every state label, in insertion order.\n
    """
    return list(self._states)

  def expected_steps_to_absorption(self) -> dict[Label, float]:
    """
      For each state, the expected number of steps to reach any absorbing\n
      state, solving E = 1 + P E over the transient states by Gaussian\n
      elimination. Handles cyclic chains where a reverse sweep cannot.\n
    """
    # index the transient states; absorbing states are pinned to 0 separately.
    transient: list[Label] = [
      label for label in self._states if not self._states[label].is_absorbing
    ]
    position: dict[Label, int] = {
      label: index for index, label in enumerate(transient)
    }
    size: int = len(transient)

    # build (I - Q) E = 1 over transient states; absorbing successors carry
    # E = 0 and so drop out of each row.
    matrix: list[list[float]] = [[0.0 for _ in range(size)] for _ in range(size)]
    constants: list[float] = [1.0 for _ in range(size)]
    for label in transient:
      row: int = position[label]
      matrix[row][row] = 1.0
      for target, probability in self._states[label].transitions.items():
        if target in position:
          matrix[row][position[target]] -= probability

    # solve, then scatter the answers back; absorbing states keep their 0.
    solution: list[float] = _solve_linear_system(matrix, constants)
    expected: dict[Label, float] = {label: 0.0 for label in self._states}
    for label in transient:
      expected[label] = solution[position[label]]
    return expected

  def hitting_probability(self, target_set: set[Label]) -> dict[Label, float]:
    """
      For each state, the probability of ever reaching a state in\n
      `target_set`. Targets are pinned to 1; other absorbing states to 0;\n
      transient states solve P[v] = sum_u p(v, u) P[u] by Gaussian\n
      elimination.\n
    """
    # transient states exclude targets (pinned to 1) and other absorbers (0).
    transient: list[Label] = [
      label
      for label in self._states
      if not self._states[label].is_absorbing and label not in target_set
    ]
    position: dict[Label, int] = {
      label: index for index, label in enumerate(transient)
    }
    size: int = len(transient)

    # build P[v] = sum_u p(v, u) P[u]: a step into a target adds its known 1
    # to the constant; a step into another transient stays on the left.
    matrix: list[list[float]] = [[0.0 for _ in range(size)] for _ in range(size)]
    constants: list[float] = [0.0 for _ in range(size)]
    for label in transient:
      row: int = position[label]
      matrix[row][row] = 1.0
      for target, probability in self._states[label].transitions.items():
        if target in target_set:
          constants[row] += probability
        elif target in position:
          matrix[row][position[target]] -= probability

    # pin the boundary: targets hit with probability 1, other absorbers 0.
    solution: list[float] = _solve_linear_system(matrix, constants)
    probability_of_hit: dict[Label, float] = {}
    for label in self._states:
      if label in target_set:
        probability_of_hit[label] = 1.0
      elif self._states[label].is_absorbing:
        probability_of_hit[label] = 0.0

    # fill in the solved transient probabilities.
    for label in transient:
      probability_of_hit[label] = solution[position[label]]
    return probability_of_hit

def _solve_linear_system(
  matrix: list[list[float]], constants: list[float]
) -> list[float]:
  """
    Solve `matrix @ solution = constants` by Gaussian elimination with partial\n
    pivoting. Returns the solution vector; an empty system yields an empty\n
    vector.\n
  """
  size: int = len(constants)
  if size == 0:
    return []

  # work on an augmented copy so the caller's data is untouched.
  augmented: list[list[float]] = [
    matrix[row][:] + [constants[row]] for row in range(size)
  ]

  for column in range(size):
    # partial pivot: bring the largest-magnitude entry into the pivot row.
    pivot_row: int = max(
      range(column, size),
      key=lambda candidate: abs(augmented[candidate][column]),
    )
    if abs(augmented[pivot_row][column]) < 1e-12:
      raise ValueError("singular system: no unique solution")
    augmented[column], augmented[pivot_row] = (
      augmented[pivot_row],
      augmented[column],
    )

    # eliminate this column from every other row.
    pivot: float = augmented[column][column]
    for target_row in range(size):
      if target_row == column:
        continue
      factor: float = augmented[target_row][column] / pivot
      for term in range(column, size + 1):
        augmented[target_row][term] -= factor * augmented[column][term]

  # the matrix is now diagonal: read each unknown straight off its row.
  return [augmented[row][size] / augmented[row][row] for row in range(size)]

What these two patterns really are

Digit DP is a special case of a much older idea: counting the accepting paths of a finite automaton over a bounded input. The digit-by-digit walk traces a run of a deterministic automaton whose states are the carries $s$, and the count valid numbers $\le N$ question is the automaton-intersection count between that property-automaton and the less-than-or-equal-to-$N$ automaton — the tight flag is the state of the second automaton. Once framed this way, digit DP generalizes far past base $10$: the same technique counts binary strings avoiding a pattern (a transfer-matrix / regular-language count), lattice points under algebraic constraints, and, in the theory of automatic sequences (Allouche and Shallit, Automatic Sequences, 2003), whole families of number-theoretic counting functions. For a fixed carry-automaton the count is a transfer-matrix product, so $f(N)$ can even be computed by matrix exponentiation in $O(|S{|}^3\log N)$ when the number of positions is huge.⁵

Expectation DP is standard Markov chain theory. The system $E=1+PE$ is the standard expected hitting time to an absorbing set, and the matrix that makes it solvable is the fundamental matrix$N=(I-Q{)}^{-1}$ of the transient part $Q$, whose row sums give the expected steps to absorption (Kemeny and Snell, Finite Markov Chains, 1960). The acyclic DP case is precisely when $Q$ is nilpotent (can be permuted to strictly triangular form), so the inverse is a finite sum and no linear solve is needed — the reverse topological sweep is Gaussian elimination that happens to require no back-substitution. This is the same object behind absorbing random walks on graphs, the gambler's ruin problem, and the stationary analysis that underlies Google's original PageRank (Page, Brin, Motwani, Winograd, 1999), where the relevant quantity is the stationary distribution of a chain rather than a hitting time, but the linear-algebra machinery is identical. The dividing line this lesson draws — acyclic means DP, cyclic means solve the system — is the boundary between a nilpotent and a general transient matrix.

Takeaways

• Digit DP counts integers in $[L,R]$ with a digit property by solving the prefix problem $f(N)$ and using $f(R)-f(L-1)$; it walks the $\le 19$ decimal positions of the bound rather than the $10^{18}$ values.
• The essential extra state is the tight flag: while the built prefix equals $N$'s prefix the next digit is capped at $d_i$; the first below-cap digit makes it free, after which all later digits range over $0$–$9$. State is $(\text{position},\text{carry }s,\text{tight})$, and only free states are memoized.
• Choose the carry $s$ to be exactly what the validity test reads (running sum mod $k$, previous digit, used-digit mask) — no more, no less — so the substructure holds and the table stays small.
• Expectation DP stores $E[v]$ = expected value from state $v$ and combines children with a probability-weighted sum, $E[v]=1+{\sum }_{u}p(v,u)E[u]$, by linearity of expectation; absorbing states have $E=0$.
• It is a plain DP when the state graph is acyclic (quantities that only advance) — evaluate in reverse topological order — and a linear system when transitions form cycles, solved by Gaussian elimination instead.