recurrence

「遞迴數列」。以遞迴函數得到的數列。

recursive function:       |  f(0) = 1
⎰ f(0) = 1                |  f(1) = 2 f(0)² - 4 = -2
⎱ f(n) = 2 f(n-1)² - 4    |  f(2) = 2 f(1)² - 4 = 4
                          |  f(3) = 2 f(2)² - 4 = 28
recursive sequence:       |  f(4) = 2 f(3)² - 4 = 1564
1 -2 4 28 1564 ......     |    :        :          :

數學當中，遞迴數列與遞迴函數一體兩面，同稱recurrence。計算學當中，則是各飾一角。本文為了方便講解，暫將recurrence視作遞迴數列。

linear recurrence（linear feedback shift register）

「線性遞迴數列」。以線性遞迴函數得到的數列。

⎧ f(0) = c₀
⎨ f(1) = c₁
⎪   :     :
⎩ f(n) = a₁ f(n-1) + a₂ f(n-2) + ... + aₖ f(n-k)

知名範例是Fibonacci sequence。往前兩項、係數皆是一。

recursive function:       |  f(0) = 0
⎧ f(0) = 0                |  f(1) = 1
⎨ f(1) = 1                |  f(2) = f(1) + f(0) = 1
⎩ f(n) = f(n-1) + f(n-2)  |  f(3) = f(2) + f(1) = 2
                          |  f(4) = f(3) + f(2) = 3
recursive sequence:       |  f(5) = f(4) + f(3) = 5
0 1 1 2 3 5 ...           |    :      :      :    :

演算法（dynamic programming）

從頭開始一項一項算。線性遞迴函數K項，求數列第N項，時間複雜度O(K + (N-K)K)，通常簡單寫成O(NK)。

const int N = 10; int f[N+1]; void Fibonacci_sequence() { f[0] = 0; f[1] = 1; for (int i=2; i<=N; i++) f[i] = f[i-1] + f[i-2]; }

演算法（companion matrix）

線性函數＝矩陣。線性遞迴函數＝同伴矩陣（轉置）。

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 a(x) = 2x² - 4      |
f(0) = 1                |  companion matrix exponentiation:
f(1) = a(1)             |  ⎡ f(n)   ⎤   ⎡ 0 1 0 ⎤ⁿ ⎡ f(0) ⎤
f(2) = a(a(1))          |  ⎢ f(n+1) ⎥ = ⎢ 0 0 1 ⎥  ⎢ f(1) ⎥
f(3) = a(a(a(1)))       |  ⎣ f(n+2) ⎦   ⎣ 7 8 9 ⎦  ⎣ f(2) ⎦

線性遞迴函數的第N項＝同伴矩陣（轉置）的N次方。

同伴矩陣n次方，得到第n項。
⎡ f(n)   ⎤   ⎡ 0 1 0 ⎤ⁿ ⎡ f(0) ⎤   <----- f(n) is here
⎢ f(n+1) ⎥ = ⎢ 0 0 1 ⎥  ⎢ f(1) ⎥
⎣ f(n+2) ⎦   ⎣ 7 8 9 ⎦  ⎣ f(2) ⎦

同伴矩陣n-k+1次方，稍微快一點。此處k等於3。
⎡ f(n-2) ⎤   ⎡ 0 1 0 ⎤ⁿ⁻² ⎡ f(0) ⎤
⎢ f(n-1) ⎥ = ⎢ 0 0 1 ⎥    ⎢ f(1) ⎥
⎣ f(n)   ⎦   ⎣ 7 8 9 ⎦    ⎣ f(2) ⎦   <--- f(n) is here

矩陣次方有高速演算法，如同整數次方。

線性遞迴函數K項，求數列第N項，需要O(log(N-K))次矩陣乘法，時間複雜度O(K² + K³log(N-K))，通常簡單寫成O(K³logN)。此處假設矩陣相乘是O(K³)。

