Part Ⅳ | Computer Science Study Notes

20 Tries

20.1 R-Way Tries

Tries (from retrieval, but pronounced "try"): Store characters in nodes (not keys), each node has children, one for each possible character.

Trie Search

Follow links corresponding to each character in the key.
Search hit: node where search ends has a non-null value.
Search miss: reach null link or node where search ends has null value.

Trie Delete

Find the node corresponding to key and set value to null.
If node has null value and all null links, remove that node (and recur).

import java.util.HashMap; public class RWayTrie { private static final int R = 256; private final Node root; private static class Node { private boolean isEndOfWord; private final HashMap<Character, Node> children; public Node() { isEndOfWord = false; children = new HashMap<>(); } } public RWayTrie() { root = new Node(); } public void insert(String word) { Node current = root; for (char c : word.toCharArray()) { if (!current.children.containsKey(c)) { current.children.put(c, new Node()); } current = current.children.get(c); } current.isEndOfWord = true; } public boolean search(String word) { Node current = root; for (char c : word.toCharArray()) { if (!current.children.containsKey(c)) { return false; } current = current.children.get(c); } return current.isEndOfWord; } public boolean startsWith(String prefix) { Node current = root; for (char c : prefix.toCharArray()) { if (!current.children.containsKey(c)) { return false; } current = current.children.get(c); } return true; } }

#include <iostream> #include <unordered_map> #include <ranges> constexpr int R = 26; struct Node { bool isEndOfWord; std::unordered_map<char, Node*> children; Node() : isEndOfWord(false) {} }; class RWayTrie { private: Node* root; static void deleteNode(Node* node) { if (node == nullptr) { return; } for (Node* child : node->children | std::views::values) { deleteNode(child); } delete node; } public: RWayTrie() { root = new Node(); } void insert(const std::string& word) const { Node* current = root; for (char c : word) { if (!current->children.contains(c)) { current->children[c] = new Node(); } current = current->children[c]; } current->isEndOfWord = true; } [[nodiscard]] bool search(const std::string& word) const { Node* current = root; for (char c : word) { if (!current->children.contains(c)) { return false; } current = current->children[c]; } return current->isEndOfWord; } [[nodiscard]] bool startsWith(const std::string& prefix) const { Node* current = root; for (char c : prefix) { if (!current->children.contains(c)) { return false; } current = current->children[c]; } return true; } ~RWayTrie() { deleteNode(root); } };

class Node: def __init__(self): self.isEndOfWord = False self.children = {} # Dictionary to store child nodes class RWayTrie: def __init__(self): self.root = Node() def insert(self, word): current = self.root for char in word: if char not in current.children: current.children[char] = Node() current = current.children[char] current.isEndOfWord = True def search(self, word): current = self.root for char in word: if char not in current.children: return False current = current.children[char] return current.isEndOfWord def startsWith(self, prefix): current = self.root for char in prefix: if char not in current.children: return False current = current.children[char] return True

20.2 Ternary Search Tries

Ternary Search Trees

Store characters and values in nodes (not keys).
Each node has 3 children: smaller (left), equal (middle), larger (right).

TST Search

Follow links corresponding to each character in the key.
- If less, take left link; if greater, take right link.
- If equal, take the middle link and move to the next key character.
Search hit & miss:
- Search hit: Node where search ends has a non-null value.
- Search miss: Reach null link or node where search ends has null value.

public class TernarySearchTree { private Node root; private static class Node { char data; boolean isEndOfString; Node left, equal, right; public Node(char data) { this.data = data; this.isEndOfString = false; this.left = null; this.equal = null; this.right = null; } } public TernarySearchTree() { root = null; } public void insert(String word) { root = insertRecursive(root, word, 0); } private Node insertRecursive(Node node, String word, int index) { char c = word.charAt(index); if (node == null) { node = new Node(c); } if (c < node.data) { node.left = insertRecursive(node.left, word, index); } else if (c > node.data) { node.right = insertRecursive(node.right, word, index); } else { if (index < word.length() - 1) { node.equal = insertRecursive(node.equal, word, index + 1); } else { node.isEndOfString = true; } } return node; } public boolean search(String word) { return searchRecursive(root, word, 0); } private boolean searchRecursive(Node node, String word, int index) { if (node == null) { return false; } char c = word.charAt(index); if (c < node.data) { return searchRecursive(node.left, word, index); } else if (c > node.data) { return searchRecursive(node.right, word, index); } else { if (index == word.length() - 1) { return node.isEndOfString; } else { return searchRecursive(node.equal, word, index + 1); } } } public void getWordsWithPrefix(String prefix) { Node node = getPrefixNode(root, prefix, 0); if (node != null) { traverseAndPrint(node, prefix); } } private Node getPrefixNode(Node node, String prefix, int index) { if (node == null) { return null; } char c = prefix.charAt(index); if (c < node.data) { return getPrefixNode(node.left, prefix, index); } else if (c > node.data) { return getPrefixNode(node.right, prefix, index); } else { if (index == prefix.length() - 1) { return node; } else { return getPrefixNode(node.equal, prefix, index + 1); } } } private void traverseAndPrint(Node node, String prefix) { if (node == null) { return; } if (node.isEndOfString) { System.out.println(prefix); } traverseAndPrint(node.left, prefix); traverseAndPrint(node.equal, prefix + node.data); traverseAndPrint(node.right, prefix); } }

#include <iostream> #include <string> class TernarySearchTree { private: struct Node { char data; bool isEndOfString; Node *left, *equal, *right; explicit Node(const char data) : data(data), isEndOfString(false), left(nullptr), equal(nullptr), right(nullptr) {} }; Node *root; static Node* insertRecursive(Node* node, const std::string& word, const int index) { const char c = word[index]; if (node == nullptr) { node = new Node(c); } if (c < node->data) { node->left = insertRecursive(node->left, word, index); } else if (c > node->data) { node->right = insertRecursive(node->right, word, index); } else { if (index < word.length() - 1) { node->equal = insertRecursive(node->equal, word, index + 1); } else { node->isEndOfString = true; } } return node; } static bool searchRecursive(const Node* node, const std::string& word, const int index) { if (node == nullptr) { return false; } const char c = word[index]; if (c < node->data) { return searchRecursive(node->left, word, index); } else if (c > node->data) { return searchRecursive(node->right, word, index); } else { if (index == word.length() - 1) { return node->isEndOfString; } else { return searchRecursive(node->equal, word, index + 1); } } } static Node* getPrefixNode(Node* node, const std::string& prefix, const int index) { if (node == nullptr) { return nullptr; } const char c = prefix[index]; if (c < node->data) { return getPrefixNode(node->left, prefix, index); } else if (c > node->data) { return getPrefixNode(node->right, prefix, index); } else { if (index == prefix.length() - 1) { return node; } else { return getPrefixNode(node->equal, prefix, index + 1); } } } static void traverseAndPrint(const Node* node, const std::string& prefix) { if (node == nullptr) { return; } if (node->isEndOfString) { std::cout << prefix << std::endl; } traverseAndPrint(node->left, prefix); traverseAndPrint(node->equal, prefix + node->data); traverseAndPrint(node->right, prefix); } static void deleteNodes(const Node* node) { if (node == nullptr) { return; } deleteNodes(node->left); deleteNodes(node->equal); deleteNodes(node->right); delete node; } public: TernarySearchTree() : root(nullptr) {} void insert(const std::string& word) { root = insertRecursive(root, word, 0); } [[nodiscard]] bool search(const std::string& word) const{ return searchRecursive(root, word, 0); } void getWordsWithPrefix(const std::string& prefix) const { Node* node = getPrefixNode(root, prefix, 0); if (node != nullptr) { traverseAndPrint(node, prefix); } } ~TernarySearchTree() { deleteNodes(root); } };

