Polynomial

Polynomial

「多項式」。離散與連續的橋樑。

f(x) = 8x⁰ + 4x¹ + 0x² + 2x³

多項式視作數值:底數是x,位數是多項式係數。

多項式視作函數:輸入是x,輸出是多項式求值。

多項式可能無限長:有限項多項式、無限項多項式。

多項式必須搭配數域:整數多項式、有理數多項式、……。

UVa 392 126 10719 10586 10951 463 930 10428 498 10268 10105

Polynomial資料結構

Dense Polynomial:array,一格存一項。索引值是次方,內容是係數。

 ₀ ₁ ₂ ₃ ₄ ₅ ₆ ₇ ₈ ₉
 ___________________
|8|4|0|2| | | | | | |
 ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾

Sparse Polynomial:list或array,一份存一項。內容是係數和次方。係數為零的項,統統去掉,串成一串。

  ___   ___   ___
⤷|8|0|→|4|1|→|2|3|
  ‾‾‾   ‾‾‾   ‾‾‾

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
differentiation d/dx (x² + 3x + 2) = 2x + 3
integration     ∫ (x² + 3x + 2) dx = ¹⁄₃x³ + ³⁄₂x² + 2x + c
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》

專著《Modern Computer Algebra》

三本專著詳細介紹了多項式運算。以下章節不再重複整理,只做簡單介紹。

快速多項式運算

多項式運算,普通演算法是O(N²)。

幸運的是,多項式乘法,運用複數系統、餘數系統,運用快速傅立葉轉換、快速數論轉換,時間複雜度降為O(NlogN),接近線性時間。以多項式乘法為根基,各種多項式運算皆可降至接近線性時間。多項式運算從此以後變得快速。

拓展數系的三個元件:對稱(自然數到整數)、數對(整數到有理數)、循環(有理數到複數)。多項式乘法,同時運用三者,得到傅立葉矩陣,形成高速演算法,O(N²)變O(NlogN)。箇中玄妙之處,令人拍案叫絕。

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))。

加速技巧。計算a²。自己乘自己的情況下,省下一次傅立葉轉換。原本需要三次傅立葉轉換,現在只需要兩次傅立葉轉換。

加速技巧。計算ab。改為計算(a + b𝑖)²然後擷取虛部。(a + b𝑖)² = (a² - b²) + 2ab𝑖。a放實部,b放虛部,平方之後,虛部恰是2ab。原本需要三次傅立葉轉換,現在只需要兩次傅立葉轉換。

「循環多項式Circular Polynomial」:xᵃ = xᵃ⁺ⁿ。

luogu P3803

multiplicative inverse: reciprocal ¹⁄ₓ

遞推法。長除法。O(N²)。

遞推法。Hensel Lifting=位數倍增+牛頓法。O(M(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(M(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

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)讓兩者平衡。

luogu P5373

composite inverse: inverse function ⁻¹

遞推法。Series Reversion。O(N² + M(N)N)。

公式解。Lagrange Inversion。O(N² + M(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(M(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(M(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(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

high-order finite difference Δᵏ

重複K次。O(NK)。

乘以(1-x)ᵏ。次方運算。O(M(N)logK)。

乘以(1-x)ᵏ。自然對數指數,避開次方運算。O(M(N))。

乘以(1-x)ᵏ。二項式展開,遞推法求係數。O(M(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(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。我沒有研究。

luogu P5050

interpolation

基礎知識請見本站文件「Polynomial Interpolation」。

遞增法。Lagrange 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)。

arithmetic

四則運算(除法不能有餘數)、循環版本四則運算(循環除法不會有餘數),擁有容易理解的解釋方式。

多點求值、點對點四則運算、內插。O(M(N)logN)。

多點求值用複數單位根,內插用快速傅立葉轉換。O(M(N))。

Polynomial (Resolution)

equation

多項式方程式=多項式求根=Companion Matrix特徵多項式求根=Companion Matrix求特徵值。請見最後章節。

equations

Gröbner Basis。我沒有研究。

Polynomial (Divisor)

factorization

回溯法。Kronecker's Algorithm。兩個多項式整除,所有對應點也會整除。原式的函數點的某個因數,是因式的函數點。遞迴窮舉因式的函數點,以「Newton Interpolation」將函數點變成多項式;試除原式,如果整除,就是因式。時間複雜度是指數時間。

貪心法。LLL Algorithm。我沒有研究。時間複雜度是多項式時間。

「不可約多項式Irreducible Polynomial」:沒有因式。質式。

luogu P4506

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)。

primitive root of unity

遞增法。第n項是(xⁿ - 1)除以每一個n的因數項。不必除以第n項。

「分圓多項式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。原理是Newton Interpolation,但是無需製作前K點都是零的多項式,而是沿用舊的特徵多項式。O(NK)。

加速技巧。特徵多項式,若長度足夠,則繼續沿用。O(NK)。

「最小多項式Minimal Polynomial」:內插,讓項數盡量少,最短的特徵多項式。

luogu P5487

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)。