The short answer is, no one knows. Initially people suspected this had to do with cache misses. But as my example above shows, I get the slowdown using a tiny amount of memory. Also, even after countless attempts, no one has been able to reproduce this slowdown locally. So the issue is likely something on Codeforces' side.

There are lots of alternative ways to trigger the slowdown bug:

It doesn't have to be a vector of vectors 249329395 (credit to kostia244).

You can change the size slightly, and still get the slowdown bug. From testing 39999 and 40000, as well as 7,8,9,10 can all trigger the slowdown bug. See for example 249330544. But 40000 and 7 is the only combination I've found that consistently trigger it.

Declaring the vector inside main also causes the slowdown bug 249330046.

You can still trigger the slowdown with vector<vector<char>>, but you need to increase $$$7$$$ to $$$4 \cdot 7$$$ 249331394.

Submissions using scanf instead of cin seems to be less affected by the slowdown bug. But even scanf submissions typically run around a factor of $$$2$$$ slower. For example, Radewoosh's solution to problem 1936-D - Bitwise Paradox goes from taking 2168 ms to 4336 ms.

Over the last couple of years, there has been many blogs about people getting TLE on Codeforces for no appearent reason. Here are a couple examples

Strange TLE by cin using GNU C++20 (64), and some weird subtle changes to (surprisingly) fix it by -14

Bizarre slowdown when using cin on a 64-bit compiler by defnotmee

Inconsistent times and TLE with C++17 (64bit) vs. C++17 by OleschY

Also if you dig around in these blogs, you will find comments that points to even more examples. What all of these weird TLEs seems to have in common is that

1. Submitted using 64 bit C++ on Codeforces.

2. Uses std::cin.

3. Unable to be reproduced locally.

It is for this reason that I suspect that these unexplainable TLEs all stem from the same fundamental issue on Codeforces' side.

Remark: I'm no longer to reproduce the issue in Inconsistent times and TLE with C++17 (64bit) vs. C++17. So that blog might be about a completely different issue that has now been fixed.

#define log std::__lg
#define ctz __builtin_ctz 
 
// Rooted tree
struct arborescence {
    std::vector<std::vector<int>> children;
    int root;
};
 
arborescence shallowest_decomposition_tree(std::vector<std::vector<int>> &graph, int root=0) {
    int n = (int)graph.size();
    std::vector<std::vector<int>> decomposition_tree(n), stacks(log(n) + 1);
    
    auto extract_chain = [&](int labels, int u) {
        while (labels) {
            int label = log(labels);
            labels ^= 1 << label;
            int v = stacks[label].back(); stacks[label].pop_back();
            decomposition_tree[u].push_back(v);
            u = v;
        }
    };
    std::vector<int> forbidden(n, -1);
    auto dfs = [&](int u, int p, auto&& self) -> void {
        int forbidden_once = 0, forbidden_twice = 0;
        for (auto v : graph[u]) {
            if (v != p) {
                self(v, u, self);
                forbidden_twice |= forbidden_once & (forbidden[v] + 1);
                forbidden_once |= forbidden[v] + 1;
            }
        }
        forbidden[u] = forbidden_once | ((1 << log(2*forbidden_twice  + 1)) - 1);
        int label_u = ctz(forbidden[u] + 1);
        stacks[label_u].push_back(u);
        for (int i = (int)graph[u].size() - 1; i >= 0; --i) {
            int v = graph[u][i];
            extract_chain((forbidden[v] + 1) & ((1 << label_u) - 1), u);
        }
    };
    dfs(root, -1, dfs);
    int max_label = log(forbidden[root] + 1);
    int decomposition_root = stacks[max_label].back(); stacks[max_label].pop_back();
    extract_chain((forbidden[root] + 1) & ((1 << max_label) - 1), decomposition_root);
    return {std::move(decomposition_tree), decomposition_root};
}

def ctz(x): # Count trailing zeroes
    return (x & -x).bit_length() - 1
 
def shallowest_decomposition_tree(graph, root=0):
    n = len(graph)
    forbidden = [-1] * n
    decomposition_tree = [[] for _ in range(n)]
    stacks = [[] for _ in range(n.bit_length())]
    def extract_chain(labels, u):
        while labels:
            label = labels.bit_length() - 1
            labels ^= 2**label
            v = stacks[label].pop()
            decomposition_tree[u].append(v)
            u = v
    def dfs(u, p):
        forbidden_once = forbidden_twice = 0
        for v in graph[u]:
            if v != p: 
              dfs(v, u)
              forbidden_twice  |= forbidden_once & (forbidden[v] + 1)
              forbidden_once  |= forbidden[v] + 1
        forbidden[u] = forbidden_once | (2**forbidden_twice.bit_length() - 1)    
        label_u = ctz(forbidden[u] + 1)
        stacks[label_u].append(u)
        for v in reversed(graph[u]): 
            extract_chain((forbidden[v] + 1) & (2**label_u - 1), u)
    dfs(root, -1)
    max_label = (forbidden[root] + 1).bit_length() - 1
    decomposition_root = stacks[max_label].pop()
    extract_chain((forbidden[root] + 1) & (2**max_label - 1), decomposition_root)
    return decomposition_tree, decomposition_root

Since doing deep recursion in Python is generally a really bad idea. I've also implemented a slightly less nice looking version without recursion. This is what you should actually use in practice.

def ctz(x): # Count trailing zeroes
    return (x & -x).bit_length() - 1
 
def shallowest_decomposition_tree(graph, root=0):
    n = len(graph)
    forbidden = [0] * n
    decomposition_tree = [[] for _ in range(n)]
    stacks = [[] for _ in range(n.bit_length())]
    def extract_chain(labels, u):
        while labels:
            label = labels.bit_length() - 1
            labels ^= 2**label
            v = stacks[label].pop()
            decomposition_tree[u].append(v)
            u = v
    dfs = [root]
    while dfs:
        u = dfs.pop()
        if u >= 0:
            forbidden[u] = -1
            dfs.append(~u)
            for v in graph[u]:
                if not forbidden[v]:
                    dfs.append(v)
        else:
            u = ~u
            forbidden_once = forbidden_twice = 0
            for v in graph[u]:
                forbidden_twice  |= forbidden_once & (forbidden[v] + 1)
                forbidden_once  |= forbidden[v] + 1
            forbidden[u] = forbidden_once | (2**forbidden_twice.bit_length() - 1)    
            label_u = ctz(forbidden[u] + 1)
            stacks[label_u].append(u)
            for v in graph[u]: 
                extract_chain((forbidden[v] + 1) & (2**label_u - 1), u)
    max_label = (forbidden[root] + 1).bit_length() - 1
    decomposition_root = stacks[max_label].pop()
    extract_chain((forbidden[root] + 1) & (2**max_label - 1), decomposition_root)
    return decomposition_tree, decomposition_root

#define remover(vec, v) vec.erase(find(vec.begin(), vec.end(), v))
 
// Rooted tree
struct arborescence {
    std::vector<std::vector<int>> children;
    int root;
};
 
arborescence centroid_decomposition_tree(std::vector<std::vector<int>> graph) {
    int n = graph.size();
    std::vector<int> subtree_size(n);
    auto dfs = [&](int u, auto&& self) -> int {
        for (auto v : graph[u]) {
            remover(graph[v], u);
            subtree_size[u] += self(v, self);
        }
        return ++subtree_size[u];
    };
    dfs(0, dfs);
    
    std::vector<std::vector<int>> decomposition_tree(n);
    auto centroid_reroot = [&](int u, auto&& self) -> int {
        int N = subtree_size[u];
        for (bool found = 0; !found;) {
            found = 1;
            for (auto v : graph[u]) {
                if (subtree_size[v] > N / 2) {
                    found = 0;
                    subtree_size[u] = N - subtree_size[v];
                    remover(graph[u], v);
                    graph[v].push_back(u);
                    u = v;
                    break;
                }
            }
        }
        for (auto v : graph[u]) {
            decomposition_tree[u].push_back(self(v, self));
        }
        return u;
    };
    int decomposition_root = centroid_reroot(0, centroid_reroot);
    return {std::move(decomposition_tree), decomposition_root};
}