class Node: def __init__(self, data): self.data = data self.isEndOfString = False self.left = None self.equal = None self.right = None class TernarySearchTree: def __init__(self): self.root = None def insert(self, word): self.root = self._insert_recursive(self.root, word, 0) def _insert_recursive(self, node, word, index): c = word[index] if node is None: node = Node(c) if c < node.data: node.left = self._insert_recursive(node.left, word, index) elif c > node.data: node.right = self._insert_recursive(node.right, word, index) else: if index < len(word) - 1: node.equal = self._insert_recursive(node.equal, word, index + 1) else: node.isEndOfString = True return node def search(self, word): return self._search_recursive(self.root, word, 0) def _search_recursive(self, node, word, index): if node is None: return False c = word[index] if c < node.data: return self._search_recursive(node.left, word, index) elif c > node.data: return self._search_recursive(node.right, word, index) else: if index == len(word) - 1: return node.isEndOfString else: return self._search_recursive(node.equal, word, index + 1) def get_words_with_prefix(self, prefix): node = self._get_prefix_node(self.root, prefix, 0) if node is not None: self._traverse_and_print(node, prefix) def _get_prefix_node(self, node, prefix, index): if node is None: return None c = prefix[index] if c < node.data: return self._get_prefix_node(node.left, prefix, index) elif c > node.data: return self._get_prefix_node(node.right, prefix, index) else: if index == len(prefix) - 1: return node else: return self._get_prefix_node(node.equal, prefix, index + 1) def _traverse_and_print(self, node, prefix): if node is None: return if node.isEndOfString: print(prefix) self._traverse_and_print(node.left, prefix) self._traverse_and_print(node.equal, prefix + node.data) self._traverse_and_print(node.right, prefix)

TST with branching at root: Hybrid of R-way trie and TST

Do -way branching at root.
Each of root nodes points to a TST.

TST vs. Hashing

TSTs	Hashing
Works only for strings (or digital keys)	Need to examine entire key
Only examines just enough key characters	Search hits and misses cost about the same
Search miss may involve only a few characters	Performance relies on hash function
Supports ordered symbol table operations (plus others!)	Does not support ordered symbol table operations

Implementation	Character Accesses (typical case)
Implementation	Search Hit	Search Miss	Insert	Space (references)
Red-Black BST
Hashing (linear probing)				to
R-Way Trie
TST
TST with R^2

21 Substring Search

21.1 Introduction

Goal: Find pattern of length in text of length (typically ).

Applications:

Find & replace
Computer forensics
Identify patterns indicative of spam
Electronic surveillance
Screen scraping

21.2 Brute-Force Substring Search

Disadvantages

Theoretical challenge: Linear -time guarantee (Worst case: ).
Practical challenge: Avoid backup in text stream (Brute-force algorithm needs backup for every mismatch).

public static int search (String pat, String txt) { int M = pat.length(); int N = txt.length(); int i, j; for (i = 0; i <= N - M; i++) { for (j = 0; j < M; j++) { if (txt.charAt(i + j) != pat.charAt(j)) { break; } } if (j == M) { return i; } } return N; }

public static int search(String pat, String txt) { int i, M = pat.length(); int j, N = txt.length(); for (i = 0, j = 0; i < N && j < M; i++) { if (txt.charAt(i) == pat.charAt(j)) { j++; } else { i -= j; j = 0; } } if (j == M) { return i - M; } else { return N; } }

int bruteForceSubstringSearch(const std::string& text, const std::string& pattern) { int n = text.length(); int m = pattern.length(); for (int i = 0; i <= n - m; i++) { int j; for (j = 0; j < m; j++) { if (text[i + j] != pattern[j]) { break; } } if (j == m) { return i; } } return -1; }

def brute_force_search(main_string, sub_string): len_main = len(main_string) len_sub = len(sub_string) for i in range(len_main - len_sub + 1): j = 0 while(j < len_sub): if (main_string[i + j] != sub_string[j]): break j += 1 if (j == len_sub): return i return -1

21.3 Knuth-Morris-Pratt

21.3.1 Proposition

Property: KMP substring search accesses no more than chars to search for a pattern of length in a text of length .

Proof: Each pattern char accessed once when constructing DFA; each text char accessed once (in the worst case) when simulating DFA.

Property: KMP constructs dfa[][] in time and space proportional to , where is the alphabet size and is the pattern length.

21.3.2 DFA

Deterministic Finite State Automaton (DFA) is an abstract string-search machine.

Finite number of states (including start and halt).
Exactly one transition for each char in alphabet.
Accept if sequence of transitions lead to halt state.

DFA Construction

