polynomial
polynomial
「多項式」。變數的各種次方、乘上倍率、通通相加。
8x⁰ - 4x¹ + 0x² + 2x³ 上升排列 2x³ + 0x² - 4x¹ + 8x⁰ 下降排列
polynomial資料結構
dense polynomial：array，一格存一項。索引值是次方，內容是係數。
₀ ₁ ₂ ₃ ₄ ₅ ₆ ₇ ₈ ₉ ___________________ |8|4|0|2| | | | | | | ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
sparse polynomial：list或array，一份存一項。內容是係數和次方。係數為零的項，統統去掉，串成一串。
___ ___ ___ ⤷|8|0|→|4|1|→|2|3| ‾‾‾ ‾‾‾ ‾‾‾
polynomial類型
多項式視作數值：底數是x，位數是多項式係數。
多項式視作函數：輸入是x，輸出是多項式求值。
多項式可能無限項：有限多項式、無限多項式。
多項式必須搭配數域：整數多項式、有理數多項式、……。
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 expansion (x + 1)² = x² + 2x + 1 factorization x² + 2x + 1 = (x + 1)² differentiation d/dx (x² + 3x + 2) = 2x + 3 integration ∫ (x² + 3x + 2) dx = ¹⁄₃x³ + ³⁄₂x² + 2x + c evaluation (x² + 3x + 2)|x=1 = 6 resolution (x² + 3x + 2) = 0 ⇒ x = -1 or -2
度數
最高次方。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³ deg(f) = 3
加法、減法
各項一一配對、係數相加。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³ g(x) = 1x⁰ + 2x¹ f(x) + g(x) = 9x⁰ - 2x¹ + 0x² + 2x³ f(x) - g(x) = 7x⁰ - 6x¹ + 0x² + 2x³
乘法
各項交叉配對、係數相乘、次方相加。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³
g(x) = 1x⁰ + 2x¹
f(x) g(x) = 8x⁰1x⁰ - 4x¹1x⁰ + 0x²1x⁰ + 2x³1x⁰
+ 8x⁰2x¹ - 4x¹2x¹ + 0x²2x¹ + 2x³2x¹
帶餘除法
帶餘除法的輸出是兩個多項式：商數、餘數。
長除法，一直補位、一直求餘數，直到餘數度數小於除數度數。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³ + 3x⁴ g(x) = 1x⁰ + 2x¹ + 3x² f(x) = (¹⁄₃x⁰ + 0x¹ + 1x²) g(x) + (²⁄₃x⁰ - ¹⁰⁄₃x¹)
無窮除法
無窮除法的輸出是一個多項式：商數。
長除法，一直補位、一直求餘數，商數越來越長。
1 - x
———————————— = 1 - 3x + 3x² + 3x³ + ...
1 + 2x + 3x²
知名範例是geometric series。a是任意數，例如1。
1/(1-ax) = (ax)⁰ + (ax)¹ + (ax)² + ... 1/(1-x) = x⁰ + x¹ + x² + ...
可能除得盡：餘數是零，商數是多項式。
可能除不盡：餘數不是零，商數是無窮多項式。
展開
多項式的連乘、多項式的次方，化作一個多項式。
f(x) = 2x⁰ (2x⁰ + 1x¹) (2x⁰ - 2x¹ + x²)
= 8x⁰ - 4x¹ + 0x² + 2x³
f(x) = (2x⁰ + 1x¹)²
= 4x⁰ + 4x¹ + 1x²
分解
展開的反方向。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³
= 2x⁰ (2x⁰ + 1x¹) (2x⁰ - 2x¹ + 1x²)
微分、積分
微分結果仍是多項式函數！dx竟然變不見！
積分結果仍是多項式函數！常數項是任意數，通常設成0。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³ d/dx f(x) = - 4x⁰ + 0x¹ + 6x² ∫ f(x) dx = 0x⁰ + 8x¹ - 2x² + 0x³ + ½x⁴
代入（求值）
Horner's rule。從最高次方項開始，一乘一加，得到答案。不需要次方運算。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³
= (((((2 ⋅ x) + 0) ⋅ x) - 4) ⋅ x) + 8
求根（求解）
用途不明。
f(x) = 8x⁰ - 4x¹ + 0x² + 2x³ = 0 x = -2 or 1+𝑖 or 1-𝑖
快速多項式運算
多項式運算，普通演算法是O(N²)。
幸運的是，多項式乘法，運用複數系統、餘數系統，運用快速傅立葉轉換、快速數論轉換，時間複雜度降為O(NlogN)，相當接近一次方時間O(N)。以多項式乘法為根基，各種多項式運算皆可降為接近一次方時間。多項式運算從此以後變得快速。
拓展數系的三個元件：對稱（自然數到整數）、數對（整數到有理數）、循環（有理數到複數）。多項式乘法，同時運用三者，得到傅立葉矩陣，形成高速演算法，O(N²)變O(NlogN)。箇中玄妙之處，令人拍案叫絕。
專著《Algorithmes Efficaces en Calcul Formel》。
專著《Algorithms for Computer Algebra》。
三本專著詳細介紹了快速多項式運算。以下章節不再重複整理，只做簡單介紹。
UVa 392 126 10719 10586 10951 463 930 10428 498 10268 10105
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)。
加速技巧。計算a²。自己乘自己的情況下，省下一次傅立葉轉換。原本需要三次傅立葉轉換，現在只需要兩次傅立葉轉換。
加速技巧。計算ab。改為計算(a + b𝑖)²然後擷取虛部。(a + b𝑖)² = (a² - b²) + 2ab𝑖。a放實部，b放虛部，平方之後，虛部恰是2ab。原本需要三次傅立葉轉換，現在只需要兩次傅立葉轉換。
乘法是各種運算的根基，時間複雜度簡寫成O(𝓜(N))。
「循環多項式circular polynomial」：xᵃ = xᵃ⁺ⁿ。
luogu P3803
multiplicative inverse: reciprocal ¹⁄ₓ
遞推法。長除法。O(N²)。
遞推法。Hensel lifting＝位數倍增＋牛頓法。O(𝓜(N))。
「截斷多項式truncated polynomial」：低次項N項。模xᴺ。
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(𝓜(N))。
「餘式polynomial residue」：餘數運算推廣成餘式運算。截斷多項式可以轉換到餘式系統進行處理：令模式是xⁿ。
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
infinite division /
已知多項式分式，求得商的第K項。
遞推法。長除法。O(NK)。
分治法。Bostan–Mori algorithm。利用分母的對稱函數，將分式擴分再拆分，得到兩個分式，其商恰是偶數項與奇數項。第K項位於其中一邊。O(𝓜(N)logK)。
luogu P4723
polynomial - function
evaluation ()
遞增法。逐項次方、逐項相加。O(NlogN)。
遞推法。Horner's rule。O(N)。
composition ∘
遞推法。Horner's rule。O(𝓜(AB)A)。
分治法。A對半分。O(𝓜(AB)logA)。
加速技巧。過程都在頻域。O(𝓜(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(𝓜(C)A)。
分治法。A對半分。O(𝓜(C)(AB/C)(1+log(min(A,ceil(C/B))))。簡單寫成O(𝓜(C)A)或O(𝓜(C)BlogA)，後者需C≥A。
加速技巧。Brent–Kung algorithm。B分成低次項k項、高次項B-k項。泰勒展開，只取低階C/k項。先求最低階，再以遞推法求得各階導數、冪數、分母。O(𝓜(C)(klogA + C/k))。令klogA = C/k讓時間複雜度達到最低。O(𝓜(C)sqrt(ClogA))。
把B變小（長度k），分治法開頭變快（長度到達C之前），代價是泰勒展開（長度C/k）。令k = sqrt(C/logA)讓兩者平衡。
luogu P5373
composite inverse: inverse function ⁻¹
遞推法。series reversion。O(N² + 𝓜(N)N)。
公式解。Lagrange inversion。O(N² + 𝓜(N)sqrtN)。
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()
公式解。多項式乘法視作數列卷積，拆出最後一項。O(N²)。
公式解。當輸入的常數項為1，則輸出的常數項為0——得以使用積分。當輸入的常數項不為1，則提出倍率使之為1——倍率是常數項。O(𝓜(N))。
when a > 0 ln(a) = (a-1)¹/1 - (a-1)²/2 + (a-1)³/3 - ...
when a > 0
ln(a) = b
ln′(a) a′ = b′
(1/a) a′ = b′
a′ = ab′
n aₙ xⁿ⁻¹ = (aₙ xⁿ)(n bₙ xⁿ⁻¹)
n aₙ = sum {i aₙ₋ᵢ bₙ}
1≤i≤n
aₙ = sum {i aₙ₋ᵢ bₙ} / n
1≤i≤n
aₙ - sum {i aₙ₋ᵢ bₙ} / n = a₀ bₙ
1≤i≤n-1
(aₙ - sum {i aₙ₋ᵢ bₙ} / n) / a₀ = bₙ
1≤i≤n-1
when a > 0 and a₀ = 1 ln(a) = b ln′(a) a′ = b′ (1/a) a′ = b′ ∫((1/a) a′) = b (if a₀ = 1 then b₀ = 0) when a > 0 and a₀ ≠ 1 ln(a) = ln((a/a₀) a₀) = ln(a/a₀) + ln(a₀)
luogu P4725
natural exponentiation exp()
公式解。自然對數公式解，最後不做移項。O(N²)。
遞推法。Hensel lifting。O(𝓜(N))。
bₙ = sum {i bₙ₋ᵢ aₙ} / n
1≤i≤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(𝓜(N)logK)。
公式解。O(𝓜(N))。
b = aᵏ ln(b) = k ln(a) b = exp(k ln(a))
luogu P5245 P5273
square root sqrt()
遞推法。Hensel lifting。O(𝓜(N))。
公式解。O(𝓜(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(𝓜(N))。
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𝑖
arcsin(x) = ∫ (x′ / sqrt(1 - x²)) dx arccos(x) = -∫ (x′ / sqrt(1 - x²)) dx arctan(x) = ∫ (x′ / (1 + x²)) dx
luogu P5264 P5265
high-order finite difference Δᵏ
重複K次。O(NK)。
乘以(1-x)ᵏ。次方運算。O(𝓜(N)logK)。
乘以(1-x)ᵏ。自然對數指數，避開次方運算。O(𝓜(N))。
乘以(1-x)ᵏ。二項式展開，遞推法求係數。O(𝓜(N))。
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 - evaluation
multipoint evaluation
遞推法。Horner's rule。O(NK)。
遞增法。求值a(xᵢ)化作求餘數a mod (x-xᵢ)。O(𝓜(N)K)。
分治法。項數對半分、形成二元樹。O(𝓜(N) + 𝓜(K)logK)。
一、一次多項式連乘M(x₁,x₂,...) = (x-x₁)(x-x₂)...(x-xK)。項數對半分M(all) = M(low)M(high)。形成二元樹，從葉往根乘。
二、求餘數，一回合擴展成logK回合。第一回合，模數是全部連乘M(all)。第二回合，模數是前半部M(low)、後半部M(high)。最後回合，模數是(x-xᵢ)。形成二元樹，從根往葉模。
加速技巧。Tellegen's principle。O(𝓜(N) + 𝓜(K)logK)。
luogu P5050
interpolation
基礎知識請見本站文件「polynomial interpolation」。
遞增法。Lagrange interpolation。O(𝓜(N)N)。
一、一次多項式連乘M(x₁,x₂,...) = (x-x₁)(x-x₂)...(x-xK)。
二、內插P(x₁,x₂,...)化作求商數M(x₁,x₂,...) / (x-xᵢ)。
分治法。項數對半分、形成二元樹。O(𝓜(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(𝓜(N) + 𝓜(K)logK)。
arithmetic
四則運算（除法不能有餘數）、循環版本四則運算（循環除法不會有餘數），擁有容易理解的解釋方式。
多點求值、點對點四則運算、內插。O(𝓜(N)logN)。
多點求值用複數單位根，內插用快速傅立葉轉換。O(𝓜(N))。
polynomial - resolution
root
遞推法。Bairstow's method。複數共軛根形成實數二次因式。原式除以二次因式得到一次餘式。令一次餘式是零。牛頓法求得二次因式的兩個係數。時間複雜度O(N²T)，收斂速度是二次方。
遞推法。Laguerre's method。原式是一次因式連乘。一次導數是一次因式導數連加。二次導數是一次因式平方連減。以此推導遞推公式。時間複雜度O(N²T)，收斂速度是三次方。
遞推法。companion matrix。多項式求根＝companion matrix特徵多項式求根＝companion matrix求特徵值。請見最後章節。時間複雜度O(N²T)，不必猜測初始值。MATLAB採用此方法。
equation
多項式方程式＝多項式求根。
遞迴加權總和【數學家稱作ideal】
遞迴加權總和：給定集合，窮舉各種權重組合，得到各種加權總和。
A = {a₁, a₂, a₃}
a₄ = w₁a₁ + w₂a₂ + w₃a₃
a₅ = v₁a₁ + v₂a₂ + v₃a₃ + v₄a₄
a₆ = u₁a₁ + u₂a₂ + u₃a₃ + u₄a₄ + u₅a₅
: :
數學式子展開之後，其實等同於沒有遞迴的加權總和。
span(A) = k₁a₁ + k₂a₂ + k₃a₃
多項式的遞迴加權總和：沒有差別。
A = {x³-2xy, 2x²y-y²+x, x²+1}
a₄ = 2y (x³-2xy) + x (2x²y-y²+x) + 0 (x²+1) = -3xy² - x²
a₅ = y (x³-2xy) + 0 (2x²y-y²+x) + xy (x²+1) + (x-y+1) (-3xy²-x²)
: :
遞迴餘數【尚無正式名稱】
遞迴餘數：給定模數集合，窮舉各種模數組合，得到各種餘數。
M = {m₁, m₂, m₃}
a ⤳ b (mod m₁)
b ⤳ c (mod m₃) —→ a ⤳ r (rmod M)
c ⤳ d (mod m₁)
: :
q ⤳ r (mod m₂)
遞迴餘數的性質，宛如一般餘數的性質。
a₁ ⤳ r₁ (rmod M) a₂ ⤳ r₂ (rmod M) a₁ + a₂ ⤳ r₁ + r₂ (rmod M)
多項式的遞迴餘數：不斷消滅首項，度數越來越小，最終無法再模。
項有多種排序方式。固定使用同一種，哪一種都可以。
字典順序：先以字母排序，再以度數排序。（部分變數視作常數）
x²y² + x²y + x² + xy² + xy + x + y² + y + 1
^^^^^^^^^^^^^^^ ^^^^^^^^^^^^ ^^^^^^^^^^
x² x 1
度數順序：先以度數排序，再以字母排序。
x²y² + x²y + xy² + x² + xy + y² + x + y + 1
^^^^ ^^^^^^^^^ ^^^^^^^^^^^^ ^^^^^ ^
4 3 2 1 0
Gröbner basis
給定被除數集合，實施遞迴加權總和；找到除數集合，可以遞迴整除每一個加權總和。這樣的除數集合，稱作Gröbner basis。
M is Gröbner basis of A: aᵢ ⤳ 0 (rmod M) for all aᵢ in span(A)
遞迴加權總和，總共無限多個，無法一一檢查是否整除。必須取巧。
遞推法。Buchberger's algorithm。被除數集合納入除數集合。窮舉兩兩除數，計算S-polynomial，一一檢查是否整除。如果都被整除：當前除數集合就是Gröbner basis。如果無法整除：餘數納入除數集合，重頭檢查。時間複雜度是指數時間。
遞推法。M4GB。目前最快的演算法。
lcm(lm(a),lm(b)) lcm(lm(a),lm(b))
s(a,b) = ———————————————— a - ———————————————— b
lm(a) lm(b)
lm(): leading monomial
lcm(): least common multiple
a = x³ - 2xy b = 2x²y - y² + x lm(a) = x³ lm(b) = 2x²y lcm(lm(a),lm(b)) = 2x³y s(a,b) = (2x³y / x³) a - (2x³y / 2x²y) b = -3xy² - x²
equations
Gröbner basis。我沒有研究。求解工具
polynomial - divisor
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²)。
分治法。half-GCD algorithm。位數對半分。O(𝓜(N)logN)。
factorization
回溯法。Kronecker's algorithm。兩個多項式整除，所有對應點也會整除。原式的函數點的某個因數，是因式的函數點。遞迴窮舉因式的函數點，以「Newton interpolation」將函數點變成多項式；試除原式，如果整除，就是因式。時間複雜度是指數時間。
貪心法。LLL algorithm。我沒有研究。時間複雜度是多項式時間。
「不可約多項式irreducible polynomial」：除了自己和1以外，沒有其他因式。質式。
luogu P4506
square-free factorization
擷取指數是一的因式連乘積。換句話說，捨棄指數是平方、立方、……的因式。
注意到，只求得因式連乘積，不求得各個因式。
a(x) = (1+x)¹(1-2x+2x²)¹(1-x)²(3-2x)²(2-x+x²)²x³(2x-1)³
^^^^^^^^^^^^^^^^^
square-free factor
甚至更進一步，一邊擷取指數是一的因式連乘積，一邊實施微分降低指數，求得各種指數的因式連乘積。換句話說，所有因式，依照指數，進行分類。
注意到，只求得因式連乘積暨指數，不求得各個因式。
a(x) = (1+x)¹(1-2x+2x²)¹ (1-x)²(3-2x)²(2-x+x²)² x³(2x-1)³
^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^
遞推法。gcd(a,a′)可以降低指數，而且不影響係數。時間複雜度是多項式時間。
a = a₁¹ a₂² a₃³ ... 各種指數的因式連乘積
a′ = (1 a₁¹⁻¹ a₁′) ( a₂² a₃³ ...) 微分
+ (2 a₂²⁻¹ a₂′) (a₁¹ a₃³ ...)
+ (3 a₃³⁻¹ a₃′) (a₁¹ a₂² ...)
+ ......
b = gcd(a,a′) = a₁¹⁻¹ a₂²⁻¹ a₃³⁻¹ ... 次方減一
c = a/b = a₁ a₂ a₃ ... 次方變一
d = gcd(b,c) = a₂ a₃ ... 得到指數大於一的因式連乘積
f = c/d = a₁ 得到指數為一的因式連乘積！
b/d = gcd(b,b′) = a₂²⁻² a₃³⁻² ... 開始第二回合，次方再減一
: : :
遞推法。Yun's algorithm。過程有一點類似正規正交化。時間複雜度是多項式時間。
a = a₁¹ a₂² a₃³ ... 各種指數的因式連乘積
a′ = (1 a₁¹⁻¹ a₁′) ( a₂² a₃³ ...) 微分
+ (2 a₂²⁻¹ a₂′) (a₁¹ a₃³ ...)
+ (3 a₃³⁻¹ a₃′) (a₁¹ a₂² ...)
+ ......
b = gcd(a,a′) = a₁¹⁻¹ a₂²⁻¹ a₃³⁻¹ ... 次方減一
c = a/b = a₁ a₂ a₃ ... 次方消失
p = a′/b = (1 a₁′) ( a₂ a₃ ...)
+ (2 a₂′) (a₁ a₃ ...)
+ (3 a₃′) (a₁ a₂ ...)
+ ......
q = (a/b)′ = (a₁′) ( a₂ a₃ ...)
+ (a₂′) (a₁ a₃ ...)
+ (a₃′) (a₁ a₂ ...)
+ ......
r = p - q = (1 a₂′) (a₁ a₃ ...)
+ (2 a₃′) (a₁ a₂ ...)
+ ......
gcd(c,r) = a₁
algorithm = Gram–Schmidt orthonormalization square-free = linear dot product = gcd derivative = shearing basis = factor with incremental exponent
「無平方多項式square-free polynomial」：不存在因式，其指數是二次（以上）。注意到，是指數（重複次數），不是度數（最高次方值）。
primitive root of unity
「分圓多項式cyclotomic polynomial」：單位原根連乘積。
遞推法。利用N的因數。O(𝓜(N)σ₀(N)) ≈ O(𝓜(N)sqrtN)。
(xⁿ - 1) = prod Φd(x)
d|n
Φₙ(x) = (xⁿ - 1) / prod Φd(x)
d|n
d≠n
遞推法。利用N的質因數。O(𝓜(N)ω(N)) ≈ O(𝓜(N)logN)。
n = p₁e₁ × p₂e₂ × ... × pₖeₖ f₀(x) = (x-1) fᵢ(x) = fᵢ₋₁(xpᵢ) / fᵢ₋₁(x) Φₙ(x) = fₖ(xn/(pᵢ×...×pₖ))
luogu P1520
polynomial - recurrence
evaluation
已知線性遞迴函數共N項，求數列第K項。
「特徵多項式characteristic polynomial」：線性遞迴函數，移項到等號同側，再變成生成函數。
動態規劃。依序計算每一項。O(NK)。
公式解。Fiduccia's algorithm。不斷展開最高項＝長除法除以特徵多項式。xᴷ模特徵多項式，求餘數。O(𝓜(K))。
加速技巧。一邊倍增、一邊模。O(𝓜(N)logK)。
linear recursive function (degree N) c₀f(n) + c₁f(n+1) + c₂f(n+2) + ... + cɴ₋₁f(n+(N-1)) = f(n+N) generating function c₀xⁿ + c₁xⁿ⁺¹ + c₂xⁿ⁺² + ... + cɴ₋₁xn+(N-1) = xⁿ⁺ᴺ characteristic polynomial -c₀x⁰ + -c₁x¹ + -c₂x² + ... + -cɴ₋₁xN-1 + xᴺ expansion of leading monomial ←—→ long division by characteristic polynomial F₅ = F₄ + F₃ F₅ = (F₃ + F₂) + F₃ = 2F₃ + F₂ ←—→ x⁵ mod (x² - x¹ - x⁰) F₅ = 2(F₂ + F₁) + F₂ = 3F₂ + 2F₁
分治法。Bostan–Mori algorithm。線性遞迴數列第K項化作分式第K項。O(𝓜(N)logK)。
luogu P4723 P5808 P6115
interpolation
已知數列共K項，求線性遞迴函數共N項，N盡量少。
「最小多項式minimal polynomial」：度數最少的特徵多項式。
遞推法。Berlekamp–Massey algorithm。原理是Newton interpolation，但是無需製作前N點都是零的多項式，而是沿用舊的特徵多項式。O(NK)。
加速技巧。特徵多項式，若長度足夠，則繼續沿用。O(NK)。
公式解。Brent–Gustavson–Yun algorithm。生成函數a(x)、特徵多項式頭尾顛倒rev(c(x))，兩者相乘，中間均零，頭尾非零。模xᴷ以留頭去尾。餘數即是強行計算頭N項初始值。藉由餘數的度數以推敲N是多少。餘數是gcd(a(x),xᴷ)的倍數，正確答案是a(x)所搭配的倍率。擴展輾轉相除法，逐步升高倍率的度數，直到超過餘數的度數。O(𝓜(K)logN)。
sequence (length K) a₀ a₁ a₂ ... aK-1 generating function a(x) = a₀x⁰ + a₁x¹ + a₂x² + ... + aK-1xK-1 linear recursive function (degree N) c₀f(n) + c₁f(n+1) + ... + cɴ₋₁f(n+(N-1)) = f(n+N) characteristic polynomial c(x) = -c₀x⁰ + -c₁x¹ + ... + -cɴ₋₁xN-1 + xᴺ reverse characteristic polynomial rev(c(x)) = -c₀xᴺ + -c₁xN-1 + ... + -cɴ₋₁x + x⁰ interpolation a(x) rev(c(x)) ≡ r(x) (mod xᴷ) where deg(r(x)) < deg(c(x)) Bézout's identity a(x) rev(c(x)) + xᴷ q(x) = r(x) where r(x) = gcd(a(x), xᴷ) s(x) extended Euclidean algorithm given a(x) and xᴷ find rev(c(x)) and q(x) and r(x) such that deg(c(x)) is minimized and deg(c(x)) > deg(r(x))
luogu P5487 UVa 1539
polynomial - integer
整數運算→多項式運算
整數可以視作多項式。有些整數運算可以套用多項式運算。
比方來說，整數乘法套用多項式乘法。因為3×37 = 111，所以3(3x + 7) = 9x + 21，底數x = 10。
記得認真考慮整數運算的時間複雜度。快速傅立葉轉換內含整數乘法，遞迴使用整數乘法，時間複雜度成為O(N logN loglogN)。
順帶一提，整數乘法的時間複雜度下限仍是懸案。2019年出現新演算法，時間複雜度降為O(NlogN)，很可能就是下限了。
luogu P1601 P1919 P5432
整數運算↛多項式運算
有些整數運算則不可以套用多項式運算。
比方來說，整數分解套用多項式分解，雖然111 = 3×37，但是x² + x + 1 ≠ 3(3x + 7)。我們不知道底數x應該是多少，也不知道是否有減法。
順帶一提，整數分解的時間複雜度類別仍是懸案，目前還不確定它究竟是P問題或是NP-complete問題，相當特別。
polynomial - matrix
多項式運算 → 矩陣運算
有些多項式運算可以套用矩陣運算。
convolution 乘法（卷積） ⇒ Toeplitz matrix 常對角矩陣求值 circular convolution 循環乘法 ⇔ circulant matrix 循環矩陣求值 multipoint evaluation 多點求值 ⇔ Vandermonde matrix 范德蒙矩陣求值 multipoint evaluation 多點求值特例 ⇔ Fourier matrix 傅立葉矩陣求值 coprime 互質 ⇔ Sylvester matrix 結式矩陣行列式 root 根 ⇔ companion 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
加法。常對角矩陣只有2N-1種數字。時間複雜度O(N)。
求值。填充元素成為循環矩陣。時間複雜度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)。
乘法。時間複雜度O(N²)。
反矩陣。時間複雜度O(N²)。
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)。演算法是傅立葉轉換、數論轉換。
multipoint evaluation ⇨ Vandermonde matrix
「范德蒙矩陣」。橫條是任意數的0次方到N-1次方。
⎡ 9⁰ 9¹ 9² 9³ 9⁴ ⎤ ⎢ 0⁰ 0¹ 0² 0³ 0⁴ ⎥ ⎢ 3⁰ 3¹ 3² 3³ 3⁴ ⎥ ⎢ 5⁰ 5¹ 5² 5³ 5⁴ ⎥ ⎣ 2⁰ 2¹ 2² 2³ 2⁴ ⎦
多項式多點求值＝范德蒙矩陣求值。時間複雜度O(NK)。
多項式內插＝范德蒙矩陣求解。時間複雜度O(N³)。
f(x) = c₀x⁰ + c₁x¹ + c₂x² + c₃x³ + c₄x⁴ ⎡ x₀⁰ x₀¹ x₀² x₀³ x₀⁴ ⎤ ⎡ c₀ ⎤ ⎡ f(x₀) ⎤ ⎢ x₁⁰ x₁¹ x₁² x₁³ x₁⁴ ⎥ ⎢ c₁ ⎥ ⎢ f(x₁) ⎥ ⎢ x₂⁰ x₂¹ x₂² x₂³ x₂⁴ ⎥ ⎢ c₂ ⎥ = ⎢ f(x₂) ⎥ ⎢ x₃⁰ x₃¹ x₃² x₃³ x₃⁴ ⎥ ⎢ c₃ ⎥ ⎢ f(x₃) ⎥ ⎣ x₄⁰ x₄¹ x₄² x₄³ x₄⁴ ⎦ ⎣ c₄ ⎦ ⎣ f(x₄) ⎦
multipoint evaluation ⇨ Fourier matrix
「傅立葉矩陣」可以視作范德蒙矩陣的特例：輸入們恰是單位根（複數平面單位圓的等分點），而且是方陣K = N。
f(x) = c₀x⁰ + c₁x¹ + c₂x² + c₃x³ + c₄x⁴ ⎡ x₀⁰ x₀¹ x₀² x₀³ x₀⁴ ⎤ ⎡ c₀ ⎤ ⎡ f(x₀) ⎤ ⎢ x₁⁰ x₁¹ x₁² x₁³ x₁⁴ ⎥ ⎢ c₁ ⎥ ⎢ f(x₁) ⎥ ⎢ x₂⁰ x₂¹ x₂² x₂³ x₂⁴ ⎥ ⎢ c₂ ⎥ = ⎢ f(x₂) ⎥ ⎢ x₃⁰ x₃¹ x₃² x₃³ x₃⁴ ⎥ ⎢ c₃ ⎥ ⎢ f(x₃) ⎥ ⎣ x₄⁰ x₄¹ x₄² x₄³ x₄⁴ ⎦ ⎣ c₄ ⎦ ⎣ f(x₄) ⎦ where xᵢ = e-𝑖(2π/N)⋅i
coprime ⇨ Sylvester matrix
「西爾維斯特矩陣」或「結式矩陣」。多項式係數逐步向右移位，其餘元素是0，方陣。
多項式互質＝結式矩陣行列式非零。時間複雜度O(N³)。
a(x) = a₀x⁰ + a₁x¹ + ... + aₘxᵐ
b(x) = b₀x⁰ + b₁x¹ + ... + bₙxⁿ
⎡ a₀ ... aₘ ⎤ \
⎢ a₀ ... aₘ ⎥ n rows
⎢ : ⎥ |
S = ⎢ a₀ ... aₘ ⎥ /
⎢ b₀ ... bₙ ⎥ \
⎢ b₀ ... bₙ ⎥ m rows
⎢ : ⎥ |
⎣ b₀ ... bₙ ⎦ /
0. define resultant res(a,b) = det(S)
1. gcd(a,b) is not constant polynomial iff res(a,b) = det(S) = 0
2. deg(gcd(a,b)) = m + n - rank(S)
3. exist x and y such that ax + by = res(a,b) = det(S)
where deg(x) < deg(b) and deg(y) < deg(a)
root ⇨ 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 ⎦
多項式的根＝同伴矩陣特徵值。時間複雜度O(N²(N+T))，T是遞推次數。
f(x) = c₀x⁰ + c₁x¹ + c₂x² + c₃x³ + c₄x⁴
⎡ 0 0 0 0 -c₀ ⎤
⎢ 1 0 0 0 -c₁ ⎥
C = ⎢ 0 1 0 0 -c₂ ⎥
⎢ 0 0 1 0 -c₃ ⎥
⎣ 0 0 0 1 -c₄ ⎦
root(f) = root(det(C - λI)) = eigenvalue(C)
linear recurrence evaluation ⇨ companion matrix
線性函數＝矩陣。線性遞迴函數＝同伴矩陣（轉置）。線性遞迴函數的第K項＝同伴矩陣（轉置）的K次方。
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) ⎦
線性遞迴函數共N項，求數列第K項：
一、動態規劃，求第0項到第K項。O(NK)。
二、同伴矩陣K次方。O(N³logK)。假設矩陣相乘O(N³)。
三、xᴷ模特徵多項式。O(N²logK)甚至O(NlogNlogK)。
當N是常數，一變成O(K)，二三變成O(logK)。