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》

專著《Modern 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)。