If in state (first characters of pattern have already been matched and next char c == pat.charAt(j) (next char matches), go to (now first characters of pattern have been matched).
If in state and next char c != pat.charAt(j), then the last characters of input are pat[1...j - 1], followed by c. Simulate pat[1...j - 1] on DFA and take transition c (only longest possible matched suffix now lies pat[1...j - 1] followed by c).

DFA Construction for Code

Copy dfa[][X] to dfa[][j] for mismatch case.
Set dfa[pat.charAt(j)][j] to for match case.
Update .

public class KMP { private final int R; private final int m; private final int[][] dfa; public KMP(String pat) { this.R = 256; this.m = pat.length(); dfa = new int[R][m]; dfa[pat.charAt(0)][0] = 1; for (int x = 0, j = 1; j < m; j++) { for (int c = 0; c < R; c++) dfa[c][j] = dfa[c][x]; dfa[pat.charAt(j)][j] = j + 1; x = dfa[pat.charAt(j)][x]; } } public KMP(char[] pattern, int R) { this.R = R; this.m = pattern.length; int m = pattern.length; dfa = new int[R][m]; dfa[pattern[0]][0] = 1; for (int x = 0, j = 1; j < m; j++) { for (int c = 0; c < R; c++) dfa[c][j] = dfa[c][x]; dfa[pattern[j]][j] = j + 1; x = dfa[pattern[j]][x]; } } public int search(String txt) { int n = txt.length(); int i, j; for (i = 0, j = 0; i < n && j < m; i++) { j = dfa[txt.charAt(i)][j]; } if (j == m) return i - m; return n; } public int search(char[] text) { int n = text.length; int i, j; for (i = 0, j = 0; i < n && j < m; i++) { j = dfa[text[i]][j]; } if (j == m) return i - m; return n; } }

#include <string> #include <vector> class KMP { private: const int R; const int m; std::vector<std::vector<int>> dfa; public: explicit KMP(const std::string& pat) : R(256), m(static_cast<int>(pat.length())), dfa(R, std::vector<int>(m)) { dfa[pat[0]][0] = 1; for (int x = 0, j = 1; j < m; j++) { for (int c = 0; c < R; c++) dfa[c][j] = dfa[c][x]; dfa[pat[j]][j] = j + 1; x = dfa[pat[j]][x]; } } [[nodiscard]] int search(const std::string& txt) const { const int n = static_cast<int>(txt.length()); int i, j; for (i = 0, j = 0; i < n && j < m; i++) { j = dfa[txt[i]][j]; } if (j == m) return i - m; return n; } };

class KMP: def __init__(self, pat): self.R = 256 self.m = len(pat) self.dfa = [[0] * self.m for _ in range(self.R)] self.dfa[ord(pat[0])][0] = 1 x = 0 for j in range(1, self.m): for c in range(self.R): self.dfa[c][j] = self.dfa[c][x] self.dfa[ord(pat[j])][j] = j + 1 x = self.dfa[ord(pat[j])][x] def search(self, txt): n = len(txt) i, j = 0, 0 while i < n and j < self.m: j = self.dfa[ord(txt[i])][j] i += 1 if j == self.m: return i - self.m return n

21.3.3 NFA

NFA Construction

Use pointer j for comparison.
If i = 0, next[i] = -1.
If pat[i] != pat[j], it means current state j is possible.
If pat[i] == pat[j], it means current state j is impossible, roll back.

public class KMPplus { private final String pattern; private final int[] next; public KMPplus(String pattern) { this.pattern = pattern; int m = pattern.length(); next = new int[m]; int j = -1; for (int i = 0; i < m; i++) { if (i == 0) next[i] = -1; else if (pattern.charAt(i) != pattern.charAt(j)) next[i] = j; else next[i] = next[j]; while (j >= 0 && pattern.charAt(i) != pattern.charAt(j)) { j = next[j]; } j++; } for (int i = 0; i < m; i++) System.out.println("next[" + i + "] = " + next[i]); } public int search(String text) { int m = pattern.length(); int n = text.length(); int i, j; for (i = 0, j = 0; i < n && j < m; i++) { while (j >= 0 && text.charAt(i) != pattern.charAt(j)) j = next[j]; j++; } if (j == m) return i - m; return n; } }

#include <string> #include <vector> #include <iostream> class KMPplus { private: const std::string pattern; std::vector<int> next; public: explicit KMPplus(const std::string& pattern) : pattern(pattern) { const int m = static_cast<int>(pattern.length()); next.resize(m); int j = -1; for (int i = 0; i < m; i++) { if (i == 0) next[i] = -1; else if (pattern[i] != pattern[j]) next[i] = j; else next[i] = next[j]; while (j >= 0 && pattern[i] != pattern[j]) { j = next[j]; } j++; } for (int i = 0; i < m; i++) std::cout << "next[" << i << "] = " << next[i] << std::endl; } [[nodiscard]] int search(const std::string& text) const { const int m = static_cast<int>(pattern.length()); const int n = static_cast<int>(text.length()); int i, j; for (i = 0, j = 0; i < n && j < m; i++) { while (j >= 0 && text[i] != pattern[j]) j = next[j]; j++; } if (j == m) return i - m; return n; } };

class KMPplus: def __init__(self, pattern): self.pattern = pattern m = len(pattern) self.next = [0] * m j = -1 for i in range(m): if i == 0: self.next[i] = -1 elif pattern[i] != pattern[j]: self.next[i] = j else: self.next[i] = self.next[j] while j >= 0 and pattern[i] != pattern[j]: j = self.next[j] j += 1 for i in range(m): print(f"next[{i}] = {self.next[i]}") def search(self, text): m = len(self.pattern) n = len(text) i, j = 0, 0 while i < n and j < m: while j >= 0 and text[i] != self.pattern[j]: j = self.next[j] j += 1 i += 1 if j == m: return i - m return n

21.4 Boyer-Moore

Boyer-Moore

Scan characters in pattern from right to left.
Can skip as many as text chars when finding one not in the pattern.

How much to skip?

Mismatch character not in pattern.
Mismatch character in pattern.
Mismatch character in pattern (but heuristic no help).

Property: Substring search with the Boyer-Moore mismatched character heuristic takes about character (sublinear) compares to search for a pattern of length in a text of length .

Worst Case: Can be as bad as .

Boyer-Moore variant: Can improve worst case to character compares by adding a KMP-like rule to guard against repetitive patterns.

public class BoyerMoore { private final int R; private final int[] right; private char[] pattern; private String pat; public BoyerMoore(String pat) { this.R = 256; this.pat = pat; right = new int[R]; for (int c = 0; c < R; c++) right[c] = -1; for (int j = 0; j < pat.length(); j++) right[pat.charAt(j)] = j; } public BoyerMoore(char[] pattern, int R) { this.R = R; this.pattern = new char[pattern.length]; System.arraycopy(pattern, 0, this.pattern, 0, pattern.length); right = new int[R]; for (int c = 0; c < R; c++) right[c] = -1; for (int j = 0; j < pattern.length; j++) right[pattern[j]] = j; } public int search(String txt) { int M = pat.length(); int N = txt.length(); int skip; for (int i = 0; i <= N - M; i += skip) { skip = 0; for (int j = M - 1; j >= 0; j--) { if (pat.charAt(j) != txt.charAt(i + j)) { skip = Math.max(1, j - right[txt.charAt(i + j)]); break; } } if (skip == 0) return i; } return N; } public int search(char[] text) { int M = pattern.length; int N = text.length; int skip; for (int i = 0; i <= N - M; i += skip) { skip = 0; for (int j = M - 1; j >= 0; j--) { if (pattern[j] != text[i + j]) { skip = Math.max(1, j - right[text[i + j]]); break; } } if (skip == 0) return i; } return N; } }

#include <iostream> #include <string> #include <vector> class BoyerMoore { private: int R; std::vector<int> right; std::string pat; public: explicit BoyerMoore(const std::string& pat) { this->R = 256; this->pat = pat; this->right.resize(R, -1); for (int j = 0; j < pat.size(); j++) { this->right[pat[j]] = j; } } [[nodiscard]] int search(const std::string& txt) const { const int M = static_cast<int>(pat.size()); const int N = static_cast<int>(txt.size()); int skip; for (int i = 0; i <= N - M; i += skip) { skip = 0; for (int j = M - 1; j >= 0; j--) { if (pat[j] != txt[i + j]) { skip = std::max(1, j - right[txt[i + j]]); break; } } if (skip == 0) return i; } return N; } };

class BoyerMoore: def __init__(self, pat): self.R = 256 self.pat = pat self.right = [-1] * self.R for j in range(len(pat)): self.right[ord(pat[j])] = j def search(self, txt): M = len(self.pat) N = len(txt) skip = 1 for i in range(0, N - M + 1, skip): skip = 0 for j in range(M - 1, -1, -1): if self.pat[j] != txt[i + j]: skip = max(1, j - self.right[ord(txt[i + j])]) break if skip == 0: return i return N

21.5 Rabin-Karp

Rabin-Karp (Modular Hashing)

Compute a hash of pattern characters to .
For each , compute a hash of text characters to .
If pattern hash = text substring hash, check for a match.

Modular Hashing Function: Using the notation for txt.charAt(i), we wish to compute:

M-digit, base-R integer, modulo Q.

Based on the function above, we can get:

public class RabinKarp { private final long patHash; private final int M; private final long Q; private final int R; private long RM; public RabinKarp(String pat) { M = pat.length(); R = 256; Q = longRandomPrime(); RM = 1; for (int i = 1; i <= M - 1; i++) RM = (R * RM) % Q; patHash = hash(pat, M); } private long hash(String key, int M) { long h = 0; for (int j = 0; j < M; j++) h = (R * h + key.charAt(j)) % Q; return h; } // Las Vegas version: does pat[] match txt[i..i-M+1] ? private boolean check(String txt, int i) { for (int j = 0; j < M; j++) if (patHash != hash(txt.substring(i, i + M), M)) return false; return true; } // Monte Carlo version: always return true private static long longRandomPrime() { return (1L << 31) - 1; } public int search(String txt) { int N = txt.length(); if (N < M) return N; long txtHash = hash(txt, M); if ((patHash == txtHash) && check(txt, 0)) return 0; for (int i = M; i < N; i++) { txtHash = (txtHash + Q - RM * txt.charAt(i - M) % Q) % Q; txtHash = (txtHash * R + txt.charAt(i)) % Q; int offset = i - M + 1; if ((patHash == txtHash) && check(txt, offset)) return offset; } return N; } }

#include <iostream> #include <string> class RabinKarp { private: long long patHash; int M; long long Q; int R; long long RM; std::string pat; public: explicit RabinKarp(const std::string& pat) : pat(pat) { M = static_cast<int>(pat.length()); R = 256; Q = longRandomPrime(); RM = 1; for (int i = 1; i <= M - 1; i++) RM = (R * RM) % Q; patHash = hash(pat, M); } [[nodiscard]] long long hash(const std::string& key, const int M) const { long long h = 0; for (int j = 0; j < M; j++) h = (R * h + key[j]) % Q; return h; } // Las Vegas version: does pat[] match txt[i..i-M+1] ? [[nodiscard]] bool check(const std::string& txt, const int i) const { for (int j = 0; j < M; j++) if (txt[i + j] != pat[j]) return false; return true; } // Monte Carlo version: always return true static long long longRandomPrime() { return 16777213; } [[nodiscard]] int search(const std::string& txt) const { const int N = static_cast<int>(txt.length()); if (N < M) return N; long long txtHash = hash(txt, M); if ((patHash == txtHash) && check(txt, 0)) return 0; for (int i = M; i < N; i++) { txtHash = (txtHash + Q - RM * txt[i - M] % Q) % Q; txtHash = (txtHash * R + txt[i]) % Q; int offset = i - M + 1; if ((patHash == txtHash) && check(txt, offset)) return offset; } return N; } };

def long_random_prime(): return (1 << 31) - 1 class RabinKarp: def __init__(self, pat): self.pat = pat self.M = len(pat) self.R = 256 self.Q = long_random_prime() self.RM = 1 for i in range(1, self.M): self.RM = (self.R * self.RM) % self.Q self.pat_hash = self.hash(pat, self.M) def hash(self, key, M): h = 0 for j in range(M): h = (self.R * h + ord(key[j])) % self.Q return h # Las Vegas version: does pat[] match txt[i..i-M+1] ? def check(self, txt, i): for j in range(self.M): if self.pat_hash != self.hash(txt[i:i+self.M], self.M): return False return True # Monte Carlo version: always return true def search(self, txt): N = len(txt) if N < self.M: return N txt_hash = self.hash(txt, self.M) if (self.pat_hash == txt_hash) and self.check(txt, 0): return 0 for i in range(self.M, N): txt_hash = (txt_hash + self.Q - self.RM * ord(txt[i - self.M]) % self.Q) % self.Q txt_hash = (txt_hash * self.R + ord(txt[i])) % self.Q offset = i - self.M + 1 if (self.pat_hash == txt_hash) and self.check(txt, offset): return offset return N

Cost of searching for an -character pattern in an -character text

Algorithm	Version	Operation Count		Backup in Input?	Correct?	Extra Space
Algorithm	Version	Guarantee	Typical	Backup in Input?	Correct?	Extra Space
Brute Force	-			yes	yes
Knuth-Morris-Pratt	full DFA			no	yes
Knuth-Morris-Pratt	mismatch transitions only			no	yes
Boyer-Moore	full algorithm			yes	yes
Boyer-Moore	mismatched char heuristic only			yes	yes
Rabin-Karp*	Monte Carlo			no	yes*
Rabin-Karp*	Las Vegas	*		yes	yes

*: probabilisitic guarantee, with uniform hash function

22 Regular Expressions

22.1 Regular Expressions

Pattern Searching: Find one of a specified set of strings in text.

Applications

Genomics: test for certain pattrn of base sequence
Syntax highlighting
Google code search
Scan for virus signatures
Process natural language
Specify a programming language
Access information in digital libraries
Search genome using PROSITE patterns
Filter text (spam, NetNanny, Carnivore, malware)
Validate data-entry fields (dates, email, URL, credit card)
Compile a Java program
Crawl and index the Web
Read in data stored in ad hoc input file format
Create Java documentation from Javadoc comments
...

Regular Expressions: A notation to specify a set of strings.

Operation	Order	Example RE	Matches	Does not Match
Concatenaion	3	AABAAB	AABAAB	every other string
Or	4	AA\|BAAB	AA BAAB	every other string
Closure	2	AB*A	AA ABA ABBA ABBBBBBBBA	AB ABABA
Parenthesis	1	A(A\|B)AAB	AAAAB ABAAB	every other string
Parenthesis	1	(AB)*A	A ABABABABABA	AA ABBA

Shortcuts

Operation	Example RE	Matches	Does not Match
Wildcard	.U.U.U.	CUMULUS JUGULUM	SUCCUBUS TUMULTUOUS
Character Class	[A-Za-z][a-z]*	Word Capitalized	camelCase 4illegal
At Least 1	A(BC)+DE	ABCDE ABCBCDE	ADE BCDE
Exactly k	[0-9]{5}-[0-9]{4}	08540-1321 19072-5541	111111111 166-54-111

Examples

Regular Expression	Matches	Does not Match
.SPB. (substring search)	RASPBERRY CRISPBREAD	SUBSPACE SUBSPECIES
[0-9]{3}-[0-9]{2}-[0-9]{4} (U.S. Social Security numbers)	166-11-4433 166-45-1111	11-55555555 8675309
[a-z]+@([a-z]+\.)+(edu\|com) (simplified email addresses)	wayne@princeton.edu rs@princeton.edu	spam@nowhere
[$_A-Za-z][$_A-Za-z0-9]* (Java identifiers)	ident3 PatternMatcher	3a ident#3

Regular Expression

Matches

Does not Match

.*SPB.*

(substring search)

RASPBERRY

CRISPBREAD

SUBSPACE

SUBSPECIES

[0-9]{3}-[0-9]{2}-[0-9]{4}

(U.S. Social Security numbers)

166-11-4433

166-45-1111

11-55555555

8675309

[a-z]+@([a-z]+\.)+(edu|com)

(simplified email addresses)

wayne@princeton.edu

rs@princeton.edu

spam@nowhere

[$_A-Za-z][$_A-Za-z0-9]*

(Java identifiers)

ident3

PatternMatcher

3a

ident#3

Caveat

Writing a RE is like writing a program.
Need to understand programming model.
Can be easier to write than read.
Can be difficult to debug.

22.2 REs and NFAs

Kleene's theorem

For any DFA, there exists a RE that describes the same set of strings.
For any RE, there exists a DFA that recognizes the same set of strings.

Regular-expression-matching NFA:

We assume RE enclosed in parentheses.
One state per RE character (start = 0, accept = M).
Red ε-transition (change state, but don't scan text).
Black match transition (change state and scan to next text char).
Accept if any sequence of transitions ends in accept state after scanning all text characters.

NFA Simulation

Check whether input matches all possible states that NFA could be.
Reject otherwise.

Construction

Parenthesis: Add ε-transition edge from parentheses to next state.
Closure: Add three ε-transition edges for each * operator.
Or: Add two ε-transition edges for each | operator.

NFA Construction

Left parenthesis
- Add ε-transition to next state.
- Push index of state corresponding to ( onto stack.
Alphabet symbol
- Add match transition to next state.
- Do one-character lookahead: add ε-transition if next character is *.
Or symbol
- Push index of state corresponding to | onto stack.
Right parenthesis
- Add ε-transition to next state.
- Pop correponding ( and possibly intervening |; add ε-transition edges for or.
- Do one-character lookahead: add ε-transition if next character is *.

Property: Determining whether an -character text is recognized by the NFA corresponding to an -character pattern takes time proportional to in the worst case.

Proof: For each of the text characters, we iterate through a set of states of size no more than and run DFS on the graph of ε-transitions.

Property: Building the NFA corresponding to an -character RE takes time and space proportional to .

Proof: For each of the characters in the RE, we add at most three ε-transitions and execute at most two stack operations.

Directed DFS Implementation

import java.util.List; public class DirectedDFS { private final boolean[] marked; private int count; public DirectedDFS(DirectedGraph G, int s) { marked = new boolean[G.getNumVertices()]; validateVertex(s); dfs(G, s); } public DirectedDFS(DirectedGraph G, Iterable<Integer> sources) { marked = new boolean[G.getNumVertices()]; validateVertices(sources); for (int v : sources) { if (!marked[v]) dfs(G, v); } } private void dfs(DirectedGraph G, int v) { count++; marked[v] = true; List<Integer> neighbors = G.getAdjacencyList().get(v); for (int w : neighbors) { if (!marked[w]) dfs(G, w); } } public boolean marked(int v) { validateVertex(v); return marked[v]; } public int count() { return count; } private void validateVertex(int v) { int V = marked.length; if (v < 0 || v >= V) throw new IllegalArgumentException("vertex " + v + " is not between 0 and " + (V - 1)); } private void validateVertices(Iterable<Integer> vertices) { if (vertices == null) { throw new IllegalArgumentException("argument is null"); } int vertexCount = 0; for (Integer v : vertices) { vertexCount++; if (v == null) { throw new IllegalArgumentException("vertex is null"); } validateVertex(v); } if (vertexCount == 0) { throw new IllegalArgumentException("zero vertices"); } } }

class DirectedDFS: def __init__(self, graph, source): self.marked = [False] * graph.get_num_vertices() self.count = 0 if isinstance(source, int): self._dfs(graph, source) elif isinstance(source, list): for s in source: if not self.marked[s]: self._dfs(graph, s) def _dfs(self, graph, v): self.count += 1 self.marked[v] = True for w in graph.adjacency_list[v]: if not self.marked[w]: self._dfs(graph, w) def marked_vertex(self, v): return self.marked[v] def get_count(self): return self.count

NFA Implementation

import java.util.ArrayList; import java.util.List; import java.util.Stack; public class NFA { private final DirectedGraph graph; private final String regexp; private final int m; public NFA(String regexp) { this.regexp = regexp; m = regexp.length(); Stack<Integer> ops = new Stack<>(); graph = new DirectedGraph(m + 1); for (int i = 0; i < m; i++) { int lp = i; if (regexp.charAt(i) == '(' || regexp.charAt(i) == '|') ops.push(i); else if (regexp.charAt(i) == ')') { int or = ops.pop(); if (regexp.charAt(or) == '|') { lp = ops.pop(); graph.addEdge(lp, or + 1); graph.addEdge(or, i); } else if (regexp.charAt(or) == '(') lp = or; else assert false; } if (i < m - 1 && regexp.charAt(i + 1) == '*') { graph.addEdge(lp, i + 1); graph.addEdge(i + 1, lp); } if (regexp.charAt(i) == '(' || regexp.charAt(i) == '*' || regexp.charAt(i) == ')') graph.addEdge(i, i + 1); } if (!ops.isEmpty()) throw new IllegalArgumentException("Invalid regular expression"); } public boolean recognizes(String txt) { DirectedDFS dfs = new DirectedDFS(graph, 0); List<Integer> pc = new ArrayList<>(); for (int v = 0; v < graph.getNumVertices(); v++) if (dfs.marked(v)) pc.add(v); for (int i = 0; i < txt.length(); i++) { if (txt.charAt(i) == '*' || txt.charAt(i) == '|' || txt.charAt(i) == '(' || txt.charAt(i) == ')') throw new IllegalArgumentException("text contains the metacharacter '" + txt.charAt(i) + "'"); List<Integer> match = new ArrayList<>(); for (int v : pc) { if (v == m) continue; if ((regexp.charAt(v) == txt.charAt(i)) || regexp.charAt(v) == '.') match.add(v + 1); } if (match.isEmpty()) continue; dfs = new DirectedDFS(graph, match); pc = new ArrayList<>(); for (int v = 0; v < graph.getNumVertices(); v++) if (dfs.marked(v)) pc.add(v); if (pc.isEmpty()) return false; } for (int v : pc) if (v == m) return true; return false; } }

from DirectedGraph import DirectedGraph from DirectedDFS import DirectedDFS class NFA: def __init__(self, regexp): self.regexp = regexp self.m = len(regexp) self.graph = DirectedGraph(self.m + 1) ops = [] for i in range(self.m): lp = i if regexp[i] == '(' or regexp[i] == '|': ops.append(i) elif regexp[i] == ')': or_op = ops.pop() if regexp[or_op] == '|': lp = ops.pop() self.graph.add_edge(lp, or_op + 1) self.graph.add_edge(or_op, i) elif regexp[or_op] == '(': lp = or_op else: assert False if i < self.m - 1 and regexp[i + 1] == '*': self.graph.add_edge(lp, i + 1) self.graph.add_edge(i + 1, lp) if regexp[i] == '(' or regexp[i] == '*' or regexp[i] == ')': self.graph.add_edge(i, i + 1) if ops: raise ValueError("Invalid regular expression") def recognizes(self, txt): dfs = DirectedDFS(self.graph, 0) pc = [v for v in range(self.graph.get_num_vertices()) if dfs.marked_vertex(v)] for i in range(len(txt)): if txt[i] in ['*', '|', '(', ')']: raise ValueError(f"text contains the metacharacter '{txt[i]}'") match = [] for v in pc: if v == self.m: continue if (self.regexp[v] == txt[i]) or self.regexp[v] == '.': match.append(v + 1) if not match: continue dfs = DirectedDFS(self.graph, match) pc = [v for v in range(self.graph.get_num_vertices()) if dfs.marked_vertex(v)] if not pc: return False for v in pc: if v == self.m: return True return False

23 Data Compression

23.1 Data Compression Introduction

Application

Generic file compression
- Files: GZIP, BZIP, 7z
- Archivers: PKZIP
- File systems: NTFS, HFS+, ZFS
Multimedia
- Images: GIF, JPEG
- Sound: MP3
- Video: MPEG, DivX™, HDTV
Communication
- ITU-T T4 Group 3 Fax
- V.42bis modem
- Skype

23.2 Run-Length Coding

Example

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1

4-bit counts to represent alternating runs of 0s and 1s: 15 0s, then 7 1s, then 7 0s, then 11 1s.

Applications: JPEG, ITU-T T4 Group 3 Fax, ...

import java.io.ByteArrayInputStream; import java.io.ByteArrayOutputStream; import java.io.IOException; public class RunLength { private static final int R = 256; private static final int LG_R = 8; private RunLength() { } public static byte[] expand(byte[] compressed) throws IOException { ByteArrayInputStream in = new ByteArrayInputStream(compressed); ByteArrayOutputStream out = new ByteArrayOutputStream(); boolean b = false; while (in.available() > 0) { int run = in.read(); for (int i = 0; i < run; i++) { out.write(b ? 1 : 0); } b = !b; } return out.toByteArray(); } public static byte[] compress(byte[] input) throws IOException { ByteArrayOutputStream out = new ByteArrayOutputStream(); int run = 0; boolean old = false; for (byte bVal : input) { boolean b = bVal != 0; if (b != old) { out.write(run); run = 1; old = !old; } else { if (run == R - 1) { out.write(run); run = 0; out.write(run); } run++; } } out.write(run); return out.toByteArray(); } private static void printByteArray(byte[] arr) { for (byte b : arr) { System.out.print(b + " "); } System.out.println(); } }

#include <iostream> #include <vector> #include <sstream> constexpr int R = 256; constexpr int LG_R = 8; std::vector<unsigned char> expand(const std::vector<unsigned char>& compressed) { std::vector<unsigned char> expanded; bool b = false; for (const unsigned char run : compressed) { for (int i = 0; i < run; ++i) { expanded.push_back(b ? 1 : 0); } b = !b; } return expanded; } std::vector<unsigned char> compress(const std::vector<unsigned char>& input) { std::vector<unsigned char> compressed; int run = 0; bool old = false; for (const unsigned char bVal : input) { bool b = bVal != 0; if (b != old) { compressed.push_back(run); run = 1; old = !old; } else { if (run == R - 1) { compressed.push_back(run); run = 0; compressed.push_back(run); } run++; } } compressed.push_back(run); return compressed; } void printByteArray(const std::vector<unsigned char>& arr) { for (const unsigned char b : arr) { std::cout << static_cast<int>(b) << " "; } std::cout << std::endl; }

R = 256 LG_R = 8 def expand(compressed): expanded = [] b = False for run in compressed: expanded.extend([1 if b else 0] * run) b = not b return expanded def compress(input_data): compressed = [] run = 0 old = False for b_val in input_data: b = b_val != 0 if b != old: compressed.append(run) run = 1 old = not old else: if run == R - 1: compressed.append(run) run = 0 compressed.append(run) run += 1 compressed.append(run) return compressed def print_byte_array(arr): print(*arr)

23.3 Huffman Coding

Inorder to produce prefix-free code, we need to ensure that no codeword is a prefix of another.

23.4 LZW Coding

24 Reductions

24.1 Introduction

Reduction: Problem reduces to problem if you can use an algorithm that solves to help solve .

Cost of solving = total cost of solving + cost of reduction

Example 1: Finding the median reduces to sorting

Find the media of N items

Sort items.
Return item in the middle.

Cost of solving this problem:

Example 2: Element distinctness reduces to sorting

Element distinctness on N items

Sort items.
Check adjacent pairs for equality.

Cost of solving this problem:

24.2 Designing Algorithms

Examples

3-collinear reduces to sorting.
Finding the median reduces to sorting.
Element distinctness reduces to sorting.
CPM reduces to topological sort.
Arbitrage reduces to shortest paths.
Burrows-Wheeler transform reduces to suffix sort.

For more examples on algorithm designing using reductions, please visit convex hull or shortest path.

24.3 Establishing Lower Bounds

Very difficult to establish lower bounds from scratch => Use reductions

Definition: Problem linear-time reduces to problem if can be solved with:

Linear number of standard computational steps.
Constant number of calls to .

Property: In quadratic decision tree model, any algorithm for sorting integers requires steps. Sorting linear-time reduces to convex hull.

Proof

Sorting instance: , , ...,
Convex hull instance: , , ...,
Region is convex => all points are on hull.
Starting at point with most negative , counterclockwise order of hull points yields integers in ascending order.

24.4 Classifying Problems

Goal: Classyify problem with algorithm that matches lower bound.

Integer Arithmetic Reductions

Goal: Given two -bit integers, compute their product.

Problem	Arithmetic	Order of Growth
Integer Multiplication		M(N)
Integer Division	,	M(N)
Integer Square		M(N)
Integer Square Root		M(N)

Integer arithmetic problems with the same complexity as integer multiplication.

Matrix Multiplication

Goal: Given two -by- matrices, compute their product.

Problem	Linear Algebra	Order of Growth
Matrix Multiplication		MM(N)
Matrix Inversion		MM(N)
Determinant		MM(N)
System of Linear Equations		MM(N)
LU Decomposition		MM(N)
Least Squares	min	MM(N)

Numerical linear algebra problems with the same complexity as matrix multiplication.

25 Linear Programming

Linear Programming: A problem-solving model for optimal allocation of scarce resources.

Agriculture: Diet problem
Computer science: Compiler register allocation, data mining
Electrical engineering: VLSI design, optimal clocking
Energy: Blending petroleum products
Economics: Equilibrium theory, two-person zero-sum games
Environment: Water quality management
Finance: Portfolio optimization
Logistics: Supply-chain management
Management: Hotel yield management
Marketing: Direct mail advertising
Manufacturing: Production line balancing, cutting stock
Medicine: Radioactive seed placement in cancer treatment
Operations research: Airline crew assignment, vehicle routing
Physics: Ground states of 3-D Ising spin glasses
Telecommunication: Network design, Internet routing
Sports: Scheduling ACC basketball, handicapping horse races

25.1 Brewer's Problem

Brewer's problem: Small brewery produces ale and beer, Production limited by scarce resources: corn, hops, barley malt.

$13 profit per barrel: 5 pounds corn, 4 ounces hops, 35 pounds halt
$23 profit per barrel: 15 pounds corn, 4 ounces hops, 20 pounds halt
corn: 480 lbs
hops: 160 oz
malt: 1190 lbs

Now we have:

Standard form linear program

Extreme point property: If there exists an optimal solution to (P), then there exists one that is an extreme point.

25.2 Simplex Algorithm

Generic Implementation

Start at some extreme point.
Pivot from one extreme point to an adjacent one.
Repeat until optimal.

Simplex Algorithm - Basic feasible solution

Set nonbasic variables to , solve for remaining variables.
Solve equations in unknowns.
If unique and feasible => BFS
BFS <=> extreme point

Recall:

Initialize:

Pivot 1: Substitute y = (1/15) (480 – 5A – SC) and add y into the basis (rewrite 2nd equation, eliminate y in 1st, 3rd, and 4th equations)
basis = {B, S _H, S _M}
x = S _C = 0
Z = 736
y = 32
S _H = 32
S _M = 550
Q: Why pivot on column 2 (corresponding to variable )?
A:
- Its objective function coefficient is positive.
  (each unit increase in from increases objective value by )
- Pivoting on column 1 (corresponding to ) also OK.
Q: Why pivot on row 2?
A:
- Preserves feasibility by ensuring .
- Minimum ratio rule: min {, , }
Pivot 2: Substitute x = (3/8) (32 + (4/15) S _C – S _H) and add x into the basis (rewrite 3rd equation, eliminate x in 1st, 2nd, and 4th equations)
basis = { x, y, S _M}
S _C = S _H = 0
Z = 800
y = 28
x = 12
S _M = 110
Q: When to stop pivoting?
A: When no objective function coefficient is positive.
Q: Why is resulting solution optimal?
A: Any feasible solution satisfies current system of equations.
- In particular:
- Thus, optimal objective value since .
- Current BFS has value => optimal.

25.3 Implementation

Bland's Rule: Find entering column using Bland's rule: index of first column whose objective function coefficient is positive.

Min-ratio Rule: Find leaving row using min ratio rule. (Bland's rule: if a tie, choose first such row)

Pivot: Pivot on element row , column .

Property: In typical practical applications, simplex algorithm terminates after at most pivots.

public class LinearProgramming { private static final double EPSILON = 1.0E-10; private final double[][] a; private final int m; private final int n; private final int[] basis; public LinearProgramming(double[][] A, double[] b, double[] c) { m = b.length; n = c.length; for (int i = 0; i < m; i++) { if (!(b[i] >= 0)) { throw new IllegalArgumentException("RHS must be nonnegative"); } } a = new double[m + 1][n + m + 1]; for (int i = 0; i < m; i++) { System.arraycopy(A[i], 0, a[i], 0, n); } for (int i = 0; i < m; i++) { a[i][n + i] = 1.0; } System.arraycopy(c, 0, a[m], 0, n); for (int i = 0; i < m; i++) { a[i][m + n] = b[i]; } basis = new int[m]; for (int i = 0; i < m; i++) { basis[i] = n + i; } solve(); assert check(A, b, c); } private void solve() { while (true) { int q = bland(); if (q == -1) { break; } int p = minRatioRule(q); if (p == -1) { throw new ArithmeticException("Linear program is unbounded"); } pivot(p, q); basis[p] = q; } } private int bland() { for (int j = 0; j < m + n; j++) { if (a[m][j] > 0) { return j; } } return -1; } private int minRatioRule(int q) { int p = -1; for (int i = 0; i < m; i++) { if (a[i][q] <= EPSILON) { continue; } else if (p == -1) { p = i; } else if ((a[i][m + n] / a[i][q]) < (a[p][m + n] / a[p][q])) { p = i; } } return p; } private void pivot(int p, int q) { for (int i = 0; i <= m; i++) { for (int j = 0; j <= m + n; j++) { if (i != p && j != q) { a[i][j] -= a[p][j] * (a[i][q] / a[p][q]); } } } for (int i = 0; i <= m; i++) { if (i != p) { a[i][q] = 0.0; } } for (int j = 0; j <= m + n; j++) { if (j != q) { a[p][j] /= a[p][q]; } } a[p][q] = 1.0; } public double value() { return -a[m][m + n]; } public double[] primal() { double[] x = new double[n]; for (int i = 0; i < m; i++) { if (basis[i] < n) { x[basis[i]] = a[i][m + n]; } } return x; } public double[] dual() { double[] y = new double[m]; for (int i = 0; i < m; i++) { y[i] = -a[m][n + i]; if (y[i] == -0.0) { y[i] = 0.0; } } return y; } private boolean isPrimalFeasible(double[][] A, double[] b) { double[] x = primal(); for (int j = 0; j < x.length; j++) { if (x[j] < -EPSILON) { System.out.println("x[" + j + "] = " + x[j] + " is negative"); return false; } } for (int i = 0; i < m; i++) { double sum = 0.0; for (int j = 0; j < n; j++) { sum += A[i][j] * x[j]; } if (sum > b[i] + EPSILON) { System.out.println("not primal feasible"); System.out.println("b[" + i + "] = " + b[i] + ", sum = " + sum); return false; } } return true; } private boolean isDualFeasible(double[][] A, double[] c) { double[] y = dual(); for (int i = 0; i < y.length; i++) { if (y[i] < -EPSILON) { System.out.println("y[" + i + "] = " + y[i] + " is negative"); return false; } } for (int j = 0; j < n; j++) { double sum = 0.0; for (int i = 0; i < m; i++) { sum += A[i][j] * y[i]; } if (sum < c[j] - EPSILON) { System.out.println("not dual feasible"); System.out.println("c[" + j + "] = " + c[j] + ", sum = " + sum); return false; } } return true; } private boolean isOptimal(double[] b, double[] c) { double[] x = primal(); double[] y = dual(); double value = value(); double value1 = 0.0; for (int j = 0; j < x.length; j++) { value1 += c[j] * x[j]; } double value2 = 0.0; for (int i = 0; i < y.length; i++) { value2 += y[i] * b[i]; } if (Math.abs(value - value1) > EPSILON || Math.abs(value - value2) > EPSILON) { System.out.println("value = " + value + ", cx = " + value1 + ", yb = " + value2); return false; } return true; } private boolean check(double[][] A, double[] b, double[] c) { return isPrimalFeasible(A, b) && isDualFeasible(A, c) && isOptimal(b, c); } private static void test(double[][] A, double[] b, double[] c) { LinearProgramming lp; try { lp = new LinearProgramming(A, b, c); } catch (ArithmeticException e) { System.out.println(e); return; } System.out.println("value = " + lp.value()); double[] x = lp.primal(); for (int i = 0; i < x.length; i++) { System.out.println("x[" + i + "] = " + x[i]); } double[] y = lp.dual(); for (int j = 0; j < y.length; j++) { System.out.println("y[" + j + "] = " + y[j]); } } }

#include <cassert> #include <iostream> #include <vector> #include <cmath> #include <random> class LinearProgramming { private: static constexpr double EPSILON = 1.0E-10; std::vector<std::vector<double>> a; int m; int n; std::vector<int> basis; public: LinearProgramming(const std::vector<std::vector<double>>& A, const std::vector<double>&b, const std::vector<double>& c) { m = static_cast<int>(b.size()); n = static_cast<int>(c.size()); for (int i = 0; i < m; i++) { if (!(b[i] >= 0)) { throw std::runtime_error("RHS must be nonnegative"); } } a.resize(m + 1, std::vector<double>(n + m + 1)); for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { a[i][j] = A[i][j]; } } for (int i = 0; i < m; i++) { a[i][n + i] = 1.0; } for (int j = 0; j < n; j++) { a[m][j] = c[j]; } for (int i = 0; i < m; i++) { a[i][m + n] = b[i]; } basis.resize(m); for (int i = 0; i < m; i++) { basis[i] = n + i; } solve(); assert(check(A, b, c)); } private: void solve() { while (true) { const int q = bland(); if (q == -1) { break; } int p = minRatioRule(q); if (p == -1) { throw std::runtime_error("Linear program is unbounded"); } pivot(p, q); basis[p] = q; } } [[nodiscard]] int bland() const { for (int j = 0; j < m + n; j++) { if (a[m][j] > 0) { return j; } } return -1; } [[nodiscard]] int minRatioRule(const int q) const { int p = -1; for (int i = 0; i < m; i++) { if (a[i][q] <= EPSILON) { continue; } else if (p == -1) { p = i; } else if ((a[i][m + n] / a[i][q]) < (a[p][m + n] / a[p][q])) { p = i; } } return p; } void pivot(const int p, const int q) { for (int i = 0; i <= m; i++) { for (int j = 0; j <= m + n; j++) { if (i != p && j != q) { a[i][j] -= a[p][j] * (a[i][q] / a[p][q]); } } } for (int i = 0; i <= m; i++) { if (i != p) { a[i][q] = 0.0; } } for (int j = 0; j <= m + n; j++) { if (j != q) { a[p][j] /= a[p][q]; } } a[p][q] = 1.0; } public: [[nodiscard]] double value() const { return -a[m][m + n]; } [[nodiscard]] std::vector<double> primal() const { std::vector<double> x(n); for (int i = 0; i < m; i++) { if (basis[i] < n) { x[basis[i]] = a[i][m + n]; } } return x; } [[nodiscard]] std::vector<double> dual() const { std::vector<double> y(m); for (int i = 0; i < m; i++) { y[i] = -a[m][n + i]; if (y[i] == -0.0) { y[i] = 0.0; } } return y; } private: [[nodiscard]] bool isPrimalFeasible(const std::vector<std::vector<double>>& A, const std::vector<double>& b) const { const std::vector<double> x = primal(); for (int j = 0; j < x.size(); j++) { if (x[j] < -EPSILON) { std::c false; } } for (int i = 0; i < m; i++) { double sum = 0.0; for (int j = 0; j < n; j++) { sum += A[i][j] * x[j]; } if (sum > b[i] + EPSILON) { std::cout << "not primal feasible" << std::endl; std::cout << "b[" << i << "] = " << b[i] << ", sum = " << sum << std::endl; return false; } } return true; } [[nodiscard]] bool isDualFeasible(const std::vector<std::vector<double>>& A, const std::vector<double>& c) const { std::vector<double> y = dual(); for (size_t i = 0; i < y.size(); i++) { if (y[i] < -EPSILON) { std::cout << "y[" << i << "] = " << y[i] << " is negative" << std::endl; return false; } } for (int j = 0; j < n; j++) { double sum = 0.0; for (int i = 0; i < m; i++) { sum += A[i][j] * y[i]; } if (sum < c[j] - EPSILON) { std::cout << "not dual feasible" << std::endl; std::cout << "c[" << j << "] = " << c[j] << ", sum = " << sum << std::endl; return false; } } return true; } [[nodiscard]] bool isOptimal(const std::vector<double>& b, const std::vector<double>&c) const { const std::vector<double> x = primal(); const std::vector<double> y = dual(); const double value = this->value(); double value1 = 0.0; for (size_t j = 0; j < x.size(); j++) value1 += c[j] * x[j]; double value2 = 0.0; for (size_t i = 0; i < y.size(); i++) value2 += y[i] * b[i]; if (std::abs(value - value1) > EPSILON || std::abs(value - value2) > EPSILON) { std::cout << "value = " << value << ", cx = " << value1 << ", yb = " << value2 << std::endl; return false; } return true; } [[nodiscard]] bool check(const std::vector<std::vector<double>>& A, const std::vector<double>& b, const std::vector<double>& c) const { return isPrimalFeasible(A, b) && isDualFeasible(A, c) && isOptimal(b, c); } }; void test(const std::vector<std::vector<double>>& A, const std::vector<double>&b, const std::vector<double>& c) { try { const LinearProgramming lp(A, b, c); std::cout << "value = " << lp.value() << std::endl; const std::vector<double> x = lp.primal(); for (size_t i = 0; i < x.size(); i++) std::cout << "x[" << i << "] = " << x[i] << std::endl; const std::vector<double> y = lp.dual(); for (size_t j = 0; j < y.size(); j++) std::cout << "y[" << j << "] = " << y[j] << std::endl; } catch (const std::runtime_error& error) { std::cerr << error.what() << std::endl; } }

import numpy as np EPSILON = 1e-10 class LinearProgramming: def __init__(self, A, b, c): self.m = len(b) self.n = len(c) if not all(b_i >= 0 for b_i in b): raise ValueError("RHS must be nonnegative") self.a = np.zeros((self.m + 1, self.n + self.m + 1)) self.a[:self.m, :self.n] = A np.fill_diagonal(self.a[:self.m, self.n:self.n + self.m], 1.0) self.a[self.m, :self.n] = c self.a[:self.m, self.n + self.m] = b self.basis = list(range(self.n, self.n + self.m)) self.solve() assert self.check(A, b, c) def solve(self): while True: q = self.bland() if q == -1: break p = self.min_ratio_rule(q) if p == -1: raise ArithmeticError("Linear program is unbounded") self.pivot(p, q) self.basis[p] = q def bland(self): for j in range(self.m + self.n): if self.a[self.m, j] > 0: return j return -1 def min_ratio_rule(self, q): p = -1 min_ratio = float('inf') for i in range(self.m): if self.a[i, q] > EPSILON: ratio = self.a[i, self.n + self.m] / self.a[i, q] if ratio < min_ratio: min_ratio = ratio p = i return p def pivot(self, p, q): self.a[p] /= self.a[p, q] for i in range(self.m + 1): if i != p: self.a[i] -= self.a[i, q] * self.a[p] def value(self): return -self.a[self.m, self.n + self.m] def primal(self): x = np.zeros(self.n) for i in range(self.m): if self.basis[i] < self.n: x[self.basis[i]] = self.a[i, self.n + self.m] return x def dual(self): y = -self.a[self.m, self.n:self.n + self.m] return y def is_primal_feasible(self, A, b): x = self.primal() if any(x_j < -EPSILON for x_j in x): print("Primal infeasible: x contains negative values.") return False Ax = A @ x if any(Ax_i > b_i + EPSILON for Ax_i, b_i in zip(Ax, b)): print("Primal infeasible: Ax > b") return False return True def is_dual_feasible(self, A, c): y = self.dual() if any(y_i < -EPSILON for y_i in y): print("Dual infeasible: y contains negative values.") return False yA = y @ A if any(yA_j < c_j - EPSILON for yA_j, c_j in zip(yA, c)): print("Dual infeasible: yA < c") return False return True def is_optimal(self, b, c): x = self.primal() y = self.dual() value = self.value() value1 = c @ x value2 = y @ b if abs(value - value1) > EPSILON or abs(value - value2) > EPSILON: print(f"Not optimal: value = {value}, cx = {value1}, yb = {value2}") return False return True def check(self, A, b, c): return (self.is_primal_feasible(A, b) and self.is_dual_feasible(A, c) and self.is_optimal(b, c)) def test(A, b, c): try: lp = LinearProgramming(A, b, c) print("value =", lp.value()) print("x =", lp.primal()) print("y =", lp.dual()) except ArithmeticError as e: print(e)

30 Catalan Number

30.1 Properties and Formulas

30.2 Applications

It is the number of expressions containing pairs of parentheses which are correctly matched.
For , for example:
((())), (()()), (())(), ()(()), ()()().
It is the number of different ways factors can be completely parenthesized (or the number of ways of associating applications of a binary operator, as in the matrix chain multiplication problem).
For , for example:
((ab)c)d, (a(bc))d, (ab)(cd), a((bc)d), a(b(cd)).
It is the number of full binary trees with leaves , or, equivalently, with a total of internal nodes.
A full binary tree is a tree in which every node has either 0 or 2 children.
For, for example:
It is the number of structurally unique BSTs (binary search trees) which has exactly nodes of unique values from to .
For, for example:
It is the number of Dyck words of length . A Dyck word is a string consisting of X's and Y's such that no initial segment of the string has more Y's than X's.
For example, Dyck words for :
XXXYYY XYXXYY XYXYXY XXYYXY XXYXYY
It is the number of monotonic lattice paths along the edges of a grid with square cells, which do not pass above the diagonal.
- A monotonic path is one which starts in the lower left corner, finishes in the upper right corner, and consists entirely of edges pointing rightwards or upwards.
- Counting such paths is equivalent to counting Dyck words: X stands for "move right" and Y stands for "move up".
A convex polygon with sides can be cut into triangles by connecting vertices with non-crossing line segments (a form of polygon triangulation). The number of triangles formed is and it is the number of different ways that this can be achieved.

30.3 Implementation

public static BigInteger catalan(int n) { BigInteger res = BigInteger.ONE; for (int i = 0; i < n; i++) { res = res.multiply(BigInteger.valueOf(2L * n - i)); res = res.divide(BigInteger.valueOf(i + 1)); } return res.divide(BigInteger.valueOf(n + 1)); }

unsigned long int binomialCoeff(unsigned int n, unsigned int k) { if (k > n) return 0; if (k == 0 || k == n) return 1; unsigned long int res = 1; for (int i = 0; i < k; i++) { res *= (n - i); res /= (i + 1); } return res; } unsigned long int catalan(unsigned int n) { unsigned long int c = binomialCoeff(2 * n, n); return c / (n+1); }

import math def catalan_number(n): return math.comb(2 * n, n) // (n + 1)