def centroid_decomposition_tree(graph):
    n = len(graph)
    graph = [c[:] for c in graph] # copy
    bfs = [0]
    for node in bfs:
        bfs += graph[node]
        for nei in graph[node]:
            graph[nei].remove(node)
    size = [0] * n
    for node in reversed(bfs):
        size[node] = 1 + sum(size[child] for child in graph[node])
    decomposition_tree = [[] for _ in range(n)]
    def centroid_reroot(u):
        N = size[u]
        while True:
            for v in graph[u]:
                if size[v] > N // 2:
                    size[u] = N - size[v]
                    graph[u].remove(v)
                    graph[v].append(u)
                    u = v
                    break
            else: # u is the centroid
                decomposition_tree[u] = [centroid_reroot(v) for v in graph[u]]
                return u
    decomposition_root = centroid_reroot(0)
    return decomposition_tree, decomposition_root

Suppose you are given a tree, where there is a treasure hidden at some node on this tree. Initially, you do not know where the treasure is hidden, but you are allowed to make guesses of where the treasure is hidden. If you guess that the treasure is at some node $$$u$$$, then

• if the treasure is at node $$$u$$$, then you have found the treasure and you are done.
• if the treasure is not at node $$$u$$$, then you will be told the direction of the treasure, i.e. which neighbour of u is in the direction of the treasure.

The goal is to find the treasure using the fewest possible number of guesses (in the worst case).

This tree is an example where centroid decomposition has depth $$$4$$$, but there exists a decomposition tree with depth $$$3$$$.

The red node marks the centroid of the tree. Note that the centroid guessing strategy requires $$$4$$$ guesses in the worst case to find the treasure. However, by guessing the blue node first, you only need 3 guesses to find the treasure.

Construct a tree of $$$n$$$ nodes consisting of several stars, where $$$n$$$ is some power of two minus one. Make the first star out of $$$(n + 1)/2$$$ nodes, the second star out of $$$(n + 1)/4$$$ nodes, the 3rd star out of $$$(n + 1)/8$$$ nodes, etc... Finally, the smallest star should have exactly $$$2$$$ nodes in it. There will now be one node left over. Connect that node to the centers of all of the stars.

Here is an example of the starry tree for $$$n=15$$$. It consists of $$$3$$$ stars and a top node connecting the centers of all of the stars.

Note that the centroid decomposition will split the stars in the order of largest star to smallest star. So the centroid decomposition tree has depth $$$O(\log n)$$$. But if you first split at the top node, you will get a decomposition of depth $$$3$$$. So this graph is an example of where the centroid decomposition can be beaten by a factor of $$$O(\log(n))$$$.

To create this tree, start with a single node. Then from this starting node, build out paths of odd lengths $$$(3,5,7, ...)$$$.

Here is an example with 4 paths going out from the starting node (which is placed at the top).

Note that the center of this graph is the right most child of the starting node. This results in a really unbalanced decomposition tree with a depth of $$$O(\sqrt{n})$$$, which is far from optimal. In comparison, the centroid decomposition tree has a depth of $$$O(\log{n})$$$.

Given a tree. A labeling of the nodes in the tree is called valid if

1. Each node is labeled by a non-negative integer.

2. If two nodes $$$u$$$ and $$$v$$$ have the same label $$$L = \text{label}(u) = \text{label}(v)$$$, then there must exist some node $$$w$$$ on the path between $$$u$$$ and $$$v$$$ such that $$$\text{label}(w) > L$$$.

The goal is to as few distinct labels as possible.

Note that the version of this problem up on Codeforces (321C - Ciel the Commander) simply requires the usage of at most $$$26$$$ different labels, making it so centroid decomposition is sufficient to solve the problem. But the goal for this blog is to instead solve this labeling problem optimally, using as few distinct labels as possible.

Start by picking some arbitrary node as the root. Then label the nodes bottom up, starting from the leaves. Label each node greedily by picking the smallest possible allowed label. I.e. give $$$u$$$ the smallest possible label such that there is no contradiction with the labels of the descendants of $$$u$$$.

Remark: The greedy labeling does depend on the choice of the root, so the choice of the root is not completely arbitrary. However, as we will see in Section 3.2, the number of distinct labels used is independent of the choice of the root. So the choice of the root does not matter in the end.

Recall the DP-algorithm used to compute forbidden(u) in the case of greedy labeling. It computes which labels are forbidden for $$$u$$$ based on the values of forbidden(v) + 1 over all children $$$v$$$ of $$$u$$$. Note that this DP-algorithm can be modified to compute forbidden(u) for any labeling (and not just the greedy labeling). To do this, instead of always adding + 1, we need to add some other number $$$b_v \geq 1$$$. The value of $$$b_v$$$ depends only on the label of $$$v$$$ and the value of forbidden(v).

Here is a generalized DP-algorithm that given a root and a labeling computes forbidden(u) for all nodes $$$u$$$:

def compute_b(label_v, forbidden_v):
    assert 2**label_v & forbidden_v == 0
    b_v = ((forbidden_v | 2**label_v) & ~(2**label_v - 1)) - forbidden_v
    assert b_v >= 1
    # Note: if v was labeled greedily, then b_v will be 1
    # A non-greedy choice of labeling v will make b_v > 1.
    return b_v

forbidden = [0] * n
def dfs(u, p):
    forbidden_once = forbidden_twice = 0
    for v in graph[u]:
        if v != p:
            dfs(v, u)
            forbidden_by_v = forbidden[v] + compute_b(label[v], forbidden[v])
            forbidden_twice |= forbidden_once & forbidden_by_v
            forbidden_once |= forbidden_by_v
    forbidden[u] = forbidden_once | (2**forbidden_twice.bit_length() - 1)
dfs(root, -1)

Now all that remains to prove the Lemma is to show that forbidden(u) is a monotone multivariate function as a function of forbidden(v_1) + b_{v_1}, ...,forbidden(v_m) + b_{v_m} over the children $$$v_1$$$,...,$$$v_m$$$ of $$$u$$$, as monotonicity would imply that the greedy labeling (picking all $$$b_v$$$'s to be $$$1$$$) minimizes forbidden(u). So we want to prove that the following multivariate function $$$f(x_1, ..., x_m)$$$ is monotone

def f(X):
    forbidden_once = forbidden_twice = 0
    for x in X:
        forbidden_twice |= forbidden_once & x
        forbidden_once |= x
    return forbidden_once | (2**forbidden_twice.bit_length() - 1)

I don't know of a simple argument to show this. However, one way to see that $$$f$$$ is monotone think of the input variables x_1, ..., x_m as being 0/1 vectors v_1, ..., v_m in a real vector space. Additionally, also think of $$$f(X) - 1$$$ as a 0/1 vector. In this view, $$$f(X) - 1$$$ is the lexographically smallest 0/1 vector $$$\geq v_1 + ... + v_m$$$. From this it is possible to draw the conclusion that $$$f$$$ must be monotone. If anyone has a cleaner argument for this, then please post it below.

Consider the following tree. Here I've picked the top node to be the root, and then I've greedily labeled the nodes from bottom to top.

In the following picture, I've drawn all of the different chains found in this tree (I've used one color per chain). There are five chains in total. Every node that has an ancestor with a higher label than its label, is part of the chain that starts at its lowest higher labeled ancestor. Each chain is then sorted in decreasing label order.

Using the chains above, we can easily build the shallowest decomposition tree of the original tree. In the following figure, showing the shallowest decomposition tree, I've kept the same colors as I used for the chains in the previous figure. As you can see, the edges in the decomposition tree have a 1-to-1 correspondence with the chain edges.

In the previous example, the root turned out to be the highest labeled node. But this is not always the case. In the following tree, the root (at the top) has label $$$0$$$, but the maximum label in the tree has label $$$3$$$.