int Fibonacci_sequence(int n) { if (n < 0) return -n & 1; // 初始向量。O(K)。 // f(0) = 0 , f(1) = 1 Matrix v0(2, 1); v0[0][0] = 0; v0[0][1] = 1; // 初始向量已經有第N項 if (n < 2) return v0[0][n]; // 同伴矩陣。O(K²)。 // f(n+2) = 1 ⋅ f(n) + 1 ⋅ f(n+1) Matrix A(2, 2); A[0][0] = 0; A[0][1] = 1; A[1][0] = 1; A[1][1] = 1; // 同伴矩陣n次方。O(K³logN)。 // Matrix v(2, 1); // v = (A ^ n) * v0; // retrun v[0][0]; // 同伴矩陣n-k+1次方，稍微快一點。O(K³log(N-K))。 Matrix v(2, 1); v = (A ^ (n-1)) * v0; return v[0][1]; }

UVa 10229 10518 10655 10743 10754 10870 12653

演算法（characteristic polynomial）

特殊的生成函數。f(n)變成xⁿ。

linear recursive function:
f(n+3) = 9 f(n+2) + 8 f(n+1) + 7 f(n)

generating function:
xⁿ⁺³ = 9xⁿ⁺² + 8xⁿ⁺¹ + 7xⁿ

特殊的特徵多項式。移項、刪除變數n。

characteristic polynomial:
x³ - 9x² - 8x¹ - 7x⁰

展開最高項＝長除法求得餘數。

expand leading term:
  f(5)
= 9 f(4) + 8 f(3) + 7 f(2)
= 89 f(3) + 79 f(2) + 63 f(1)
= 880 f(2) + 775 f(1) + 623 f(0)

divide characteristic polynomial:
  (x⁵) 
⇒ (x⁵)                 ÷ (x³-9x²-8x-7) = x² ... 9x⁴ + 8x³ + 7x²
⇒ (9x⁴ + 8x³ + 7x²)    ÷ (x³-9x²-8x-7) = 9x ... 89x³ + 79x² + 63x¹
⇒ (89x³ + 79x² + 63x¹) ÷ (x³-9x²-8x-7) = 89 ... 880x² + 775x¹ + 623x⁰

展開第n項＝xⁿ模特徵多項式。最後再代入實際數值。

線性遞迴函數K項，求數列第N項，使用多項式除法，時間複雜度O(K(N-K) + K)，通常簡單寫成O(NK)。

多項式除法還可以更快。請見本站文件「polynomial」。

nonhomogeneous linear recurrence

線性遞迴函數有常數項。其實可以把常數消掉。

第n式與第n-1式，相減移項。最後多了一項。

⎰ f(n)   = a₁ f(n-1) + a₂ f(n-2) + ... + aₖ f(n-k)   + b
⎱ f(n-1) = a₁ f(n-2) + a₂ f(n-3) + ... + aₖ f(n-k-1) + b

⎰ f(n)   - a₁ f(n-1) - a₂ f(n-2) - ... - aₖ f(n-k)   = b
⎱ f(n-1) - a₁ f(n-2) - a₂ f(n-3) - ... - aₖ f(n-k-1) = b

  f(n)   - a₁ f(n-1) - a₂ f(n-2) - ... - aₖ f(n-k)
= f(n-1) - a₁ f(n-2) - a₂ f(n-3) - ... - aₖ f(n-k-1)

  f(n)   = a₁ f(n-1) + a₂ f(n-2) + ... + aₖ f(n-k)
         +    f(n-1) - a₁ f(n-2) - ... - aₖ₋₁ f(n-k) - aₖ f(n-k-1)

f(n) = (a₁+1) f(n-1) + (a₂-a₁) f(n-2) + (a₃-a₂) f(n-3) + ...
       ... + (aₖ-aₖ₋₁) f(n-k) - aₖ f(n-k-1)

linear recurrence system

多個線性遞迴函數，共享變數。

⎰ f(n+1) = a₀ f(n) + a₁ f(n-1) + ... + b₀ g(n) + b₁ g(n-1) + ...
⎱ g(n+1) = c₀ f(n) + c₁ f(n-1) + ... + d₀ g(n) + d₁ g(n-1) + ...

似乎可以把方程組化作方程式。【尚待確認】

steady state

「穩態」。數列最終趨近的一個固定數字。可能不存在。

recurrence:
f(n) = a₁ f(n-1) + a₂ f(n-2) + ... + aₖ f(n-k) + b

steady state:
f̄ = f(n)   when n→∞

求得穩態：遞迴函數的變數，一律換成穩態，解方程式。

判斷穩態是否存在：分母不能等於零。諸如此類的。

steady state:
f̄ = f(n) = f(n-1) = f(n-2) = ... = f(n-k)   when n→∞

equation:
f̄ = a₁ f̄ + a₂ f̄ + ... + aₖ f̄ + b
(1 - a₁ - a₂ - ... - aₖ) f̄ = b
f̄ = b / (1 - a₁ - a₂ - ... - aₖ)

steady state exists iff 1 - a₁ - a₂ - ... - aₖ ≠ 0

Fibonacci sequence

0 1 1 2 3 5 8 ...

⎧ f(0) = 0
⎨ f(1) = 1
⎩ f(n) = f(n-1) + f(n-2)

UVa 495 11582 luogu P3986

generating function / characteristic equation

(0 1 1 2 3 5 8 ...)  ←—→  x/(1-x-x²)

G(x) = 0x⁰ + 1x¹ + 1x² + 2x³ + 3x⁴ + 5x⁵ + 8x⁶ + ...
G(x) = f(0)x⁰ + f(1)x¹ + f(2)x² + f(3)x³ + ...
G(x) = f(0)x⁰ + f(1)x¹ + [f(0) + f(1)]x² + [f(1) + f(2)]x³ + ...
G(x) = f(0)x⁰ + f(1)x¹ + f(0)x² + f(1)x² + f(1)x³ + f(2)x³ + ...
G(x) = f(0)x⁰ + f(1)x¹ + [f(0)x² + f(1)x³ + ...] + [f(1)x² + f(2)x³ + ...]
G(x) = f(0)x⁰ + f(1)x¹ + [G(x)x²] + [G(x)x¹ - f(0)x¹]
G(x) =      0 +     x¹ + [G(x)x²] + [G(x)x¹ -      0]
G(x) = x + G(x)x + G(x)x²
G(x) - G(x)x - G(x)x² = x
G(x)(1-x-x²) = x
G(x) = x/(1-x-x²)

f(n) = f(n-1) + f(n-2)       when n ≥ 2
f(n) - f(n-1) - f(n-2) = 0   when n ≥ 2
(G(x) - f(0)x⁰ - f(1)x¹)x⁰ - (G(x) - f(0)x⁰)x¹ - G(x)x² = 0
G(x) - G(x)x - G(x)x² = x
G(x) = x/(1-x-x²)

Binet formula

       αⁿ - βⁿ    1  ⎛⎛  1+√5̅  ⎞ⁿ  ⎛  1-√5̅  ⎞ⁿ⎞
f(n) = ——————— = ——— ⎜⎜ —————— ⎟ - ⎜ —————— ⎟ ⎟
         √5̅      √5̅  ⎝⎝    2   ⎠   ⎝    2   ⎠ ⎠

  0x⁰ + 1x¹ + 1x² + 2x³ + 3x⁴ + 5x⁵ + 8x⁶ + ...

    x
= ——————
  1-x-x²

        x               quadratic factoring
= ————————————          ⎰ α = (1 + √5̅) / 2
  (1-αx)(1-βx)          ⎱ β = (1 - √5̅) / 2

  1/(α-β)   1/(β-α)
= ——————— + ———————     partial fraction expansion
   1-αx      1-βx  

   1    1     1    1 
= ——— ———— - ——— ————
  √5̅  1-αx   √5̅  1-βx

   1                                   1
= ——— ((αx)⁰ + (αx)¹ + (αx)² + ...) - ——— ((βx)⁰ + (βx)¹ + (βx)² + ...)
  √5̅                                  √5̅

  α⁰-β⁰      α¹-β¹      α²-β²
= ————— x⁰ + ————— x¹ + ————— x² + ...
   √5̅         √5̅         √5̅  
  ^^^^^      ^^^^^      ^^^^^
  f(0)       f(1)       f(2)

luogu P4451

Fibonacci number test

golden ratio

         f(n)    1 + √5̅
φ = lim —————— = —————— = 1.618...
    n→∞ f(n-1)     2

f(n) = f(n-1) + f(n-2)
f(n) - f(n-1) - f(n-2) = 0
f(n)/f(n-1) - 1 - f(n-2)/f(n-1) = 0
φ - 1 - 1/φ = 0     when n → ∞
φ² - φ - 1 = 0
φ = (1 ± √5̅) / 2
φ = (1 + √5̅) / 2     since φ > 0

additive law

爬樓梯問題：每步一階或兩階，爬上第n階有幾種步法。

c(n) = c(n-1) + c(n-2)

爬上第n+m階，分成兩種情況：沒踩第n階、有踩第n階。

c(n+m) = c(n-1)c(m-1) + c(n)c(m)

爬樓梯問題的答案，右移一項，即是費氏數列。

f(n+m+1) = f(n)f(m) + f(n+1)f(m+1)

變數代換，n+1換成n。

f(n+m) = f(n-1)f(m) + f(n)f(m+1)

GCD distributive law

gcd(f(n), f(n-1)) = 1
fib(gcd(a,b)) = gcd(fib(a),fib(b))
f(n)|f(m) when n|m

Fibonacci base（Zeckendorf's theorem）

UVa 763 948

微分不動點/輾轉相除法

微分不動點是0數列。微分之後左移一項是倍增數列。微分之後左移兩項是Fibonacci數列。移位的傅立葉轉換是角加速度。

Pisano period

https://en.wikipedia.org/wiki/Pisano_period

luogu P4000

Benford's law

Fibonacci convolution

斐波那契卷積可以用來解決哥德巴赫猜想與黎曼猜想。我已找到一個精妙的證明，而硬碟沒有足夠的空位寫下。

Fibonacci sequence

Vajda's identity：f(n+i)f(n+j) - f(n)f(n+i+j) = (-1)ⁿf(i)f(j)

Catalan's identity：f(n)² - f(n-i)f(n+i) = (-1)ⁿ⁻ⁱf(i)

Cassini's identity：f(n-1)f(n+1) - f(n)² = (-1)ⁿ

Pell sequence

Lucas sequence

Ulam sequence

Skolem–Mahler–Lech theorem

linear differential equation

「一次微分方程式」。

變數從離散數列推廣成連續函數。

c₀ x(t) + c₁ d/dt x(t) + c₂ d²/dt² x(t) = 0

polynomial -> root -> partial fraction decomposition -> exp
https://physiology.weill.cornell.edu/faculty/skrabanek/lab/qbio/resources_2010/2010_4.2%20Stability%20and%20Linearization%20of%20ODEs.pdf

luogu P6613

finite difference

微分改成有限差分。微分方程式改成差分方程式。

備忘

polynomial sequence

discrete derivative

高階微分可以改變底數。

difference operator:
Δa(n) = a(n+1) - a(n)

Δᵏa(n) = sum (-1)ⁱ C(k,i) a(n+k-i) = sum (-1)ᵏ⁻ⁱ C(k,i) a(n+i) 
        i=0⋯k                       i=0⋯k

Δ(x³) = (x+1)³ - x³ = 3x² + 3x + 1     let a(x) = x³

Δᵏ(xⁿ) = sum (-1)ⁱ C(k,i) (x+k-i)ⁿ
        i=0⋯k

derangement number:
D(n) = sum (-1)ⁿ⁻ⁱ C(n,i) i!
      i=0⋯n

D(n) = Δⁿ(x!) at x = 0
     = Δⁿ(0!)

shift operator: E

continuous derivative

d/dx(xⁿ) = n xⁿ⁻¹

L'Hôpital's rule: a/b = a′/b′
https://en.wikipedia.org/wiki/L'Hôpital's_rule

Hasse–Teichmüller derivative

Dᵏ(xⁿ) = C(n,k) xⁿ⁻ᵏ
https://en.wikipedia.org/wiki/Hasse_derivative

Umbral derivative

Dᵏ(Bₙ(x)) = sum Bₙ₋ₖ(x) C(n,k) xᵏ
           k=0⋯n

Appell sequence        d/dx pₙ(x) = n pₙ₋₁(x)
Bernoulli polynomial   ΔBₙ(x) = n xⁿ⁻¹
Euler polynomial       ΔEₙ(x) = 2(xⁿ - Eₙ(x))

https://en.wikipedia.org/wiki/Appell_sequence
https://en.wikipedia.org/wiki/Bernoulli_polynomials#Another_explicit_formula
https://planetmath.org/eulerpolynomial
https://en.wikipedia.org/wiki/Umbral_calculus
https://blog.csdn.net/EI_Captain/article/details/118317406

https://en.wikipedia.org/wiki/Bernoulli_number
https://en.wikipedia.org/wiki/Bernoulli_polynomials
http://luschny.de/math/zeta/The-Bernoulli-Manifesto.html

polar derivative

Dₖ P(x) = n P(x) + (k-x) P′(x)

Apolar polynomials
⎰ p(x) = a₀x⁰ + ... + aₙxⁿ
⎱ q(x) = b₀x⁰ + ... + bₙxⁿ
where        a₀bₙ           a₁bₙ₋₁                aₙb₀
      (-1)⁰ —————— + (-1)¹ —————— + ... + (-1)ⁿ —————— = 0
            C(n,0)         C(n,1)               C(n,n)

https://mathoverflow.net/questions/378463/