Polynomial
Polynomial
「多項式」。離散與連續的橋樑。
f(x) = 2x³ + 3x² - 0x¹ + 1x⁰
多項式視作數值：底數是x，位數是多項式係數。
多項式視作函數：輸入是x，輸出是多項式求值。
多項式可能無限長：有限項多項式、無限項多項式。
多項式必須搭配數域：整數多項式、有理數多項式、……。
Polynomial資料結構
Dense Polynomial：array，一格存一項。索引值是次方，內容是係數。
₀ ₁ ₂ ₃ ₄ ₅ ₆ ₇ ₈ ₉ ___________________ |8|4|0|2| | | | | | | ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
Sparse Polynomial：list或array，一份存一項。內容是係數和次方。係數為零的項，統統去掉，串成一串。
___ ___ ___ ⤷|2|3|→|4|1|→|8|0| ‾‾‾ ‾‾‾ ‾‾‾
Polynomial運算
加、減、乘、除、積、微、展開、分解、代入（求值）、求根（求解）。
addition (x² + 3x + 2) + (x + 2) = x² + 4x + 4 subtraction (x² + 3x + 2) - (x + 2) = x² + 2x multiplication (x² + 3x + 2) × (x + 2) = 2x³ + 5x² + 8x + 4 division (x² + 3x + 2) ÷ (x + 2) = x + 1 integration ∫ (x² + 3x + 2) dx = ¹⁄₃x³ + ³⁄₂x² + 2x + c differentiation d/dx (x² + 3x + 2) = 2x + 3 expansion (x + 1)² = x² + 2x + 1 factorization x² + 2x + 1 = (x + 1)² evaluation (x² + 3x + 2)|x=1 = 6 resolution (x² + 3x + 2) = 0 ⇒ x = -1 or -2
專著《Algorithmes Efficaces en Calcul Formel》。
專著《Algorithms for Computer Algebra》。
三本專著詳細介紹了多項式運算。以下章節不再重複整理，只做簡單介紹。
UVa 392 126 10719 10586 10951 463 930 10428 498 10268 10105
快速多項式運算
多項式運算，普通演算法是O(N²)。幸運的是，多項式乘法，運用複數系統、餘數系統，時間複雜度降為O(NlogN)，接近線性時間。以多項式乘法為根基，各種多項式運算皆可降至接近線性時間。多項式運算從此以後變得快速。
拓展數系的三個元件：對稱（自然數到整數）、數對（整數到有理數）、循環（有理數到複數）。多項式乘法，同時運用三者，得到傅立葉矩陣，形成高速演算法，O(N²)變O(NlogN)。箇中玄妙之處，令人拍案叫絕。
整數運算→多項式運算
整數可以視作多項式。有些整數運算可以套用多項式運算。
比方來說，整數乘法套用多項式乘法。因為3×37 = 111，所以3(3x + 7) = 9x + 21，底數x = 10。
記得認真考慮整數運算的時間複雜度。快速傅立葉轉換內含整數乘法，遞迴使用整數乘法，時間複雜度成為O(N logN loglogN)。
順帶一提，整數乘法的時間複雜度下限仍是懸案。2019年出現新演算法，時間複雜度降至O(NlogN)，很可能就是下限了。
https://en.wikipedia.org/wiki/Fürer's_algorithm
luogu P1601 P1919 P5432
整數運算↛多項式運算
有些整數運算則不可以套用多項式運算。
比方來說，整數分解套用多項式分解，雖然111 = 3×37，但是x² + x + 1 ≠ 3(3x + 7)。我們不知道底數x應該是多少，也不知道是否有減法。
順帶一提，整數分解的時間複雜度類別仍是懸案，目前還不確定它究竟是P問題或是NP-complete問題，相當特別。
Polynomial (Data Structure)
data
read / write
截斷版本。
scan / print
luogu P1067
實作細節
內建資料結構。vector自帶各種建構子、解構子，不用擔心生成、消滅、拷貝、指派導致的記憶體漏失。
建構子。只宣告、不定義。另行賦值。
左值函式。調整容量無法採用左值函式，這是C++設計哲學。調整度數採用左值函式，增加可讀性。
運算子重載。簡化程式碼，增加可讀性。順帶一提，operator[]無法放在struct外面，這是C++設計哲學。
實作細節
計算機當中，項數優於度數。項數的優點：陣列容量不必加一、迴圈判斷式不必寫等於。度數的優點：乘法是度數相加、除法是度數相減。雖然我喜歡項數，但是本篇文章採用度數。
不自動調整度數、容量。不消滅首項的零（減法），不消滅超過度數的數值（截斷）。優點：節省時間。缺點：使用者必須自行調整度數、容量，程式容易生蟲。
函式回傳值納入參數。不新增區域變數。優點：節省空間。缺點：使用者必須自行預知度數、容量，程式容易生蟲。
函式變數改成靜態。只新增一次區域變數。優點：節省時間。缺點：使用者必須自行限制度數、容量，程式容易生蟲。
同一函式重複使用記憶體。使用別名，增加可讀性。
不同函式重複使用全域記憶體。試著想像一下，不同函式的迴圈索引值，嘗試改成全域變數。人工追蹤函式呼叫軌跡，容易出紕漏，缺乏可維護性。因此本篇文章改為採用靜態區域變數。
參數邊界防禦。優點：使用者少寫許多程式碼。
參數相同防禦。缺點：度數必須挪到最後修改。
順帶一提，參數相同，可以使用inplace演算法，重複利用記憶體。優點：得以inplace的情況下，呼叫者不必自行籌備記憶體，少寫許多程式碼。缺點：無法inplace的情況下，被呼叫者必須額外建立記憶體，代價轉嫁。
命名。自訂函式跟內建函式可以撞名，function overloading。區域變數跟全域變數可以撞名，operator::。區域變數跟成員變數可以撞名，this。
Polynomial (Arithmetic)
addition +
遞增法。對應項相加。O(N)。
multiplication ×
遞增法。兩兩相乘、兩層迴圈。O(N²)。
分治法。補零、循環乘法。O(NlogN)。簡寫成O(M(N))。
加速技巧。自己乘自己的情況下，省下一次傅立葉轉換。
加速技巧。計算ab，原本需要三次傅立葉轉換。(a + b𝑖)² = (a² - b²) + 2ab𝑖。a放實部，b放虛部，平方之後，虛部恰是2ab。計算(a + b𝑖)²然後擷取虛部，只需要兩次傅立葉轉換。
「循環多項式Circular Polynomial」：xᵃ = xᵃ⁺ⁿ。
luogu P3803
multiplicative inverse: reciprocal ¹⁄ₓ
遞推法。長除法。O(N²)。
遞推法。Hensel Lifting＝位數倍增＋牛頓法。O(M(N))。
「截斷多項式Truncated Polynomial」：低次項N項。模xᴺ。
http://picks.logdown.com/archives
given f(a) ≡ 0 (mod xⁿ) and aꜛ ≡ a (mod xⁿ) find f(aꜛ) ≡ 0 (mod x²ⁿ) f(aꜛ) ≡ f(a)(aꜛ-a)⁰ + f′(a)(aꜛ-a)¹ + f″(a)(aꜛ-a)²/2! + ... (mod x²ⁿ) f(aꜛ) ≡ f(a) + f′(a)(aꜛ-a) (mod x²ⁿ) ^^^^^^^ aꜛ ≡ a - f(a) / f′(a) (mod x²ⁿ) double the zeros
ab ≡ 1 (mod xⁿ) b ≡ 1/a (mod xⁿ) a ≡ 1/b (mod xⁿ) a - 1/b ≡ 0 (mod xⁿ) let f(b) = a - 1/b , f′(b) ≡ 1/b² bꜛ ≡ b - f(b) / f′(b) (mod x²ⁿ) bꜛ ≡ b - (a - 1/b) / (1/b²) (mod x²ⁿ) bꜛ ≡ b - (a - 1/b) b² (mod x²ⁿ) bꜛ ≡ b - a b² + b (mod x²ⁿ) bꜛ ≡ b (2 - a b) (mod x²ⁿ)
b (mod x¹) := 1/a (mod x¹)
b (mod x²) := (2 - ab) b (mod x¹)
b (mod x⁴) := (2 - ab) b (mod x²)
b (mod x⁸) := (2 - ab) b (mod x⁴)
: :
luogu P4238 P4721
euclidean division ÷
遞推法。長除法。O(N²)。
頭尾顛倒，引入餘式系統，除法改成乘以倒數。O(M(N))。
「餘式Polynomial Residue」：餘數運算推廣成餘式運算。截斷多項式可以轉換到餘式系統進行處理：令模式是xⁿ。
https://www.csd.uwo.ca/~mmorenom/CS424/Lectures/FastDivisionAndGcd.html/node3.html
a = bq + r rev(a) = rev(bq + r) rev(a) = rev(b) rev(q) + rev(r) xⁿ where n = deg(a)-deg(b)+1 rev(a) ≡ rev(b) rev(q) (mod xⁿ) rev(a) / rev(b) ≡ rev(q) (mod xⁿ)
luogu P4512
Polynomial (Function)
evaluation ()
遞增法。逐項次方、逐項相加。O(NlogN)。
遞推法。Horner's Rule。O(N)。
composition ∘
遞推法。Horner's Rule。O(M(AB)A)。
分治法。A對半分。O(M(AB)logA)。
加速技巧。過程都在頻域。O(M(AB))。
順帶一提，進位轉換可以套用composition。
a₀x⁰ + a₁x¹ + a₂x² + a₃x³ = (a₀x⁰ + a₁x¹) + (a₂x⁰ + a₃x¹)x² where x = b₀x⁰ + b₁x¹ + b₂x²
令計算結果是截斷多項式。
遞推法。隨時模xᴺ。O(M(C)A)。
分治法。A對半分。O(M(C)(AB/C)(1+log(min(A,ceil(C/B))))。簡單寫成O(M(C)A)或O(M(C)BlogA)，後者需C≥A。
加速技巧。Brent-Kung Algorithm。B分成低次項k項、高次項B-k項。泰勒展開，只取低階C/k項。先求最低階，再以遞推法求得各階導數、冪數、分母。O(M(C)(klogA + C/k))。令klogA = C/k讓時間複雜度達到最低。O(M(C)sqrt(ClogA))。
把B變小（長度k），分治法開頭變快（長度到達C之前），代價是泰勒展開（長度C/k）。令k = sqrt(C/logA)讓兩者平衡。
https://www.luogu.com.cn/blog/yurzhang/solution-p5373
luogu P5373
composite inverse: inverse function ⁻¹
我沒有研究。
https://en.wikipedia.org/wiki/Lagrange_inversion_theorem
https://math.stackexchange.com/questions/210720/
luogu P5809
decomposition
我沒有研究。
Kozen-Landau Algorithm。O(N²)。
Polynomial (Calculus)
differentiation d/dx
逐項微分。O(N)。
a = a₀x⁰ + a₁x¹ + a₂x² + a₃x³ d/dx a = 1 a₁x¹⁻¹ + 2 a₂x²⁻¹ + 3 a₃x³⁻¹
integration ∫ dx
逐項積分。O(N)。
a = a₀x⁰ + a₁x¹ + a₂x² + a₃x³ ∫a dx = ¹⁄₁ a₀x⁰⁺¹ + ¹⁄₂ a₁x¹⁺¹ + ¹⁄₃ a₂x²⁺¹ + ¹⁄₄ a₃x³⁺¹
natural logarithm ln()
方便起見，討論截斷多項式。
公式解。當輸入的常數項為1，則輸出的常數項為0——得以使用微積分求解。當輸入的常數項不為1，則提出倍率使之為1——倍率是常數項。整數多項式不適用此方法。O(M(N))。
when a > 0 ln(a) = (a-1)¹/1 - (a-1)²/2 + (a-1)³/3 - ...
when a₀ = 1 and a > 0 ln(a) = b (if a₀ = 1 then b₀ = 0) ln′(a) a′ = b′ (1/a) a′ = b′ ∫((1/a) a′) = b when a₀ ≠ 1 and a > 0 ln(a) = ln((a/a₀) a₀) = ln(a/a₀) + ln(a₀)
luogu P4725
natural exponentiation exp()
遞推法。Hensel Lifting。O(M(N))。
b ≡ exp(a) (mod xⁿ) ln(b) ≡ a (mod xⁿ) ln(b) - a ≡ 0 (mod xⁿ) let f(b) = ln(b) - a , f′(b) ≡ 1/b bꜛ ≡ b - f(b) / f′(b) (mod x²ⁿ) bꜛ ≡ b - (ln(b) - a) / (1/b) (mod x²ⁿ) bꜛ ≡ b - (ln(b) - a) b (mod x²ⁿ) bꜛ ≡ b (1 - ln(b) + a) (mod x²ⁿ)
luogu P4726 P6613
exponentiation pow()
分治法。次方對半分。O(M(N)logK)。
公式解。O(M(N))。
b = aᵏ ln(b) = k ln(a) b = exp(k ln(a))
luogu P5245 P5273
square root sqrt()
遞推法。Hensel Lifting。O(M(N))。
公式解。O(M(N))。
b ≡ sqrt(a) (mod xⁿ) b² ≡ a (mod xⁿ) b² - a ≡ 0 (mod xⁿ) let f(b) = b² - a , f′(b) ≡ 2b bꜛ ≡ b - f(b) / f′(b) (mod x²ⁿ) bꜛ ≡ b - (b² - a) / (2b) (mod x²ⁿ) bꜛ ≡ b - b/2 + a/b/2 (mod x²ⁿ) bꜛ ≡ b/2 + a/b/2 (mod x²ⁿ) bꜛ ≡ (b + a/b) / 2 (mod x²ⁿ)
b = a½ ln(b) = ½ ln(a) b = exp(½ ln(a))
luogu P5205 P5277
trigonometric functions sin() cos() tan()
公式解。O(M(N))。
https://en.wikipedia.org/wiki/Trigonometric_functions#Relationship_to_exponential_function_(Euler's_formula)
exp(𝑖x) = cos(x) + 𝑖 sin(x) —→ cos(x) = (exp(𝑖x) + exp(-𝑖x)) / 2 exp(-𝑖x) = cos(x) - 𝑖 sin(x) sin(x) = (exp(𝑖x) - exp(-𝑖x)) / 2𝑖
https://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Derivatives_of_inverse_trigonometric_functions
arcsin(x) = ∫ (x′ / sqrt(1 - x²)) dx arccos(x) = -∫ (x′ / sqrt(1 - x²)) dx arctan(x) = ∫ (x′ / (1 + x²)) dx
luogu P5264 P5265
finite difference Δ
重複K次。O(NK)。
乘以(1-x)ᵏ。次方運算。O(M(N)logK)。
乘以(1-x)ᵏ。自然對數指數，避開次方運算。O(M(N))。
乘以(1-x)ᵏ。二項式展開，遞推法求係數。O(M(N))。
https://www.luogu.com.cn/blog/Soulist/solution-p5488
finite difference: | cumulative sum: b = a (1-x)ᵏ | b = a / (1-x)ᵏ b = a exp(k ln(1-x)) | b = a / exp(k ln(1-x))
luogu P5488
Polynomial (Interpolation)
求值
求值與內插、係值轉換、Unisolvence Theorem唯一解定理。
多項式運算：加減乘除乘方開方；簡易演算法：先多點求值、再內插。時間複雜度稍高，但是容易理解。
多項式循環乘除（循環卷積反卷積）亦然。多點求值是單位根，內插是快速傅立葉轉換。
multipoint evaluation
遞推法。Horner's Rule。O(NK)。
遞增法。求值a(xᵢ)化作求餘數a mod (x-xᵢ)。O(M(N)K)。
分治法。項數對半分、形成二元樹。O(M(N) + M(K)logK)。
一、一次多項式連乘(x-x₁)(x-x₂)...(x-xK)。
二、求餘數，一回合化作logK回合。第一回合，模數是全部連乘(x-x₁)(x-x₂)...(x-xK)。第二回合，模數是前半部、後半部，形成二元樹。最後回合，模數是(x-xᵢ)。
三、加速技巧Tellegen's Principle。
https://cs.stackexchange.com/questions/60239/
https://www.luogu.com.cn/blog/EntropyIncreaser/solution-p5050
luogu P5050
interpolation
基礎知識請見本站文件「Polynomial Interpolation」。
遞增法。拉格朗日內插。O(M(N)N)。
一、一次多項式連乘M(x₁,x₂,...) = (x-x₁)(x-x₂)...(x-xK)。
二、內插P(x₁,x₂,...)化作求商數M(x₁,x₂,...) / (x-xᵢ)。
分治法。項數對半分、形成二元樹。O(M(N)logN)。
一、一次多項式連乘：M(all) = M(low)M(high)。
二、係數：L'Hôpital's Rule。
三、內插：P(all) = P(low)M(high) + P(high)M(low)。
luogu P5158 P4781
monic linear congruences f(x) ≡ rᵢ(x) (mod mᵢ(x))
分治法。式子對半分。O(M(N) + M(K)logK)。
中國餘數定理本質就是內插。
Polynomial (Root Finding)（Under Construction!）
求根
求根與乘法、係根轉換、Fundamental Theorem of Algebra代數基本定理。
factorization
回溯法。Kronecker's Algorithm。兩個多項式整除，所有對應點也會整除。原式的函數點的某個因數，是因式的函數點。遞迴窮舉因式的函數點，以「Newton Interpolation」將函數點變成多項式；試除原式，如果整除，就是因式。時間複雜度是指數時間。
「不可約多項式Irreducible Polynomial」：沒有因式。質式。
greatest common divisor
多項式搭配不同數域，多項式最大公因數產生不同定義。
一、有理數多項式：輾轉相除，總是讓商數的首項係數是一。
「首一多項式Monic Polynomial」：首項係數為一。
二、整數多項式：係數事先提出最大公因數，以保證唯一解。輾轉相除，總是讓被除數乘上足夠倍率，以保證商數是整數。
「本原多項式Primitive Polynomial」：係數最大公因數為一。
a = 72x³ - 12x² - 6x¹ - 12
b = -6x³ + 22x² - 24x¹ + 8
| rational polynomial | content | primitive part
---------| ------------------------| --------| -----------------
a | 72(x-⅔)(x²+½x+¼) | 1 | 72(x-⅔)(x²+½x+¼)
b | -6(x-⅔)(x-1)(x-2) | 1 | -6(x-⅔)(x-1)(x-2)
gcd(a,b) | (x-⅔) | 1 | (x-⅔)
| integer polynomial | content | primitive part
---------| ------------------------| --------| -----------------
a | (2)(3)(3x-2)(4x²+2x+1) | (2)(3) | (3x-2)(4x²+2x+1)
b | (-1)(2)(3x-2)(x-1)(x-2) | (-1)(2) | (3x-2)(x-1)(x-2)
gcd(a,b) | (2)(3x-2) | (2) | (3x-2)
a = a₀x⁰ + a₁x¹ + a₂x² + ... polynomial a = cont(a) pp(a) integer polynomial cont(a) = gcd(a₀,a₁,a₂,...) content pp(a) = a / cont(a) primitive part cont(gcd(a,b)) = gcd(cont(a),cont(b)) commutativity pp(gcd(a,b)) = gcd(pp(a),pp(b)) commutativity gcd(a,b) = cont(gcd(a,b)) pp(gcd(a,b)) definition gcd(a,b) = gcd(cont(a),cont(b)) gcd(pp(a),pp(b)) algorithm a lc(b)^(deg(a)-deg(b)+1) = bq + r pseudo division pquo(a,b) = q pseudo quotient prem(a,b) = r pseudo remainder
遞歸法。輾轉相除。O(N²)。
加速技巧。Lehmer's GCD Algorithm。O(N²)。
分治法。位數對半分。O(M(N)logN)。
https://www.csd.uwo.ca/~mmorenom/CS424/Lectures/FastDivisionAndGcd.html/node6.html
equation
luogu P4506
equations
http://homepages.laas.fr/henrion/aime@cz11/kukelova.pdf
Gröbner Basis
Polynomial (Group)
primitive root of unity
遞增法。第n項等於(xⁿ - 1)除以其所有因數項。
「分圓多項式Cyclotomic Polynomial」：單位原根。
(xⁿ - 1) = pro cyc(d)
d|n
luogu P1520
Polynomial (Recurrence)
evaluation
已知線性遞迴函數共N項，求數列第K項。
動態規劃。求第0項到第N項。O(NK)。
遞推法。不斷展開最高項＝長除法除以特徵多項式。xᴷ模特徵多項式，求值化作求餘數。O(M(K))。
分治法。一邊倍增、一邊模。O(M(N)logK)。
「特徵多項式Characteristic Polynomial」：線性遞迴函數，移項到等號同側，再變成生成函數。
expand leading monomial ⟷ mod characteristic polynomial F₅ = F₄ + F₃ F₅ = (F₃ + F₂) + F₃ = 2F₃ + F₂ ⟷ x⁵ mod (x² - x¹ - x⁰) F₅ = 2(F₂ + F₁) + F₂ = 3F₂ + 2F₁
luogu P4723 P5808 P6115
interpolation
已知數列共K項，求線性遞迴函數共N項，項數盡量少。
遞推法。Berlekamp-Massey Algorithm＝牛頓內插。無需製作前K點都是零的多項式，而是沿用舊的特徵多項式。加速技巧：舊的特徵多項式，若長度足夠，則繼續沿用。O(NK)。
「最小多項式Minimal Polynomial」：內插，讓項數盡量少，最短的特徵多項式。
luogu P5487
Polynomial (Matrix)（Under Construction!）
多項式運算 → 矩陣運算
多項式可以視作向量、併成矩陣。有些多項式運算可以套用矩陣運算，然而矩陣運算時間複雜度較高。
唯一例外：多項式求根的演算法，利用了矩陣運算。目前沒有更快的演算法。
convolution 乘法（卷積） ⇒ Toeplitz matrix 常對角矩陣 circular convolution 循環乘法 ⇔ circulant matrix 循環矩陣 interpolation 內插 ⇔ Vandermonde Matrix 內插矩陣 Fourier transform 傅立葉轉換 ⇔ Fourier Matrix 傅立葉矩陣 greatest common divisor 最大公因數 ⇔ Sylvester Matrix 結式矩陣 linear recurrence 線性遞迴函數 ⇔ companion matrix 相伴矩陣
multiplication: Toeplitz Matrix
「常對角矩陣」。每一條左上右下斜線是相同元素的矩陣。
[ 1 5 3 2 ] [ 2 1 5 3 ] [ 7 2 1 5 ] [ 0 7 2 1 ] [ 8 0 7 2 ] [ 9 8 0 7 ]
多項式乘法（數列卷積）化作常對角矩陣求值。反之則不然。
[ 0 ]
[ 8 7 6 0 0 0 0 0 0 ] [ 0 ]
[ 0 8 7 6 0 0 0 0 0 ] [ 1 ]
(1 2 3 4 5) [ 0 0 8 7 6 0 0 0 0 ] [ 2 ]
∗ ———> [ 0 0 0 8 7 6 0 0 0 ] [ 3 ]
(6 7 8) [ 0 0 0 0 8 7 6 0 0 ] [ 4 ]
[ 0 0 0 0 0 8 7 6 0 ] [ 5 ]
[ 0 0 0 0 0 0 8 7 6 ] [ 0 ]
[ 0 ]
convolution Toeplitz matrix evaluation
[ 6 0 0 0 0 ]
[ 7 6 0 0 0 ] [ 1 ]
(1 2 3 4 5) [ 8 7 6 0 0 ] [ 2 ]
∗ ———> [ 0 8 7 6 0 ] [ 3 ]
(6 7 8) [ 0 0 8 7 6 ] [ 4 ]
[ 0 0 0 8 7 ] [ 5 ]
[ 0 0 0 0 8 ]
convolution Toeplitz matrix evaluation
加法O(N)。常對角矩陣只有2N-1種數字。
求值O(NlogN)。填充元素成為循環矩陣。
[ 1 5 3 7 2 ] [ 1 ]
[ 1 5 3 ] [ 1 ] [ 2 1 5 3 7 ] [ 2 ]
[ 2 1 5 ] [ 2 ] ———> [ 7 2 1 5 3 ] [ 3 ]
[ 7 2 1 ] [ 3 ] [ 3 7 2 1 5 ] [ 0 ]
[ 5 3 7 2 1 ] [ 0 ]
Toeplitz matrix circulant matrix
求解O(N²)。也有時間複雜度O(Nlog²N)的演算法。
https://en.wikipedia.org/wiki/Levinson_recursion
乘法O(N²)。
http://stackoverflow.com/questions/15889521/
反矩陣O(N²)。
http://www.sciencedirect.com/science/article/pii/S0893965907000535
circular multiplication: Circulant Matrix
「循環矩陣」。直向橫向皆循環的方陣。
[ 5 9 8 7 6 ] [ 6 5 9 8 7 ] [ 7 6 5 9 8 ] [ 8 7 6 5 9 ] [ 9 8 7 6 5 ]
多項式循環乘法（數列循環卷積）＝循環矩陣求值。
[ 5 9 8 7 6 ] [ 1 ]
(1 2 3 4 5) [ 6 5 9 8 7 ] [ 2 ]
⊛ <——> [ 7 6 5 9 8 ] [ 3 ]
(5 6 7 8 9) [ 8 7 6 5 9 ] [ 4 ]
[ 9 8 7 6 5 ] [ 5 ]
circular convolution circulant matrix evaluation
求值、求解、乘法、反矩陣，都是O(NlogN)。演算法是傅立葉轉換、數論轉換。
interpolation: Vandermonde Matrix
interpolation: Fourier Matrix
greatest common divisor: Sylvester Matrix
https://cp-algorithms.com/algebra/polynomial.html
a: a0 x^0 + ... + am x^m
b: b0 x^0 + ... + bn x^n
[ a0...am ] n
[ a0...am ]
S = [ a0...am ]
[ b0...bn ] m
[ b0...bn ]
ax + by = res(a,b) = det(S)
res(a,b) = det(S) = 0 iff gcd(a,b) exist
root: Companion Matrix
多項式的根＝Companion Matrix的特徵值。
時間複雜度是多項式時間。
linear recurrence evaluation: Companion Matrix
「相伴矩陣」。下方次對角線是1，最右直條是任意數，其餘元素是0的方陣。
[ 0 0 0 0 6 ] [ 1 0 0 0 7 ] [ 0 1 0 0 8 ] [ 0 0 1 0 9 ] [ 0 0 0 1 5 ]
線性函數＝矩陣。線性遞迴函數＝相伴矩陣（轉置）。
linear recursive function: f(n+3) = 7 f(n) + 8 f(n+1) + 9 f(n+2) companion matrix: [ f(n+1) ] [ 0 1 0 ] [ f(n) ] [ f(n+2) ] = [ 0 0 1 ] [ f(n+1) ] [ f(n+3) ] [ 7 8 9 ] [ f(n+2) ]
遞迴函數＝函數複合許多次。線性遞迴函數＝線性函數複合許多次＝相伴矩陣（轉置）次方。
recursive function: | linear recursive function:
{ f(0) = 1 | { f(0) = 0
{ f(n) = 2 f(n-1)² - 4 | { f(1) = 1
| { f(2) = 2
function compositions: | { f(n+3) = 7 f(n) + 8 f(n+1) + 9 f(n+2)
let g(x) = 2x² - 4 |
f(0) = 1 | companion matrix exponentiation:
f(1) = g(1) | [ f(n) ] [ 0 1 0 ]ⁿ [ f(0) ]
f(2) = g(g(1)) | [ f(n+1) ] = [ 0 0 1 ] [ f(1) ]
f(3) = g(g(g(1))) | [ f(n+2) ] [ 7 8 9 ] [ f(2) ]
線性遞迴函數的第K項＝相伴矩陣（轉置）的K次方。
矩陣次方有高速演算法，如同整數次方。
線性遞迴函數N項，求數列第K項，需要O(log(K-N))次矩陣乘法，時間複雜度O(N² + N³log(K-N))，通常簡單寫成O(N³logK)。此處假設矩陣相乘是O(N³)。
相伴矩陣k次方，得到第k項。 [ f(k) ] [ 0 1 0 ]ᵏ [ f(0) ] <----- f(k) is here [ f(k+1) ] = [ 0 0 1 ] [ f(1) ] [ f(k+2) ] [ 7 8 9 ] [ f(2) ] 相伴矩陣k-n+1次方，稍微快一點。此處k等於3。 [ f(k-2) ] [ 0 1 0 ]ᵏ⁻² [ f(0) ] [ f(k-1) ] = [ 0 0 1 ] [ f(1) ] [ f(k) ] [ 7 8 9 ] [ f(2) ] <--- f(k) is here
總結。線性遞迴函數N項，求數列第K項：
一、動態規劃，求第0項到第K項。O(NK)。
二、相伴矩陣K次方。O(N³logK)。假設矩陣相乘O(N³)。
三、多項式乘法與餘數。O(N²logK)甚至O(NlogNlogK)。
當N是常數，一變成O(K)，二三變成O(logK)。
UVa 10229 10518 10655 10743 10754 10870 12653