Routes that visit every node once

def count_hamiltonian(n, edges, source, destination):
    out = [[] for _ in range(n)]
    for u, v in edges:
        out[u].append(v)
    full = (1 << n) - 1
    # dp[mask][v]: paths from source covering exactly mask, ending at v
    dp = [[0] * n for _ in range(1 << n)]
    dp[1 << source][source] = 1
    for mask in range(1 << n):
        if not mask >> source & 1:
            continue                    # every path contains the source
        for v in range(n):
            ways = dp[mask][v]
            if ways == 0 or v == destination and mask != full:
                continue                # the destination must be the last node
            for w in out[v]:
                if not mask >> w & 1:
                    dp[mask | 1 << w][w] += ways
    return dp[full][destination]

#include <vector>

long long solve(int n, const std::vector<std::vector<int>>& edges, int source, int destination) {
    std::vector<std::vector<int>> out(n);
    for (const auto& e : edges) out[e[0]].push_back(e[1]);
    const int full = (1 << n) - 1;
    // dp[mask][v]: paths from source covering exactly mask, ending at v
    std::vector<std::vector<long long>> dp(1 << n, std::vector<long long>(n, 0));
    dp[1 << source][source] = 1;
    for (int mask = 0; mask <= full; ++mask) {
        if (!(mask >> source & 1)) continue;
        for (int v = 0; v < n; ++v) {
            long long ways = dp[mask][v];
            if (ways == 0 || (v == destination && mask != full)) continue;   // destination comes last
            for (int w : out[v])
                if (!(mask >> w & 1)) dp[mask | 1 << w][w] += ways;
        }
    }
    return dp[full][destination];
}

import java.util.*;

class Solution {
    public long countHamiltonian(int n, int[][] edges, int source, int destination) {
        List<List<Integer>> out = new ArrayList<>();
        for (int i = 0; i < n; i++) out.add(new ArrayList<>());
        for (int[] e : edges) out.get(e[0]).add(e[1]);
        int full = (1 << n) - 1;
        long[][] dp = new long[1 << n][n];             // paths covering mask, ending at v
        dp[1 << source][source] = 1;
        for (int mask = 0; mask <= full; mask++) {
            if ((mask >> source & 1) == 0) continue;
            for (int v = 0; v < n; v++) {
                long ways = dp[mask][v];
                if (ways == 0 || (v == destination && mask != full)) continue;   // destination comes last
                for (int w : out.get(v))
                    if ((mask >> w & 1) == 0) dp[mask | 1 << w][w] += ways;
            }
        }
        return dp[full][destination];
    }
}

function countHamiltonian(n, edges, source, destination) {
  const out = Array.from({ length: n }, () => []);
  for (const [u, v] of edges) out[u].push(v);
  const full = (1 << n) - 1;
  // dp[mask * n + v]: paths from source covering exactly mask, ending at v
  const dp = new Array((1 << n) * n).fill(0);
  dp[(1 << source) * n + source] = 1;
  for (let mask = 0; mask <= full; mask++) {
    if (!(mask >> source & 1)) continue;
    for (let v = 0; v < n; v++) {
      const ways = dp[mask * n + v];
      if (ways === 0 || (v === destination && mask !== full)) continue;   // destination comes last
      for (const w of out[v]) if (!(mask >> w & 1)) dp[(mask | 1 << w) * n + w] += ways;
    }
  }
  return dp[full * n + destination];
}

fn solve(n: i64, edges: &[[i64; 2]], source: i64, destination: i64) -> i64 {
    let (n, s, d) = (n as usize, source as usize, destination as usize);
    let mut out = vec![Vec::new(); n];
    for e in edges {
        out[e[0] as usize].push(e[1] as usize);
    }
    let full = (1usize << n) - 1;
    // dp[mask][v]: paths from source covering exactly mask, ending at v
    let mut dp = vec![vec![0i64; n]; 1 << n];
    dp[1 << s][s] = 1;
    for mask in 0..=full {
        if mask >> s & 1 == 0 {
            continue;
        }
        for v in 0..n {
            let ways = dp[mask][v];
            if ways == 0 || (v == d && mask != full) {
                continue;                                  // destination comes last
            }
            for &w in &out[v] {
                if mask >> w & 1 == 0 {
                    dp[mask | 1 << w][w] += ways;
                }
            }
        }
    }
    dp[full][d]
}

import Data.Bits (setBit, testBit, (.|.), shiftL)
import qualified Data.Map.Strict as M

solve :: Int -> [[Int]] -> Int -> Int -> Int
solve n edges source dest = M.findWithDefault 0 (full, dest) table
  where
    full = (1 `shiftL` n) - 1 :: Int
    out = M.fromListWith (++) [(u, [v]) | [u, v] <- edges]
    -- table: (mask, v) -> number of paths from source covering exactly mask and ending at v
    table = foldl step (M.singleton (setBit 0 source, source) 1) [0 .. full]
    step t mask = foldl (extend mask) t [0 .. n - 1]
    extend mask t v = case M.lookup (mask, v) t of
      Just ways | not (v == dest && mask /= full) ->
        foldl (\acc w -> M.insertWith (+) (mask .|. setBit 0 w, w) ways acc) t
              [w | w <- M.findWithDefault [] v out, not (testBit mask w)]
      _ -> t