Because the root doesn't have the highest label. On top of the chains going "downwards" (as in the previous example), there will also be a chain going upwards, consisting of the set of all nodes that are missing a higher labeled ancestor.

Finally this is the shallowest decomposition tree of the tree.

Each chain is a subset of nodes with distinct labels. A node $$$w$$$ is part of a chain starting at a node $$$u$$$, if $$$u \neq w$$$ and either

Case 1. $$$u$$$ is the lowest higher labeled ancestor of $$$w$$$.

Case 2. $$$u$$$ is the highest labeled node and $$$w$$$ is missing a higher labeled ancestor.

Furthermore, if two nodes share the same lowest higher labeled ancestor $$$u$$$ (Case 1). Then they are part of the same chain if and only if they are descendants of the same child $$$v$$$ of $$$u$$$. On the other hand, if two nodes are missing a higher labeled ancestor (Case 2), then they are always part of the same chain.

def ctz(x): # Count trailing zeros
    return (x & -x).bit_length() - 1

def shallowest_decomposition_tree(graph, root=0):
    n = len(graph)
  
    # Create a greedy labeling
    forbidden = [0] * n
    def dfs(u, p):
        forbidden_once = forbidden_twice = 0
        for v in graph[u]:
            if v != p:
                dfs(v, u)
                forbidden_by_v = forbidden[v] + 1
                forbidden_twice |= forbidden_once & forbidden_by_v
                forbidden_once |= forbidden_by_v
        forbidden[u] = forbidden_once | (2**forbidden_twice.bit_length() - 1)
    dfs(root, -1)
    labels = [ctz(forbidden[u] + 1) for u in range(n)]  
  
    # Extract the shallowest decomposition tree in O(n) time
    decomposition_tree = [[] for _ in range(n)]
    stacks = [[] for _ in range(max(labels) + 1)]
    def extract_chain(labels, u):
        # Extract a chain starting at node u
        while labels:
            label = labels.bit_length() - 1
            labels ^= 2**label
            v = stacks[label].pop()
            decomposition_tree[u].append(v)
            u = v
    def dfs(u, p):
        stacks[labels[u]].append(u)
        for v in graph[u]:
            if v != p: 
                dfs(v, u)
                extract_chain((forbidden[v] + 1) & (2**labels[u] - 1), u)
    dfs(root, -1)
    decomposition_root = stacks[-1].pop()
    extract_chain((forbidden[root] + 1) & (2**labels[decomposition_root] - 1), decomposition_root)
    return decomposition_tree, decomposition_root

def dfs(node, parent=-1):
  nodeDP = 0
  for nei in graph[node]:
    if nei != parent:
      nodeDP = nodeDP + max(dfs(nei, node), 0)
  return nodeDP + 2 * color[node] - 1
print(dfs(u))

def treeDP(root, graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
  DP = [0] * len(graph)
  def dfs(node, parent=-1):
    nodeDP = default[node]
    for eind, nei in enumerate(graph[node]):
      if nei != parent:
        neiDP = dfs(nei, node)
        nodeDP = combine(nodeDP, neiDP, node, eind)
    parent_eind = -1 if parent == -1 else graph[node].index(parent)
    DP[node] = finalize(nodeDP, node, parent_eind)
    return DP[node]
  
  dfs(root)
  return DP

default = [0] * n

def combine(nodeDP, neiDP, node, eind):
  return nodeDP + max(neiDP, 0)

def finalize(nodeDP, node, eind):
  return nodeDP + 2 * color[node] - 1

print(treeDP(u, graph, default, combine, finalize)[u])

default = [0] * n

def combine(nodeDP, neiDP, node, eind):
  return nodeDP + max(neiDP, 0)

def finalize(nodeDP, node, eind):
  return nodeDP + 2 * color[node] - 1

rootDP, _, _ = rerooter(graph, default, combine, finalize)
print(*rootDP)

default = [0] * n

def combine(nodeDP, neiDP, node, eind):
  return nodeDP + max(neiDP + 2 * edge_color[node][eind] - 1, 0)

def finalize(nodeDP, node, eind):
  return nodeDP + 2 * color[node] - 1

rootDP, _, _ = rerooter(graph, default, combine, finalize)
print(*rootDP)

def treeDP(root, graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
  DP = [0] * len(graph)
  def dfs(node, parent=-1):
    nodeDP = default[node]
    for eind, nei in enumerate(graph[node]):
      if nei != parent:
        neiDP = dfs(nei, node)
        nodeDP = combine(nodeDP, neiDP, node, eind)
    parent_eind = -1 if parent == -1 else graph[node].index(parent)
    DP[node] = finalize(nodeDP, node, parent_eind)
    return DP[node]
  
  dfs(root)
  return DP

def rerooter(graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
  n = len(graph)
  rootDP = []
  edgeDP = {}
  for root in range(n):
    DP = treeDP(root, graph, default, combine, finalize)
    rootDP.append(DP[root])
    for nei in graph[root]:
      edgeDP[root, nei] = DP[nei]
    
  forwardDP = [[edgeDP[node, nei] for nei in graph[node]] for node in range(n)]
  reverseDP = [[edgeDP[nei, node] for nei in graph[node]] for node in range(n)]

  return rootDP, forwardDP, reverseDP 

def exclusive(A, zero, combine, node):
    n = len(A)

    exclusiveA = [zero] * n
    for i in range(n):
        for j in range(n):
            if i != j:
                exclusiveA[i] = combine(exclusiveA[i], A[j], node, j)
    return exclusiveA

def treeDP(root, graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
  DP = [0] * len(graph)
  def dfs(node, parent=-1):
    nodeDP = default[node]
    for eind, nei in enumerate(graph[node]):
      if nei != parent:
        neiDP = dfs(nei, node)
        nodeDP = combine(nodeDP, neiDP, node, eind)
    parent_eind = -1 if parent == -1 else graph[node].index(parent)
    DP[node] = finalize(nodeDP, node, parent_eind)
    return DP[node]
  
  dfs(root)
  return DP

def rerooter(graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
    n = len(graph)
    rootDP = [0] * n
    forwardDP = [None] * n
    reverseDP = [None] * n
 
    # Compute DP for root=0
    DP = treeDP(0, graph, default, combine, finalize)
     
    # Do rerooting using the exclusive function
    def dfs(node, parent=-1, parent_eind=-1):
        forwardDP[node] = [DP[nei] for nei in graph[node]]
        rerootDP = exclusive(forwardDP[node], default[node], combine, node)
        reverseDP[node] = [finalize(nodeDP, node, eind) for eind, nodeDP in enumerate(rerootDP)]
        nodeDP = combine(rerootDP[0], forwardDP[node][0], node, 0) if n > 1 else default[node]
        rootDP[node] = finalize(nodeDP, node, -1)
        for eind, nei in enumerate(graph[node]):
            if nei != parent:
                DP[node] = reverseDP[node][eind]
                dfs(nei, node, eind)
    dfs(0)
    return rootDP, forwardDP, reverseDP

Each node in the exclusive segment tree should accumulate everything outside of its interval. For example, in the following figure, the node marked with a red cross should accumulate the 4 values in the bottom right.

Likewise, the node marked here with a red cross should accumulate 6 values.

def exclusive(A, zero, combine, node):
    n = len(A)
    exclusiveA = [zero] * n # Exclusive segment tree
 
    # Build exclusive segment tree
    for bit in range(n.bit_length())[::-1]:
        for i in range(n)[::-1]:
            # Propagate values down the segment tree    
            exclusiveA[i] = exclusiveA[i // 2]
        for i in range(n & ~(bit == 0)):
            # Fold A[i] into exclusive segment tree
            ind = (i >> bit) ^ 1
            exclusiveA[ind] = combine(exclusiveA[ind], A[i], node, i)
    return exclusiveA

def treeDP(root, graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
  DP = [0] * len(graph)
  def dfs(node, parent=-1):
    nodeDP = default[node]
    for eind, nei in enumerate(graph[node]):
      if nei != parent:
        neiDP = dfs(nei, node)
        nodeDP = combine(nodeDP, neiDP, node, eind)
    parent_eind = -1 if parent == -1 else graph[node].index(parent)
    DP[node] = finalize(nodeDP, node, parent_eind)
    return DP[node]
  
  dfs(root)
  return DP

def rerooter(graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
    n = len(graph)
    rootDP = [0] * n
    forwardDP = [None] * n
    reverseDP = [None] * n
 
    # Compute DP for root=0
    DP = treeDP(0, graph, default, combine, finalize)
     
    # Do rerooting using the exclusive function
    def dfs(node, parent=-1, parent_eind=-1):
        forwardDP[node] = [DP[nei] for nei in graph[node]]
        rerootDP = exclusive(forwardDP[node], default[node], combine, node)
        reverseDP[node] = [finalize(nodeDP, node, eind) for eind, nodeDP in enumerate(rerootDP)]
        nodeDP = combine(rerootDP[0], forwardDP[node][0], node, 0) if n > 1 else default[node]
        rootDP[node] = finalize(nodeDP, node, -1)
        for eind, nei in enumerate(graph[node]):
            if nei != parent:
                DP[node] = reverseDP[node][eind]
                dfs(nei, node, eind)
    dfs(0)
    return rootDP, forwardDP, reverseDP

def exclusive(A, zero, combine, node):
    n = len(A)
    exclusiveA = [zero] * n # Exclusive segment tree
 
    # Build exclusive segment tree
    for bit in range(n.bit_length())[::-1]:
        for i in range(n)[::-1]:
            # Propagate values down the segment tree    
            exclusiveA[i] = exclusiveA[i // 2]
        for i in range(n & ~(bit == 0)):
            # Fold A[i] into exclusive segment tree
            ind = (i >> bit) ^ 1
            exclusiveA[ind] = combine(exclusiveA[ind], A[i], node, i)
    return exclusiveA
 
def rerooter(graph, default, combine, finalize = lambda nodeDP,node,eind: nodeDP):
    n = len(graph)
    rootDP = [0] * n
    forwardDP = [None] * n
    reverseDP = [None] * n
 
    # Compute DP for root=0
    DP = [0] * n
    bfs = [0]
    P = [0] * n
    for node in bfs:
        for nei in graph[node]:
            if P[node] != nei:
                P[nei] = node
                bfs.append(nei)
 
    for node in reversed(bfs):
        nodeDP = default[node]
        for eind, nei in enumerate(graph[node]):
            if P[node] != nei:
                nodeDP = combine(nodeDP, DP[nei], node, eind)
        DP[node] = finalize(nodeDP, node, graph[node].index(P[node]) if node else -1)
    # DP for root=0 done
    
    # Use the exclusive function to reroot 
    for node in bfs:
        DP[P[node]] = DP[node]
        forwardDP[node] = [DP[nei] for nei in graph[node]]
        rerootDP = exclusive(forwardDP[node], default[node], combine, node)
        reverseDP[node] = [finalize(nodeDP, node, eind) for eind, nodeDP in enumerate(rerootDP)]
        rootDP[node] = finalize((combine(rerootDP[0], forwardDP[node][0], node, 0) if n > 1 else default[node]), node, -1)
        for nei, dp in zip(graph[node], reverseDP[node]):
            DP[nei] = dp
    return rootDP, forwardDP, reverseDP

Note that \begin{equation} A(x) = \frac{1}{2} \left(1 + \frac{x^{n/2}}{\sqrt{c}}\right) A(x) + \frac{1}{2} \left(1 — \frac{x^{n/2}}{\sqrt{c}}\right) A(x). \end{equation}

Let $$$Q^-(x)$$$ denote the quotient of $$$A(x)$$$ divided by $$$(x^n/2 - \sqrt{c})$$$ and let $$$Q^+(x)$$$ denote the quotient of $$$A(x)$$$ divided by $$$(x^n/2 + \sqrt{c})$$$. Then

and

With this we have shown that

Here $$$A(x)$$$ is expressed as remainder + quotient times $$$(x^n - c)$$$. So we have proven the formula.

"""
Calculates P(x) * Q(x) % (x^n - c) in O(n log n) time

Input:
  n: Integer, needs to be power of 2
  c: Non-zero complex floating point number
  P: A list of length n representing a polynomial P(x) % (x^n - c)
  Q: A list of length n representing a polynomial Q(x) % (x^n - c)
Output:
  A list of length n representing the polynomial P(x) * Q(x) % (x^n - c)
"""
def fast_polymult_mod(P, Q, n, c):
    assert len(P) == n and len(Q) == n
    
    # Base case
    if n == 1:
        return [P[0] * Q[0]]

    assert n % 2 == 0
    import cmath
    sqrtc = cmath.sqrt(c)

    # Calulate P_minus := P mod (x^(n/2) - sqrt(c))
    #          Q_minus := Q mod (x^(n/2) - sqrt(c))

    P_minus = [p1 + sqrtc * p2 for p1,p2 in zip(P[:n//2], P[n//2:])]
    Q_minus = [q1 + sqrtc * q2 for q1,q2 in zip(Q[:n//2], Q[n//2:])]

    # Calulate P_plus := P mod (x^(n/2) + sqrt(c))
    #          Q_plus := Q mod (x^(n/2) + sqrt(c))

    P_plus = [p1 - sqrtc * p2 for p1,p2 in zip(P[:n//2], P[n//2:])]
    Q_plus = [q1 - sqrtc * q2 for q1,q2 in zip(Q[:n//2], Q[n//2:])]

    # Recursively calculate PQ_minus := P * Q % (x^n/2 - sqrt(c)) 
    #                       PQ_plus  := P * Q % (x^n/2 + sqrt(c))
    
    PQ_minus = fast_polymult_mod(P_minus, Q_minus, n//2, sqrtc)
    PQ_plus  = fast_polymult_mod(P_plus,  Q_plus,  n//2, -sqrtc)

    # Calculate PQ mod (x^n - c) using PQ_minus and PQ_plus
    PQ = [(m + p)/2         for m,p in zip(PQ_minus, PQ_plus)] +\
         [(m - p)/(2*sqrtc) for m,p in zip(PQ_minus, PQ_plus)]
    
    return PQ

"""
Calculates P(x) * Q(x)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x)
"""
def fast_polymult(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    n = 1
    while n < res_len:
        n *= 2

    # Pad with extra 0s to reach length n
    P = P + [0] * (n - n1)
    Q = Q + [0] * (n - n2)
    
    # Pick non-zero c arbitrarily =)
    c = 123.24

    # Calculate P*Q mod x^n - c
    PQ = fast_polymult_mod(P, Q, n, c)

    # Remove extra 0 padding and return
    return PQ[:res_len]

We don't actually need the assumption that $$$n$$$ is a power of 2. If $$$n$$$ ever becomes odd during the recursion, then we have two choices: Either fall back to a $$$O(n^2)$$$ algorithm or fall back to the unmodded $$$O(n \log{n})$$$ Polynomial multiplication algorithm.

Let us discuss the run time of falling back to the $$$O(n^2)$$$ algorithm when $$$n$$$ becomes odd. Assume that $$$n = a \cdot 2^b$$$, where $$$a$$$ is an odd integer and $$$b$$$ is an integer. Think of the recursive algorithm as having layers, one layer for each possible value of $$$n$$$. The first $$$b$$$ layers will all take $$$O(n)$$$ time each. In the $$$(b+1)$$$-th layer the value of $$$n$$$ is $$$a$$$. Using the $$$O(n^2)$$$ polynomial multiplication algorithm leads to this layer taking $$$O(n/a \cdot a^2) = O(n \cdot a)$$$ time. The final time complexity comes out to be $$$O((a + b) \, n)$$$.

"""
Calculates P(x) * Q(x) % (x^n - c) in O((a + b) * n) time, where n = a*2^b.

Input:
  n: Integer
  c: Non-zero complex floating point number
  P: A list of length n representing a polynomial P(x) % (x^n - c)
  Q: A list of length n representing a polynomial Q(x) % (x^n - c)
Output:
  A list of length n representing the polynomial P(x) * Q(x) % (x^n - c)
"""
def fast_polymult_mod2(P, Q, n, c):
    assert len(P) == n and len(Q) == n
    
    # Base case (n is odd)
    if n & 1:
        # Calculate the answer in O(n^2) time
        res = [0] * (2*n)
        for i in range(n):
            for j in range(n):
                res[i + j] += P[i] * Q[j]
        return [r1 + c * r2 for r1,r2 in zip(res[:n], res[n:])]

    assert n % 2 == 0
    import cmath
    sqrtc = cmath.sqrt(c)

    # Calulate P_minus := P mod (x^(n/2) - sqrt(c))
    #          Q_minus := Q mod (x^(n/2) - sqrt(c))

    P_minus = [p1 + sqrtc * p2 for p1,p2 in zip(P[:n//2], P[n//2:])]
    Q_minus = [q1 + sqrtc * q2 for q1,q2 in zip(Q[:n//2], Q[n//2:])]

    # Calulate P_plus := P mod (x^(n/2) + sqrt(c))
    #          Q_plus := Q mod (x^(n/2) + sqrt(c))

    P_plus = [p1 - sqrtc * p2 for p1,p2 in zip(P[:n//2], P[n//2:])]
    Q_plus = [q1 - sqrtc * q2 for q1,q2 in zip(Q[:n//2], Q[n//2:])]

    # Recursively calculate PQ_minus := P * Q % (x^n/2 - sqrt(c)) 
    #                       PQ_plus  := P * Q % (x^n/2 + sqrt(c))
    
    PQ_minus = fast_polymult_mod2(P_minus, Q_minus, n//2, sqrtc)
    PQ_plus  = fast_polymult_mod2(P_plus,  Q_plus,  n//2, -sqrtc)

    # Calculate PQ mod (x^n - c) using PQ_minus and PQ_plus
    PQ = [(m + p)/2         for m,p in zip(PQ_minus, PQ_plus)] +\
         [(m - p)/(2*sqrtc) for m,p in zip(PQ_minus, PQ_plus)]
    
    return PQ

The reason why this is super useful is that it allows us to speed up the fast unmodded polynomial multiplication algorithm. As long as we are fine with $$$a$$$ being less than say $$$10$$$, then we might be able to choose a significantly smaller $$$n$$$ compared to what would be possible if we were allowed to only choose powers of two. This trick has the potential of making the fast unmodded polynomial multiplication algorithm run twice as fast.

"""
Calculates P(x) * Q(x)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x)
"""
def fast_polymult2(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    b = 0
    alim = 10
    while alim * 2**b < res_len:
        b += 1
    a = (res_len - 1) // 2**b + 1
    n = a * 2**b

    # Pad with extra 0s to reach length n
    P = P + [0] * (n - n1)
    Q = Q + [0] * (n - n2)
    
    # Pick non-zero c arbitrarily =)
    c = 123.24

    # Calculate P*Q mod x^n - c
    PQ = fast_polymult_mod2(P, Q, n, c)

    # Remove extra 0 padding and return
    return PQ[:res_len]

Suppose that $$$P(x)$$$ and $$$Q(x)$$$ are two real polynomial, and that we want to calculate $$$P(x) \, Q(x)$$$. As discussed earlier, we can calculate the unmodded polynomial product by picking $$$n$$$ large enough such that $$$(P(x) \, Q(x)) \% (x^n - c) = P(x) \, Q(x)$$$ (here $$$c$$$ is any non-zero complex number), and then running the divide and conquer algorithm. But it turns out there is something smarter that we can do.

If we use $$$c = \text{i}$$$ (the imaginary unit) as the inital value of $$$c$$$, then this will allow us to pick an even smaller value for $$$n$$$. The reason for this is that if we get "overflow" from $$$n$$$ being too small, then that overflow will be placed into the imaginary part of the result $$$(P(x) \, Q(x)) \% (x^n - \text{i})$$$. This means that by using $$$c = \text{i}$$$ we are allowed to to pick $$$n$$$ as half the size compared to if we weren't using $$$c=\text{i}$$$. So this trick speeds the fast unmodded polynomial multiplication algorithm up by exactly a factor of 2.

There is somewhat well known technique called "reweighting" that allows us to switch between working with $$$\% (x^n - c)$$$ and working with $$$\% (x^n - 1)$$$. I've previously written a blog explaining this technique, see here. The technique uses the (2^b)-th root of $$$c$$$ to switch variables from $$$\frac{x^n}{c}$$$ to $$$y^n$$$, where $$$n = a \cdot 2^b$$$.

So why would we be interested in switching from $$$\% (x^n - c)$$$ to $$$\% (x^n - 1)$$$? The reason is that by using $$$c=1$$$, we don't need to bother with multiplying or dividing with $$$c$$$ or $$$\sqrt{c}$$$ anywhere, since $$$c=\sqrt{c}=1$$$. Additionally, if $$$c=-1$$$ or $$$c=\text{i}$$$ or $$$c=\text{-i}$$$, then multiplying or dividing by $$$c$$$ can be done very efficiently. So whenever $$$c$$$ becomes something other than $$$1,-1,\text{i}$$$ or $$$-\text{i}$$$, then it makes sense to use the reweight trick to switch back to $$$c=1$$$. This will significantly reduce the number of floating point operations used by the algorithm. Fewer floating point operations means that the algorithm both has the potential to be faster and more nummerically stable. So reweighting is definitely something to consider if you want to create a heavily optimized polynomial multiplication implementation.

In the tail of the recursion (i.e. when $$$n$$$ reaches 1), you are calculating $$$P(x) \, Q(x) \% (x - c)$$$, for some non-zero complex number $$$c$$$. This is infact the same thing as evaluating the polynomial $$$P(x) \, Q(x)$$$ at $$$x=c$$$. Furthermore, if you initially started with $$$c=1$$$, then the $$$c$$$ in the tail will be some $$$n$$$-th root of unity. If you analyze it more carefully, then you will see that each tail corresponds to a different $$$n$$$-th root of unity. So what the algorithm is actually doing is evaluating $$$P(x) \, Q(x)$$$ in all possible $$$n$$$-th roots of unity.

The $$$n$$$-th order FFT of a polynomial is defined as the polynomial evaluated in all $$$n$$$-th roots of unity. This means that the algorithm is infact an FFT algorithm. However, if you want to use it to calculate FFT, then make sure you order the $$$n$$$-th roots of unity according to the standard order used for FFT algorithms. The standard order is $$$\exp{(\frac{2 \pi \text{i}}{n} 0)}, \exp{(\frac{2 \pi \text{i}}{n} 1)}, ..., \exp{(\frac{2 \pi \text{i}}{n} (n-1))}$$$. To get the ordering correct, you will probably need to do a "bit reversal" at the end.

The Cooley-Tukey algorithm is the standard algorithm for calculating FFT. It is for exmple used in this blog [Tutorial] FFT by -is-this-fft-. The idea behind the algorithm is to split up the polynomial $$$P(x)$$$ into an even part $$$P_{\text{even}}(x^2)$$$ and an odd part $$$x \, P_{\text{odd}}(x^2)$$$. You can calculate the FFT of $$$P(x)$$$ using the FFTs of $$$P_{\text{even}}(x)$$$ and $$$P_{\text{odd}}(x)$$$. So Cooley-Tukey is a $$$O(n \log{n})$$$ divide and conquer algorithm that repeatedly splits up the polynomial into odd and even parts.

The wiki article for Cooley-Tukey has a nice description of the algorithm

If you compare this to calculating FFT using the divide and conquer polynomial mod method you instead get

where $$$c$$$ is an $$$n$$$-th root of unity that is independent of $$$k$$$. This is very different compared to Cooley-Tukey since $$$c$$$ doesn't have a dependence on $$$k$$$ unlike $$$e^{- \frac{2 \pi \text{i}}{n} k}$$$. Infact, $$$c$$$ being constant means that the polynomial mod method has the potential to be faster than Cooley-Tukey.

Here is an FFT implementation. A codegolfed version of the same code can be found on Pyrival.

"""
Calculates FFT(P) in O(n log n) time.

This implementation is based on the 
polynomial modulo multiplication algorithm.

Input:
  P: A list of length n representing a polynomial P(x).
     n needs to be a power of 2.
Output:
  A list of length n representing the FFT of the polynomial P,
     i.e. the list [P(exp(2j pi / n * i) for i in range(n)]
"""
rt = [1] # List used to store roots of unity
def FFT(P):
    n = len(P)
    # Assert n is a power of 2
    assert n and (n - 1) & n == 0
    # Make a copy of P to not modify original P
    P = P[:] 
    
    # Precalculate the roots
    while 2 * len(rt) < n:
        # 4*len(rt)-th root of unity
        import cmath
        root = cmath.exp(2j * cmath.pi / (4 * len(rt)))
        rt.extend([r * root for r in rt])

    # Transform P
    k = n
    while k > 1:
        for i in range(n//k):
            r = rt[i]
            for j1 in range(i*k, i*k + k//2):
                j2 = j1 + k//2
                z = r * P[j2]
                P[j2] = P[j1] - z
                P[j1] += z
        k //= 2
    
    # Bit reverse P before returning
    rev = [0] * n
    for i in range(1, n):
        rev[i] = rev[i // 2] // 2 + (i & 1) * n // 2

    return [P[r] for r in rev]

# Inverse of FFT(P) using a standard trick
def inverse_FFT(fft_P):
    n = len(fft_P)
    return FFT([fft_P[-i]/n for i in range(n)])

"""
Calculates P(x) * Q(x)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x)
"""
def fast_polymult_using_FFT(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    n = 1
    while n < res_len:
        n *= 2

    # Pad with extra 0s to reach length n
    P = P + [0] * (n - n1)
    Q = Q + [0] * (n - n2)
    
    # Transform P and Q
    fft_P = FFT(P)
    fft_Q = FFT(Q)

    # Calculate FFT of P*Q
    fft_PQ = [p*q for p,q in zip(fft_P,fft_Q)]

    # Inverse FFT
    PQ = inverse_FFT(fft_PQ)
    
    # Remove padding and return
    return PQ[:res_len]
"""
Calculates P(x) * Q(x)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x)
"""
def fast_polymult_using_FFT(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    n = 1
    while n < res_len:
        n *= 2

    # Pad with extra 0s to reach length n
    P = P + [0] * (n - n1)
    Q = Q + [0] * (n - n2)
    
    # Transform P and Q
    fft_P = FFT(P)
    fft_Q = FFT(Q)

    # Calculate FFT of P*Q
    fft_PQ = [p*q for p,q in zip(fft_P,fft_Q)]

    # Inverse FFT
    PQ = inverse_FFT(fft_PQ)
    
    # Remove padding and return
    return PQ[:res_len]

Here is an FTT implementation. It is coded in the same style as in KACTL.

typedef complex<double> C;
typedef vector<double> vd;
void fft(vector<C>& a) {
	int n = sz(a);
	static vector R{1.L + 0il};
	static vector rt{1. + 0i};
	for (static int k = 2; k < n; k *= 2) {
		R.resize(n/2); rt.resize(n/2);
                auto x = pow(1il, 2./k);
		rep(i,k/2,k) rt[i] = R[i] = R[i-k/2] * x;;
	}
	for (int k = n; k > 1; k /= 2) rep(i,0,n/k) rep(j,i*k,i*k + k/2) {
		C &u = a[j], &v = a[j+k/2], &r = rt[i];
		C z(v.real()*r.real() - v.imag()*r.imag(), 
		    v.real()*r.imag() + v.imag()*r.real());
		v = u - z;
		u = u + z;
	}
	vi rev(n);
	rep(i,0,n) rev[i] = rev[i / 2] / 2 + (i & 1) * n / 2;
	rep(i,0,n) if (i < rev[i]) swap(a[i], a[rev[i]]);
}

vd conv(const vd& a, const vd& b) {
	if (a.empty() || b.empty()) return {};
	vd res(sz(a) + sz(b) - 1);
	int L = 32 - __builtin_clz(sz(res)), n = 1 << L;
	vector<C> in(n), out(n);
	copy(all(a), begin(in));
	rep(i,0,sz(b)) in[i].imag(b[i]);
	fft(in);
	for (C& x : in) x *= x;
	rep(i,0,n) out[i] = in[-i & (n - 1)] - conj(in[i]);
	fft(out);
	rep(i,0,sz(res)) res[i] = imag(out[i]) / (4 * n);
	return res;
}

This algorithm requires three properties. Firstly it needs to be possible to divide by $$$2$$$, and secondly $$$\sqrt{c}$$$ needs to exist, and thirdly it needs to be possible to divide by $$$\sqrt{c}$$$. This means that we don't technically need complex numbers, we could also use other number systems (like working modulo a prime or modulo an odd composite number).

Primes that work nicely for this purpose are called "NTT primes", which means that the prime — 1 is divisible by a large power of $$$2$$$. Common examples of NTT primes are: $$$998244353 = 119 \cdot 2^{23} + 1$$$, $$$167772161 = 5 \cdot 2^{25} + 1$$$ and $$$469762049 = 7 \cdot 2^{26} + 1$$$.

One very important remark is that it is possible to calculate $$$P(x) \, Q(x) \% (x^n - c)$$$ by multiplying $$$P(x) \, Q(x)$$$ unmodded, and afterwards modding it with $$$(x^n - c)$$$. The reason why this is useful is that if $$$\sqrt{c}$$$ doesn't exist, then this trick allows us to go back to having $$$c = 1$$$. This can be costly timewise, but it opens up the possibility of using "bad NTT primes"! For example $$$13 \cdot 2^{10} + 1$$$ is a pretty bad NTT prime, but it can still be used by this algorithm to efficiently multiply polynomials.

One of the things I dislike about NTT is that for NTT to be defined, there needs to exist a $$$n$$$-th root of unity. Usually problems involving NTT are designed so that this is never an issue. But if you want to use NTT where it hasn't been designed to magically work, then this is a really big issue. The NTT can become undefined!

Note that this algorithm does not exactly share the same drawback of being undefined. The reason for this is that if $$$\sqrt{c}$$$ doesn't exist, then the algorithm can simply choose to either switch over to a $$$O(n^2)$$$ polynomial multiplication algorithm, or fall back to fast unmodded polynomial multiplication. The implications from this is that this algorithm can do fast modded polynomial multiplication even if it is given a relatively bad NTT prime. I just find this property to be really cool!

A good example of when NTT becomes undefined is this yosup problem convolution_mod_large. Here the NTT mod is $$$998244353 = 119 \cdot 2^{23}$$$. The tricky thing about the problem is that $$$n=2^{24}$$$. Since $$$998244353 = 119 \cdot 2^{23} + 1$$$ there wont exist any $$$n$$$-th root of unity, so the NTT of length $$$n$$$ is undefined. However, the divide and conquer approach from this blog can easily solve the problem by falling back to the $$$O(n^2)$$$ algorithm.

At the time of writing this, I currently hold the record for the fastest solution to Convolution (Large) 2572 ms. But even a simple recursive version is enough to get AC 9292 ms.

Here is an NTT implementation. A codegolfed version of the same code can be found on Pyrival.

# Mod used for NTT
# Requirement: Any odd integer > 2
# It is important that MOD - 1 is
# divisible by lots of 2s
MOD = (119 << 23) + 1
assert MOD & 1

# Precalc non-quadratic_residue (used by the NTT)
non_quad_res = 2
while pow(non_quad_res, MOD//2, MOD) != MOD - 1:
    non_quad_res += 1
rt = [1]

"""
Calculates NTT(P) in O(n log n) time.

This implementation is based on the 
polynomial modulo multiplication algorithm.

Input:
  P: A list of length n representing a polynomial P(x).
     n needs to be a power of 2.
Output:
  A list of length n representing the NTT of the polynomial P,
     i.e. the list [P(root**i) % MOD for i in range(n)]
     where root is an n-th root of unity mod MOD
"""
def NTT(P):
    n = len(P)
    # Assert n is a power of 2
    assert n and (n - 1) & n == 0

    # Check that NTT is defined for this n
    assert (MOD - 1) % n == 0

    # Make a copy of P to not modify original P
    P = P[:] 
    
    # Precalculate the roots
    while 2 * len(rt) < n:
        # 4*len(rt)-th root of unity
        root = pow(non_quad_res, MOD//(4 * len(rt)), MOD)
        rt.extend([r * root % MOD for r in rt])

    # Transform P
    k = n
    while k > 1:
        for i in range(n//k):
            r = rt[i]
            for j1 in range(i*k, i*k + k//2):
                j2 = j1 + k//2
                z = r * P[j2]
                P[j2] = (P[j1] - z) % MOD
                P[j1] = (P[j1] + z) % MOD
        k //= 2
    
    # Bit reverse P before returning
    rev = [0] * n
    for i in range(1, n):
        rev[i] = rev[i // 2] // 2 + (i & 1) * n // 2

    return [P[r] for r in rev]

# Inverse of NTT(P) using a standard trick
def inverse_NTT(ntt_P):
    n = len(ntt_P)
    n_inv = pow(n, -1, MOD) # Requires Python 3.8
    # The following works in any Python version, but requires MOD to be prime
    # n_inv = pow(n, MOD - 2, MOD)
    assert n * n_inv % MOD == 1
    return NTT([ntt_P[-i] * n_inv % MOD for i in range(n)])

"""
Calculates P(x) * Q(x) (where the coeffiecents are returned % MOD)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x) (with coeffients % MOD)
"""
def fast_polymult_using_NTT(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    n = 1
    while n < res_len:
        n *= 2

    # Pad with extra 0s to reach length n
    P = P + [0] * (n - n1)
    Q = Q + [0] * (n - n2)
    
    # Transform P and Q
    ntt_P = NTT(P)
    ntt_Q = NTT(Q)

    # Calculate NTT of P*Q
    ntt_PQ = [p * q % MOD for p,q in zip(ntt_P,ntt_Q)]

    # Inverse NTT
    PQ = inverse_NTT(ntt_PQ)
    
    # Remove padding and return
    return PQ[:res_len]

Here is an NTT implementation. It is coded in the same style as in KACTL.

const ll mod = (119 << 23) + 1;// = 998244353
// For p < 2^30 there is also e.g. 5 << 25, 7 << 26, 479 << 21
// and 483 << 21 The last two are > 10^9.
typedef vector<ll> vl;

#include "../number-theory/ModPow.h"

void ntt(vl &a) {
	int n = sz(a);
	static ll r = 3;
	while(modpow(r, mod/2) + 1 < mod) ++r;
	static vl rt{1};
	for (static int k = 2; k < n; k *= 2) {
		rt.resize(n/2);
                ll x = modpow(r, mod/2/k);
		rep(i,k/2,k) rt[i] = rt[i-k/2] * x % mod;
	}
	for (int k = n; k > 1; k /= 2) rep(i,0,n/k) rep(j,i*k,i*k + k/2) {
		ll &u = a[j], &v = a[j+k/2], z = rt[i] * v % mod;
		v = u - z + (u < z ? mod : 0);
		u = u + z - (u + z >= mod ? mod : 0);
	}
	vi rev(n);
	rep(i,0,n) rev[i] = rev[i / 2] / 2 + (i & 1) * n / 2;
	rep(i,0,n) if (i < rev[i]) swap(a[i], a[rev[i]]);
}
vl conv(vl a, vl b) {
	
	if (a.empty() || b.empty()) return {};
	int s = sz(a) + sz(b) - 1, B = 32 - __builtin_clz(s), n = 1 << B;
	int inv = modpow(n, mod - 2);
	vl out(n);
	a.resize(n); b.resize(n);
	ntt(a), ntt(b);
	rep(i,0,n) out[-i & (n - 1)] = (ll)a[i] * b[i] % mod * inv % mod;
	ntt(out);
	return {out.begin(), out.begin() + s};
}

At the time of writing this, I currently hold the record for the fastest solution to Convolution (Large) 2572 ms. If what you want is a super fast C++ implementation then I would suggest using this implementation.

"""
Calculates P(x) * Q(x) % (x^n - c) in O(n log n) time

Input:
  n: Integer, needs to be power of 2
  c: Non-zero complex floating point number
  P: A list of length 2*n representing a polynomial P(x) % (x^2n - c^2)
  Q: A list of length 2*n representing a polynomial Q(x) % (x^2n - c^2)
Output:
  A list of length n representing the polynomial P(x) * Q(x) % (x^n - c)
"""
def fast_polymult_mod3(P, Q, n, c):
    assert len(P) == 2*n and len(Q) == 2*n

    # Mod P and Q by (x^n - c)
    P = [p1 + c * p2 for p1,p2 in zip(P[:n], P[n:])]
    Q = [q1 + c * q2 for q1,q2 in zip(Q[:n], Q[n:])]
    
    # Base case
    if n == 1:
        return [P[0] * Q[0]]

    assert n % 2 == 0
    import cmath
    sqrtc = cmath.sqrt(c)

    # Recursively calculate PQ_minus := P * Q % (x^n/2 - sqrt(c)) 
    #                       PQ_plus  := P * Q % (x^n/2 + sqrt(c))
    
    PQ_minus = fast_polymult_mod3(P, Q, n//2, sqrtc)
    PQ_plus  = fast_polymult_mod3(P, Q, n//2, -sqrtc)

    # Calculate PQ mod (x^n - c) using PQ_minus and PQ_plus
    PQ = [(m + p)/2         for m,p in zip(PQ_minus, PQ_plus)] +\
         [(m - p)/(2*sqrtc) for m,p in zip(PQ_minus, PQ_plus)]
    
    return PQ
"""
Calculates P(x) * Q(x)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x)
"""
def fast_polymult3(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    n = 1
    while n < res_len:
        n *= 2

    # Pad with extra 0s to reach length 2*n
    P = P + [0] * (2*n - n1)
    Q = Q + [0] * (2*n - n2)
    
    # Pick non-zero c arbitrarily =)
    c = 123.24

    # Calculate P*Q mod x^n - c
    PQ = fast_polymult_mod3(P, Q, n, c)

    # Remove extra 0 padding and return
    return PQ[:res_len]

"""
Calculates P(x) * Q(x) % (x^n - c) in O(n log n) time

Input:
  n: Integer, needs to be power of 2
  c: Non-zero complex floating point number
  P: A list of length 2*n representing a polynomial P(x) % (x^2n - c^2)
  Q: A list of length 2*n representing a polynomial Q(x) % (x^2n - c^2)
Output:
  A list of length n representing the polynomial P(x) * Q(x) % (x^n - c)
"""
def fast_polymult_mod4(P, Q, n, c):
    assert len(P) == 2*n and len(Q) == 2*n

    # Mod P and Q by (x^n - c)
    P = [p1 + c * p2 for p1,p2 in zip(P[:n], P[n:])]
    Q = [q1 + c * q2 for q1,q2 in zip(Q[:n], Q[n:])]
    
    # Base case (n is odd)
    if n & 1:
        # Calculate the answer in O(n^2) time
        res = [0] * (2*n)
        for i in range(n):
            for j in range(n):
                res[i + j] += P[i] * Q[j]
        return [r1 + c * r2 for r1,r2 in zip(res[:n], res[n:])]
    
    assert n % 2 == 0
    import cmath
    sqrtc = cmath.sqrt(c)

    # Recursively calculate PQ_minus := P * Q % (x^n/2 - sqrt(c)) 
    #                       PQ_plus  := P * Q % (x^n/2 + sqrt(c))
    
    PQ_minus = fast_polymult_mod4(P, Q, n//2, sqrtc)
    PQ_plus  = fast_polymult_mod4(P, Q, n//2, -sqrtc)

    # Calculate PQ mod (x^n - c) using PQ_minus and PQ_plus
    PQ = [(m + p)/2         for m,p in zip(PQ_minus, PQ_plus)] +\
         [(m - p)/(2*sqrtc) for m,p in zip(PQ_minus, PQ_plus)]
    
    return PQ

"""
Calculates P(x) * Q(x)

Input:
  P: A list representing a polynomial P(x)
  Q: A list representing a polynomial Q(x)
Output:
  A list representing the polynomial P(x) * Q(x)
"""
def fast_polymult4(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    b = 0
    alim = 10
    while alim * 2**b < res_len:
        b += 1
    a = (res_len - 1) // 2**b + 1
    n = a * 2**b

    # Pad with extra 0s to reach length 2*n
    P = P + [0] * (2*n - n1)
    Q = Q + [0] * (2*n - n2)
    
    # Pick non-zero c arbitrarily =)
    c = 123.24

    # Calculate P*Q mod x^n - c
    PQ = fast_polymult_mod4(P, Q, n, c)

    # Remove extra 0 padding and return
    return PQ[:res_len]

"""
Calculates P(x) * Q(x) % (x^n - c) in O(n log n) time

Input:
  n: Integer, needs to be power of 2
  c: Non-zero complex floating point number
  P: A list of length 2*n representing a polynomial P(x) % (x^2n - c^2)
  Q: A list of length 2*n representing a polynomial Q(x) % (x^2n - c^2)
Output:
  A list of length n representing the polynomial P(x) * Q(x) % (x^n - c)
"""
def fast_polymult_mod4(P, Q, n, c):
    assert len(P) == 2*n and len(Q) == 2*n

    # Mod P and Q by (x^n - c)
    P = [p1 + c * p2 for p1,p2 in zip(P[:n], P[n:])]
    Q = [q1 + c * q2 for q1,q2 in zip(Q[:n], Q[n:])]
    
    # Base case (n is odd)
    if n & 1:
        # Calculate the answer in O(n^2) time
        res = [0] * (2*n)
        for i in range(n):
            for j in range(n):
                res[i + j] += P[i] * Q[j]
        return [r1 + c * r2 for r1,r2 in zip(res[:n], res[n:])]
    
    assert n % 2 == 0
    import cmath
    sqrtc = cmath.sqrt(c)

    # Recursively calculate PQ_minus := P * Q % (x^n/2 - sqrt(c)) 
    #                       PQ_plus  := P * Q % (x^n/2 + sqrt(c))
    
    PQ_minus = fast_polymult_mod4(P, Q, n//2, sqrtc)
    PQ_plus  = fast_polymult_mod4(P, Q, n//2, -sqrtc)

    # Calculate PQ mod (x^n - c) using PQ_minus and PQ_plus
    PQ = [(m + p)/2         for m,p in zip(PQ_minus, PQ_plus)] +\
         [(m - p)/(2*sqrtc) for m,p in zip(PQ_minus, PQ_plus)]
    
    return PQ

"""
Calculates P(x) * Q(x) of two real polynomials

Input:
  P: A list representing a real polynomial P(x)
  Q: A list representing a real polynomial Q(x)
Output:
  A list representing the real polynomial P(x) * Q(x)
"""
def fast_polymult5(P, Q):
    # Calculate length of the list representing P*Q
    n1 = len(P)
    n2 = len(Q)
    res_len = n1 + n2 - 1
    
    # Pick n sufficiently big
    b = 1
    alim = 10
    while alim * 2**b < res_len:
        b += 1
    a = (res_len - 1) // 2**b + 1
    n = a * 2**b
    
    # Pick c = i (imaginary unit)
    c = 1j
    # and decrease the size of n by a factor of 2
    n //= 2

    # Pad with extra 0s to reach length 2*n
    P = P + [0] * (2*n - n1)
    Q = Q + [0] * (2*n - n2)
    
    # Calculate P*Q mod x^n - i
    PQ = fast_polymult_mod4(P, Q, n, c)

    # The imaginary part contains the "overflow"
    PQ = [pq.real for pq in PQ] + [pq.imag for pq in PQ]

    # Remove extra 0 padding and return
    return PQ[:res_len]

Take a closer look at the code =P

The other day, I was coding some C++. My Pythonic muscle memory caused me to accedently write #import instead of #include, and to my great surprise, the code still compiled!

The only information that I was able to find about #import is from this stackoverflow comment. As I understand it, #import is a really old gcc feature that combines #include and #pragma once into one.

There are arguably some other positives to using #import instead of #include. #import is one character shorter than #include, and as someone that mainly uses Python, it is easier for me to remember import than include.

def m_reduce(ab):
  m = ab * Pinv % r
  y = (ab + m * P) // r
  return y if y < P else y - P

1. Switching to floats

What I usually try is switching from using integers to floats. This works if all integers are $$$\leq 2^{52}$$$:

53441394 (Fat TLE and almost MLE) vs 53456822 (AC by a mile)

2. Calculating a * b % MOD

Many problems on CF involve doing multiplication modulo $$$10^9 + 7$$$. This is super slow because of big integers. After a long time of trying around with different methods, I've found that this is the best solution. Today I always use it whenever I do multiplication mod $$$10^9 + 7$$$. But sadly, while it is much faster than using big integers, it is still slow compared to doing a * b % MOD in 64 bit PyPy.

3. When all else fails

Another trick is to use a deprecated feature that only exists in PyPy2:

114602234 (Fat TLE) vs 114611359 (AC by a mile)

A big warning here. Not only is this feature deprecated, but from my own personal experience using it means opening a can of worms. I will only use this as an absolute last resort.

#include <iostream>
using namespace std;
 
double f(double a, double b) {
    return a * a - b;
}

int main() {
  cout.precision(60);
  
  // Calculate 10*10 - 1e-15
  double ans;
  ans = f(atof("10"), atof("1e-15"));
  cout << (double)      ans << '\n';
  cout << (int)         ans << "\n\n";
  
  ans = f(atof("10"), atof("1e-15"));
  cout << (int)         ans << '\n';
  cout << (double)      ans << "\n\n";
  
  ans = f(atof("10"), atof("1e-15"));
  cout << (double)      ans << '\n';
  cout << (long double) ans << "\n\n";
  
  ans = f(atof("10"), atof("1e-15"));
  cout << (long double) ans << '\n';
  cout << (double)      ans << "\n\n";
  return 0;
}

100
100

99
100

100
100

99.99999999999999900079927783735911361873149871826171875
100

100
100

100
100

100
100

100
100

// #pragma GCC target("fpmath=sse,sse2") // Turns off excess precision
// #pragma GCC target("fpmath=387") // Turns on excess precision

#include <iostream>
#include <cstdlib>
#include <cfloat>
using namespace std;

int main() {
  cout << "This is compiled in mode "<< FLT_EVAL_METHOD << '\n';
  cout << "0 means no excess precision.\n";
  cout << "2 means there is excess precision.\n\n";
  
  cout << "The following test detects excess precision\n";
  cout << "0 if no excess precision, or 8e-17 if there is excess precision.\n";
  double a = atof("1.2345678");
  double b = a*a;
  cout << b - 1.52415765279683990130 << '\n';
  return 0;
}

#include <iostream>
#include <cmath>
using namespace std;

bool f() {
    double x;
    cin >> x;
    
    double y = x + 1.0;
    if (y >= 1.0)
        return false;
    
    int w = pow(y, 2);
    
    if (y < 1.0)
        return false;
    return y == 1.0;
}

int main() {
    if (f())
        cout << "HACKED\n";
    else
        cout << "Not hacked\n";
}