Convolution

Convolution

繁中「摺積」、簡中「卷积」。繁中翻譯錯誤,下文使用卷積。

一、加性卷積:冪級數相乘(指數相加)。

二、乘性卷積:狄利克雷級數相乘(底數相乘)。

(2 -5  1)        ⟷  2x⁰ - 5x¹ + 1x²
    ∗                       ×
(1 -1  0)        ⟷  1x⁰ - 1x¹ + 0x²
    ‖                       ‖
(2 -7  7  1  0)  ⟷  2x⁰ - 7x¹ + 7x² + 1x³ + 0x⁴
(2 -5  1)        ⟷  2⋅1ˣ - 5⋅2ˣ + 1⋅3ˣ
    ∗                       ×
(1 -1  0)        ⟷  1⋅1ˣ - 1⋅2ˣ + 0⋅3ˣ
    ‖                       ‖
(2 -7  1  5  0   ⟷  2⋅1ˣ - 7⋅2ˣ + 1⋅3ˣ + 5⋅4ˣ + 0⋅5ˣ
 -1  0  0  0)       - 1⋅6ˣ + 0⋅7ˣ + 0⋅8ˣ + 0⋅9ˣ

Circular Convolution

「循環卷積」。數列頭尾循環。(當xⁿ循環、當nˣ循環)

整數系統、實數系統不可能出現循環卷積。餘數系統、複數系統才可能出現循環卷積。

(2 -5  1)  ⟷  2x⁰ - 5x¹ + 1x²
    ⊛                 ⊗
(1 -1  0)  ⟷  1x⁰ - 1x¹ + 0x²      (when x³ = x⁰, x⁴ = x¹, ...)
    ‖                 ‖
(3 -7  7)  ⟷  3x⁰ - 7x¹ + 7x²
(2 -5  1)  ⟷  2⋅1ˣ - 5⋅2ˣ + 1⋅3ˣ
    ⊛                 ⊗
(1 -1  0)  ⟷  1⋅1ˣ - 1⋅2ˣ + 0⋅3ˣ   (when 4ˣ = 1ˣ, 5ˣ = 2ˣ, ...)
    ‖                 ‖
(7 -7  0)  ⟷  7⋅1ˣ - 7⋅2ˣ + 0⋅3ˣ

Convolution數學式子

數學式子有許多種寫法。

additive convolution:
(a∗b)(n) = sum { a(i)b(j) | i+j=n }          已知索引值n,用加法湊出n。
(a∗b)(n) = sum a(i)b(n-i) = sum a(n-i)b(i)   自身i、減出來的n-i,湊一對。
multiplicative convolution:
(a∗b)(n) = sum { a(i)b(j) | i×j=n }          已知索引值n,用乘法湊出n。
(a∗b)(n) = sum a(d)b(n/d) = sum a(n/d)b(d)   自身d、除出來的n/d,湊一對。

卷積由一個集合、三個運算組成。

三個運算目前沒有正式學術名稱。名稱是我瞎掰的。

                           | set         | match  | merge | summate
---------------------------| ------------| -------| ------| -------
additive convolution       | integer     | +      | ×     | +      
multiplicative convolution | integer     | ×      | ×     | +      
---------------------------| ------------| -------| ------| -------
dyadic convolution         | integer     | ∧/∨/⊕ | ×     | +      
infimal convolution        | function    | +      | +     | inf    
Minkowski sum              | vector sets | ×      | +     | ∪      

Convolution演算法

兩兩相乘,通通相加,時間複雜度O(N²)。

利用循環卷積,加性卷積降為O(NlogN),乘性卷積仍是謎。

                           | algorithm
---------------------------| ----------------------------------------------
additive convolution       | Fourier transform / number theoretic transform
multiplicative convolution | divide and conquer
---------------------------| ----------------------------------------------
dyadic convolution         | Walsh-Hadamard transform
infimal convolution        | ?
Minkowski sum              | Fourier transform

Additive Convolution

Additive Convolution(Cauchy Product)

「加性卷積」、「柯西乘積」。等價於冪級數乘法。

(a₀ a₁ a₂)        ⟷  a₀x⁰ + a₁x¹ + a₂x²
    ∗                         ×
(b₀ b₁ b₂)        ⟷  b₀x⁰ + b₁x¹ + b₂x²
    ‖                         ‖
(c₀ c₁ c₂ c₃ c₄)  ⟷  c₀x⁰ + c₁x¹ + c₂x² + c₃x³ + c₄x⁴
(a∗b)(n) = sum a(i)b(j) = sum a(i)b(n-i) = sum a(n-i)b(i)
          i+j=n           i≤n              i≤n

本章重點歸納成一張圖。

數列函數轉換【尚無正式名稱】

一、數列=多項式。

(a₀ a₁ a₂)  ⟷  a₀x⁰ + a₁x¹ + a₂x²

二、數列加法=多項式加法。

(a₀ a₁ a₂) + (b₀ b₁ b₂) = (a₀+b₀ a₁+b₁ a₂+b₂)
        a₀ x⁰ +       a₁ x¹ +       a₂ x²
+)      b₀ x⁰ +       b₁ x¹ +       b₂ x²
—————————————————————————————————————————
   (a₀+b₀) x⁰ +  (a₁+b₁) x¹ +  (a₂+b₂) x²

三、數列卷積=多項式乘法。

(a₀ a₁ a₂) ∗ (b₀ b₁ b₂) = (c₀ c₁ c₂ c₃ c₄)

c₀ = a₀b₀
c₁ = a₀b₁ + a₁b₀
c₂ = a₀b₂ + a₁b₁ + a₂b₀
c₃ = a₁b₂ + a₂b₁
c₄ = a₂b₂
                         a₀ x⁰ +   a₁ x¹ +   a₂ x²
×)                       b₀ x⁰ +   b₁ x¹ +   b₂ x²
——————————————————————————————————————————————————
                       a₀b₂ x² + a₁b₂ x³ + a₂b₂ x⁴
             a₀b₁ x¹ + a₁b₁ x² + a₂b₁ x³
   a₀b₀ x⁰ + a₁b₀ x¹ + a₂b₀ x²
——————————————————————————————————————————————————
     c₀ x⁰ +   c₁ x¹ +   c₂ x² +   c₃ x³ +   c₄ x⁴

四、數列循環卷積=多項式循環乘法。(當xⁿ循環)

(a₀ a₁ a₂) ⊛ (b₀ b₁ b₂) = (c₀ c₁ c₂)

c₀ = a₁b₂ + a₂b₁ + a₀b₀
c₁ = a₂b₂ + a₀b₁ + a₁b₀
c₂ = a₀b₂ + a₁b₁ + a₂b₀
                         a₀ x⁰ +   a₁ x¹ +   a₂ x²
⊗)                       b₀ x⁰ +   b₁ x¹ +   b₂ x²
—————————————————————————————————————————————————————
                       a₀b₂ x² + a₁b₂ x³ + a₂b₂ x⁴
             a₀b₁ x¹ + a₁b₁ x² + a₂b₁ x³
   a₀b₀ x⁰ + a₁b₀ x¹ + a₂b₀ x²
—————————————————————————————————————————————————————
   a₁b₂ x⁰ + a₂b₂ x¹ + a₀b₂ x²   (when x³ = x⁰, x⁴ = x¹, ...)
   a₂b₁ x⁰ + a₀b₁ x¹ + a₁b₁ x²
   a₀b₀ x⁰ + a₁b₀ x¹ + a₂b₀ x²
—————————————————————————————————————————————————————
     c₀ x⁰ +   c₁ x¹ +   c₂ x²

係值轉換【尚無正式名稱】

一、N-1次多項式,N個係數變成N個函數值。

冪級數恰是多項式函數。藉由多項式內插,N個函數值可以還原成N個係數。請見本站文件「Polynomial Interpolation」。

             a₀x⁰ + a₁x¹ + a₂x²
(a₀ a₁ a₂)  <——————————————————>  (y₀ y₁ y₂)
                 (x₀ x₁ x₂)

二、係數加法=多項式加法=函數值加法。

             a₀x⁰ + a₁x¹ + a₂x²
(a₀ a₁ a₂)  <——————————————————>  (p₀ p₁ p₂)
    +        b₀x⁰ + b₁x¹ + b₂x²       +
(b₀ b₁ b₂)  <——————————————————>  (q₀ q₁ q₂)
    ‖        c₀x⁰ + c₁x¹ + c₂x²       ‖
(c₀ c₁ c₂)  <——————————————————>  (r₀ r₁ r₂)
                 (x₀ x₁ x₂)

三、係數卷積=多項式乘法=函數值乘法。

                   a₀x⁰ + a₁x¹ + a₂x²
(a₀ a₁ a₂)        <————————————————————————————————>  (p₀ p₁ p₂)
    ∗              b₀x⁰ + b₁x¹ + b₂x²                     ×
(b₀ b₁ b₂)        <————————————————————————————————>  (q₀ q₁ q₂)
    ‖              c₀x⁰ + c₁x¹ + c₂x² + c₃x³ + c₄x⁴       ‖
(c₀ c₁ c₂ c₃ c₄)  <————————————————————————————————>  (r₀ r₁ r₂)
                              (x₀ x₁ x₂)

四、係數循環卷積=多項式循環乘法=函數值乘法。(當xⁿ循環)

             a₀x⁰ + a₁x¹ + a₂x²
(a₀ a₁ a₂)  <——————————————————>  (p₀ p₁ p₂)
    ⊛        b₀x⁰ + b₁x¹ + b₂x²       ×
(b₀ b₁ b₂)  <——————————————————>  (q₀ q₁ q₂)   (when x³ = x⁰, x⁴ = x¹, ...)
    ‖        c₀x⁰ + c₁x¹ + c₂x²       ‖
(c₀ c₁ c₂)  <——————————————————>  (r₀ r₁ r₂)
                 (x₀ x₁ x₂)

係值轉換(矩陣版本)【尚無正式名稱】

一、係值轉換視作線性函數、寫做矩陣。

A = [ x⁰ x¹ x² ]   transform

    [ a₀ ]
a = [ a₁ ]         coefficients
    [ a₂ ]

                  [ a₀ ]
Aa = [ x⁰ x¹ x² ] [ a₁ ] = a₀x⁰ + a₁x¹ + a₂x²   function value
                  [ a₂ ]

二、輸入加法=輸出加法。

A = [ x⁰ x¹ x² ]

    [ a₀ ]       [ b₀ ]
a = [ a₁ ]   b = [ b₁ ]
    [ a₂ ]       [ b₂ ]

A(a + b) = (Aa) + (Ab)

三、輸入卷積=輸出乘法。(尺寸不符,需要補項。)

à = [ x⁰ x¹ x² | x³ x⁴ ]

    [ a₀ ]       [ b₀ ]       [ c₀ ]
ã = [ a₁ ]   b̃ = [ b₁ ]   c̃ = [ c₁ ]
    [ a₂ ]       [ b₂ ]       [ c₂ ]
    [————]       [————]       [————]
    [ 0  ]       [ 0  ]       [ c₃ ]
    [ 0  ]       [ 0  ]       [ c₄ ]

Ã(ã ∗ b̃) = (Ãã) × (Ãb̃)

四、輸入循環卷積=輸出乘法。

A = [ x⁰ x¹ x² ]

    [ a₀ ]       [ b₀ ]
a = [ a₁ ]   b = [ b₁ ]
    [ a₂ ]       [ b₂ ]

A(a ⊛ b) = (Aa) × (Ab)

指數一齊添上任意倍率,運算性質仍成立。

元素一齊添上任意倍率,運算性質仍成立。

一個橫條併成多個橫條,運算性質仍成立。

exponent scalability | scalability          | row concatenation
---------------------| ---------------------| ----------------------
A = [ x⁰ x¹ x² ]     | A = [ x⁰ x¹ x² ]     |     [ 7x⁰ 7x¹  7x²   ]
A = [ x⁰ x⁵ x¹⁰ ]    | A = [ 7x⁰ 7x¹ 7x² ]  | A = [ -x⁰ -x⁵  -x¹⁰  ]
A = [ x⁰ x⁻⁷ x⁻²¹ ]  | A = [ -x⁰ -x⁵ -x¹⁰ ] |     [  x⁰  x⁻⁷  x⁻²¹ ]
A = [ x⁰ x⁰ x⁰ ]     | A = [ 0 0 0 ]        |     [  x⁰  x⁰   x⁰   ]

對偶係值轉換(矩陣版本)【尚無正式名稱】

接下來繼續補強矩陣,既滿足「輸入循環卷積=輸出乘法」,也滿足「輸入乘法=輸出循環卷積」,形成對偶運算。

從數學來看,補強性質,達成了對稱之美。從計算學來看,追加限制,產生了特殊演算法。

令原矩陣A、反矩陣A⁻¹,同時具備運算性質。一種嘗試是正規正交矩陣A⁻¹ = Aᵀ:令矩陣對稱。

實數系統,xⁿ漸增,數值漸增,無法正規正交。失敗!

               [ x⁰ x⁰ x⁰ x⁰ .. ]
               [ x⁰ x¹ x² x³ .. ]
A = A⁻¹ = Aᵀ = [ x⁰ x² x⁴ x⁶ .. ]   ✘
               [ x⁰ x³ x⁶ x⁹ .. ]
               [ :  :  :  :     ]

進一步嘗試closed set:令xⁿ循環,令數值循環。

一、複數系統傅立葉轉換:令x是N次單位根(的倒數)。

二、餘數系統數論轉換:令x是N階單位根(的倒數)。

Fourier transform: x = e-𝑖(2π/N), Nth root of unity
number theoretic transform: x = root of unity with order N

    [ x⁰ x⁰ x⁰ x⁰ .. ]             [ x⁻⁰ x⁻⁰ x⁻⁰ x⁻⁰ .. ]
    [ x⁰ x¹ x² x³ .. ]          1  [ x⁻⁰ x⁻¹ x⁻² x⁻³ .. ]
A = [ x⁰ x² x⁴ x⁶ .. ]   A⁻¹ = ——— [ x⁻⁰ x⁻² x⁻⁴ x⁻⁶ .. ]
    [ x⁰ x³ x⁶ x⁹ .. ]          N  [ x⁻⁰ x⁻³ x⁻⁶ x⁻⁹ .. ]
    [ :  :  :  :     ]             [ :   :   :   :      ]

總結:傅立葉轉換

數列、函數,互相轉換!

正向轉換:頭腦體操O(1)。

反向轉換:頭腦體操O(1)。

(a₀ a₁ a₂)  ⟷  a₀x⁰ + a₁x¹ + a₂x²

係數、函數值,互相轉換!

正向轉換:多項式求值N回合O(N²)。

反向轉換:多項式內插O(N³)。

             a₀x⁰ + a₁x¹ + a₂x²
(a₀ a₁ a₂)  <——————————————————>  (y₀ y₁ y₂)
                 (x₀ x₁ x₂)

矩陣形式。Vandermonde Matrix。

正向轉換:矩陣求值O(N²)。

反向轉換:矩陣求解O(N³)。

[ x₀⁰ x₀¹ x₀² ] [ a₀ ]   [ y₀ ]
[ x₁⁰ x₁¹ x₁² ] [ a₁ ] = [ y₁ ]
[ x₂⁰ x₂¹ x₂² ] [ a₂ ]   [ y₂ ]

高速演算法:複數系統傅立葉轉換、餘數系統數論轉換。

令矩陣共軛對稱、令矩陣正規正交、令xᵢ是單位根的每種次方、令矩陣邊長是2的次方。

正向轉換:O(NlogN)。

反向轉換:O(NlogN)。

[ x₀⁰ x₀¹ x₀² x₀³ ] [ a₀ ]   [ y₀ ]
[ x₁⁰ x₁¹ x₁² x₁³ ] [ a₁ ] = [ y₁ ]
[ x₂⁰ x₂¹ x₂² x₂³ ] [ a₂ ]   [ y₂ ]
[ x₃⁰ x₃¹ x₃² x₃³ ] [ a₃ ]   [ y₃ ]

x = (ω⁰ ω¹ ... ωᴺ⁻¹)
ω = root of unity with order N

總結:循環卷積

係值轉換:係數循環卷積=函數值乘法。

傅立葉轉換:係數循環卷積=函數值乘法。係數乘法=函數值循環卷積。

循環卷積演算法:係數循環卷積=正向傅立葉轉換&函數值乘法&逆向傅立葉轉換。

Circular Convolution演算法
Polynomial Circular Multiplication演算法

數列循環卷積的高速演算法:O(NlogN)!

正向傅立葉轉換O(NlogN),對應項相乘O(N),逆向傅立葉轉換O(NlogN)。總時間複雜度O(NlogN)。

傅立葉轉換的弱點是記憶體變兩倍(複數)、浮點數誤差。數論轉換的弱點是數值只能是整數、大於等於零、小於模數。

數列長度必須是2的次方。當數列長度不是2的次方,千萬不能直接補零到2的次方。

循環卷積計算結果:
(a₀ a₁ a₂) ⊛ (b₀ b₁ b₂) = (c₀ c₁ c₂)

c₀ = a₀b₀ + a₁b₂ + a₂b₁
c₁ = a₀b₁ + a₁b₀ + a₂b₂
c₂ = a₀b₂ + a₁b₁ + a₂b₀

補零直到長度是2的次方,計算結果完全不對:
(a₀ a₁ a₂ 0) ⊛ (b₀ b₁ b₂ 0) = (d₀ d₁ d₂ d₃)

d₀ = a₀b₀ + a₂b₂
d₁ = a₀b₁ + a₁b₀
d₂ = a₀b₂ + a₁b₁ + a₂b₀
d₃ = a₁b₂ + a₂b₁

正確方式:先補零直到不受循環影響,再補零直到長度是2的次方,最後讓輸出數列循環。

想要計算
(a₀ a₁ a₂) ⊛ (b₀ b₁ b₂) = (c₀ c₁ c₂)
先補零直到不受循環影響
(a₀ a₁ a₂ 0 0) ⊛ (b₀ b₁ b₂ 0 0) = (d₀ d₁ d₂ d₃ d₄)
再補零直到長度是2的次方
(a₀ a₁ a₂ 0 0 0 0 0) ⊛ (b₀ b₁ b₂ 0 0 0 0 0) = (d₀ d₁ d₂ d₃ d₄ 0 0 0)
最後讓輸出數列循環
c₀ = d₀ + d₃
c₁ = d₁ + d₄
c₂ = d₂

Convolution演算法
Polynomial Multiplication演算法

運用循環卷積,計算卷積。總時間複雜度O(NlogN)。

想要計算
(a₀ a₁ a₂ a₃) ∗ (b₀ b₁ b₂) = (c₀ c₁ c₂ c₃ c₄ c₅)
先補零直到不受循環影響
(a₀ a₁ a₂ a₃ 0 0) ⊛ (b₀ b₁ b₂ 0 0 0) = (c₀ c₁ c₂ c₃ c₄ c₅)
再補零直到長度是2的次方
(a₀ a₁ a₂ a₃ 0 0 0 0) ⊛ (b₀ b₁ b₂ 0 0 0 0 0) = (c₀ c₁ c₂ c₃ c₄ c₅ 0 0)
截斷輸出數列至正確長度
(c₀ c₁ c₂ c₃ c₄ c₅ 0 0) → (c₀ c₁ c₂ c₃ c₄ c₅)

Additive Convolution兩種觀點

線性函數有兩種觀點:矩陣切成直條(向量們的加權總和)、矩陣切成橫條(許多次向量點積),請見本站文件「Basis」。

卷積視作線性函數,也有兩種觀點。

(a₀ a₁ a₂ a₃ a₄) ∗ (b₀ b₁ b₂) = (c₀ c₁ c₂ c₃ c₄ c₅ c₆)

[ a₀ 0  0  0  0  0  0  ] [ b₀ ]   [ c₀ ]     c₀ = a₀b₀
[ a₁ a₀ 0  0  0  0  0  ] [ b₁ ]   [ c₁ ]     c₁ = a₀b₁ + a₁b₀
[ a₂ a₁ a₀ 0  0  0  0  ] [ b₂ ]   [ c₂ ]     c₂ = a₀b₂ + a₁b₁ + a₂b₀
[ a₃ a₂ a₁ a₀ 0  0  0  ] [ 0  ] = [ c₃ ]     c₃ = a₁b₂ + a₂b₁ + a₃b₀
[ a₄ a₃ a₂ a₁ a₀ 0  0  ] [ 0  ]   [ c₄ ]     c₄ = a₂b₂ + a₃b₁ + a₄b₀
[ 0  a₄ a₃ a₂ a₁ a₀ 0  ] [ 0  ]   [ c₅ ]     c₅ = a₃b₂ + a₄b₁
[ 0  0  a₄ a₃ a₂ a₁ a₀ ] [ 0  ]   [ c₆ ]     c₆ = a₄b₂

一、矩陣切成直條:數列移位、數列加權總和。

此觀點出現於「Signal Interpolation」。

(a₀ a₁ a₂ a₃ a₄) ∗ (b₀ b₁ b₂) = b₀ ⋅ (a₀ a₁ a₂ a₃ a₄ 0  0 )
                              + b₁ ⋅ (0  a₀ a₁ a₂ a₃ a₄ 0 )
                              + b₂ ⋅ (0  0  a₀ a₁ a₂ a₃ a₄)

二、矩陣切成橫條:數列顛倒、數列移位、數列點積。

此觀點出現於「Filter」、「Image Filtering」。

c₀ = (a₀ 0  0  0  0  0  0 ) ∙ (b₀ b₁ b₂ 0  0)
c₁ = (a₁ a₀ 0  0  0  0  0 ) ∙ (b₀ b₁ b₂ 0  0)
c₂ = (a₂ a₁ a₀ 0  0  0  0 ) ∙ (b₀ b₁ b₂ 0  0)
c₃ = (a₃ a₂ a₁ a₀ 0  0  0 ) ∙ (b₀ b₁ b₂ 0  0)
c₄ = (a₄ a₃ a₂ a₁ a₀ 0  0 ) ∙ (b₀ b₁ b₂ 0  0)
c₅ = (0  a₄ a₃ a₂ a₁ a₀ 0 ) ∙ (b₀ b₁ b₂ 0  0)
c₆ = (0  0  a₄ a₃ a₂ a₁ a₀) ∙ (b₀ b₁ b₂ 0  0)

範例:大數乘法

大數乘法即是多項式乘法!

數論轉換、傅立葉轉換得以計算大數乘法,時間複雜度從O(N²)降為O(NlogN)。

範例:兩兩和(Minkowski Sum)(X+Y Problem)

甲集合、乙集合,只有整數。甲取一個數,乙取一個數,相加,會是哪些數?

集合資料結構是循序儲存:窮舉。O(N²)。

集合資料結構是索引儲存:卷積。O(RlogR)。R為數字範圍。

多項式觀點:整數是指數,該整數出現與否(未出現0,出現>0)是係數,兩兩和是多項式乘法。

卷積觀點:整數是索引值,相加是索引值位移,兩兩和是卷積。

範例:組合計數(Combination Counting)

蓮霧2顆、椪柑3顆、芭樂5顆。欲在盤中擺放6顆水果,有幾種方式?

(x⁰+x¹+x²)(x⁰+x¹+x²+x³)(x⁰+x¹+x²+x³+x⁴+x⁵),找到x⁶的係數。加是或、乘是且、次方是數量、係數是方式。如同兩兩和。

範例:數列搜尋(Sequence Searching)

經典問題「一條陣列,尋找一個值」,解法是循序搜尋、排序之後再二分搜尋。

進階問題「一條陣列,尋找一串連續數列」,解法是字串搜尋、卷積。此處討論卷積。

兩串數列 a b
a = (1 2 3 4)
b = (1 3 5 7)

定義兩串數列 a b 的差距:對應項的差的平方的總和。
(1 - 1)² + (2 - 3)² + (3 - 5)² + (4 - 7)²

定義兩串數列 a b 的差距:對應項的差的平方的總和。
sum (a[i] - b[i])²

當差距等於零,則兩串數列相同。
sum (a[i] - b[i])² = 0  ---->  a = b

展開之後,重新整理,以數列為主角:數列每項平方和、數列點積。
  sum (a[i] - b[i])²
= sum (a[i]² + b[i]² - 2 a[i] b[i])
= sum a[i]² + sum b[i]² - 2 sum a[i] b[i]

縮寫。
(a - b)² = a² + b² - 2 (a · b)
兩串數列
(1 2 3 4 5 6)
(3 4 5)

窮舉各種對齊方式,判斷是否相符。
(1 2 3 4 5 6)  (1 2 3 4 5 6)  (1 2 3 4 5 6)  (1 2 3 4 5 6)
 | | |            | | |            | | |            | | |
(3 4 5)          (3 4 5)          (3 4 5)          (3 4 5)

換句話說,判斷差距是否等於零。
        (a - b)²              a²         b²     2(a·b)
(1-3)² + (2-4)² + (3-5)² = 1²+2²+3² + 3²+4²+5² - 2×26 = 12
(2-3)² + (3-4)² + (4-5)² = 2²+3²+4² + 3²+4²+5² - 2×38 = 3
(3-3)² + (4-4)² + (5-5)² = 3²+4²+5² + 3²+4²+5² - 2×50 = 0  <--- here!
(4-3)² + (5-4)² + (6-5)² = 4²+5²+6² + 3²+4²+5² - 2×62 = 5
                           ^^^^^^^^                ^^
                           moving window           convolution

預先計算每個數字的平方、預先計算各種對齊方式的點積(就是卷積),就可以快速求得a² + b² - 2ab。時間複雜度O(NlogN)。

搜尋(找到零)、最相似(找最小值)、最相異(找最大值)。

UVa 12298 ICPC 4671 5705 7159

Additive Convolution: Fourier Transform

Fourier Transform

「傅立葉轉換」。輸入複數數列,輸出複數數列,長度相同。

fourier(x₀ x₁ x₂ x₃ x₄) = (y₀ y₁ y₂ y₃ y₄)

「傅立葉矩陣」。傅立葉轉換的變換矩陣。

一、元素:N次單位根(的倒數)的次方。

二、直條:N次單位根(的倒數)的每種次方。

三、每個直條:次方值乘上每種倍率。

正向轉換,使用單位根的倒數e-𝑖(2π/N);逆向轉換,使用單位根e+𝑖(2π/N)。大家顛倒使用,我不知道原因。

    [ ω⁰  ω⁰  ω⁰  ω⁰ .. ] [ x₀ ]   [ y₀ ]
 1  [ ω⁰  ω¹  ω²  ω³ .. ] [ x₁ ]   [ y₁ ]
——— [ ω⁰  ω²  ω⁴  ω⁶ .. ] [ x₂ ] = [ y₂ ]   where ωᴺ = 1, ω = e-𝑖(2π/N)
√N̅  [ ω⁰  ω³  ω⁶  ω⁹ .. ] [ x₃ ]   [ y₃ ]
    [ :   :   :   :     ] [ :  ]   [ :  ]
           Fourier transform

    [ ω⁻⁰ ω⁻⁰ ω⁻⁰ ω⁻⁰ .. ] [ y₀ ]   [ x₀ ]
 1  [ ω⁻⁰ ω⁻¹ ω⁻² ω⁻³ .. ] [ y₁ ]   [ x₁ ]
——— [ ω⁻⁰ ω⁻² ω⁻⁴ ω⁻⁶ .. ] [ y₂ ] = [ x₂ ]
√N̅  [ ω⁻⁰ ω⁻³ ω⁻⁶ ω⁻⁹ .. ] [ y₃ ]   [ x₃ ]
    [ :   :   :   :      ] [ :  ]   [ :  ]
       inverse Fourier transform

正規正交

正交:直條兩兩點積是0。

正規:第一直條長度是√N̅,其餘直條長度是0。全體除以√N̅,令第一直條長度是1,其餘直條長度仍是0。

        [ ω⁰  ω⁰  ω⁰  ω⁰ .. ]             [ ω⁻⁰ ω⁻⁰ ω⁻⁰ ω⁻⁰ .. ]
     1  [ ω⁰  ω¹  ω²  ω³ .. ]          1  [ ω⁻⁰ ω⁻¹ ω⁻² ω⁻³ .. ]
F = ——— [ ω⁰  ω²  ω⁴  ω⁶ .. ]   F⁻¹ = ——— [ ω⁻⁰ ω⁻² ω⁻⁴ ω⁻⁶ .. ]
    √N̅  [ ω⁰  ω³  ω⁶  ω⁹ .. ]         √N̅  [ ω⁻⁰ ω⁻³ ω⁻⁶ ω⁻⁹ .. ]
        [ :   :   :   :     ]             [ :   :   :   :      ]

√N̅計算費時。改成正向轉換不縮放,逆向轉換除以N。

    [ x⁰  x⁰  x⁰  x⁰ .. ]             [ x⁻⁰ x⁻⁰ x⁻⁰ x⁻⁰ .. ]
    [ x⁰  x¹  x²  x³ .. ]          1  [ x⁻⁰ x⁻¹ x⁻² x⁻³ .. ]
F = [ x⁰  x²  x⁴  x⁶ .. ]   F⁻¹ = ——— [ x⁻⁰ x⁻² x⁻⁴ x⁻⁶ .. ]
    [ x⁰  x³  x⁶  x⁹ .. ]          N  [ x⁻⁰ x⁻³ x⁻⁶ x⁻⁹ .. ]
    [ :   :   :   :     ]             [ :   :   :   :      ]

Fast Fourier Transform

直接計算,時間複雜度O(N²)。數列長度是2的次方,則有高速演算法,時間複雜度O(NlogN)。此時稱作「快速傅立葉轉換」。

演算法(Cooley-Tukey Algorithm)

原理是Dynamic Programming。

一、數列長度必須是2的次方,以便對半分。

此例當中,數列長度N = 8 = 2³,總共對半分3回合。

[ ω⁰  ω⁰  ω⁰  ω⁰  ω⁰  ω⁰  ω⁰  ω⁰  ] [ x₀ ]   [ y₀ ]
[ ω⁰  ω¹  ω²  ω³  ω⁴  ω⁵  ω⁶  ω⁷  ] [ x₁ ]   [ y₁ ]
[ ω⁰  ω²  ω⁴  ω⁶  ω⁸  ω¹⁰ ω¹² ω¹⁴ ] [ x₂ ]   [ y₂ ]
[ ω⁰  ω³  ω⁶  ω⁹  ω¹² ω¹⁵ ω¹⁸ ω²¹ ] [ x₃ ] = [ y₃ ]
[ ω⁰  ω⁴  ω⁸  ω¹² ω¹⁶ ω²⁰ ω²⁴ ω²⁸ ] [ x₄ ]   [ y₄ ]
[ ω⁰  ω⁵  ω¹⁰ ω¹⁵ ω²⁰ ω²⁵ ω³⁰ ω³⁵ ] [ x₅ ]   [ y₅ ]
[ ω⁰  ω⁶  ω¹² ω¹⁸ ω²⁴ ω³⁰ ω³⁶ ω⁴² ] [ x₆ ]   [ y₆ ]
[ ω⁰  ω⁷  ω¹⁴ ω²¹ ω²⁸ ω³⁵ ω⁴² ω⁴⁹ ] [ x₇ ]   [ y₇ ]

二、直條/橫條交換,輸入元素/輸出元素隨之交換。

直條交換,輸入分成偶數項與奇數項,形成區塊。

[ ω⁰  ω⁰  ω⁰  ω⁰  | ω⁰  ω⁰  ω⁰  ω⁰  ] [ x₀ ]   [ y₀ ]
[ ω⁰  ω²  ω⁴  ω⁶  | ω¹  ω³  ω⁵  ω⁷  ] [ x₂ ]   [ y₁ ]
[ ω⁰  ω⁴  ω⁸  ω¹² | ω²  ω⁶  ω¹⁰ ω¹⁴ ] [ x₄ ]   [ y₂ ]
[ ω⁰  ω⁶  ω¹² ω¹⁸ | ω³  ω⁹  ω¹⁵ ω²¹ ] [ x₆ ]   [ y₃ ]
[———————————————————————————————————] [————] = [————]
[ ω⁰  ω⁸  ω¹⁶ ω²⁴ | ω⁴  ω¹² ω²⁰ ω²⁸ ] [ x₁ ]   [ y₄ ]
[ ω⁰  ω¹⁰ ω²⁰ ω³⁰ | ω⁵  ω¹⁵ ω²⁵ ω³⁵ ] [ x₃ ]   [ y₅ ]
[ ω⁰  ω¹² ω²⁴ ω³⁶ | ω⁶  ω¹⁸ ω³⁰ ω⁴² ] [ x₅ ]   [ y₆ ]
[ ω⁰  ω¹⁴ ω²⁸ ω⁴² | ω⁷  ω²¹ ω³⁵ ω⁴⁹ ] [ x₇ ]   [ y₇ ]

三、區塊運算,一道式子變成上下兩道式子。

左區塊是傅立葉轉換(單位根翻倍)。ω⁸ ≡ ω⁰ ≡ 1。

[ ω⁰  ω⁰  ω⁰  ω⁰  ] [ x₀ ]   [ ω⁰  ω⁰  ω⁰  ω⁰  ] [ x₁ ]   [ y₀ ]
[ ω⁰  ω²  ω⁴  ω⁶  ] [ x₂ ] + [ ω¹  ω³  ω⁵  ω⁷  ] [ x₃ ] = [ y₁ ]
[ ω⁰  ω⁴  ω⁸  ω¹² ] [ x₄ ]   [ ω²  ω⁶  ω¹⁰ ω¹⁴ ] [ x₅ ]   [ y₂ ]
[ ω⁰  ω⁶  ω¹² ω¹⁸ ] [ x₆ ]   [ ω³  ω⁹  ω¹⁵ ω²¹ ] [ x₇ ]   [ y₃ ]

[ ω⁰  ω⁸  ω¹⁶ ω²⁴ ] [ x₀ ]   [ ω⁴  ω¹² ω²⁰ ω²⁸ ] [ x₁ ]   [ y₄ ]
[ ω⁰  ω¹⁰ ω²⁰ ω³⁰ ] [ x₂ ] + [ ω⁵  ω¹⁵ ω²⁵ ω³⁵ ] [ x₃ ] = [ y₅ ]
[ ω⁰  ω¹² ω²⁴ ω³⁶ ] [ x₄ ]   [ ω⁶  ω¹⁸ ω³⁰ ω⁴² ] [ x₅ ]   [ y₆ ]
[ ω⁰  ω¹⁴ ω²⁸ ω⁴² ] [ x₆ ]   [ ω⁷  ω²¹ ω³⁵ ω⁴⁹ ] [ x₇ ]   [ y₇ ]

四、直條/橫條除以倍率,輸入元素/輸出元素隨之乘以倍率。

右區塊是傅立葉轉換。橫條分別除以ωⁿ,輸出分別乘以ωⁿ。

fft(x₀ x₂ x₄ x₆) + fft(x₁ x₃ x₅ x₇) × (ω⁰ ω¹ ω² ω³) = (y₀ y₁ y₂ y₃)
fft(x₀ x₂ x₄ x₆) + fft(x₁ x₃ x₅ x₇) × (ω⁴ ω⁵ ω⁶ ω⁷) = (y₄ y₅ y₆ y₇)
fft(xeven) + fft(xodd) × ωⁿlow  = ylow 
fft(xeven) + fft(xodd) × ωⁿhigh = yhigh

遞推過程

偶數項與奇數項分開遞歸,形成網路。

重新排列

偶數項與奇數項分開遞歸,索引值不連續。

預先重新排列,減少cache miss,達成inplace。

索引值視作二進位,高低位數顛倒。時間複雜度O(NlogN)。

建立表格。時間複雜度降為O(N),空間複雜度升為O(N)。

預先建立表格。重複執行傅立葉轉換的情況下,節省時間。

顛倒過程

顛倒過程,即是逆向轉換:

一、逆向公式(此例是最後一個步驟,位於網路右下方)。

   { x̕₆ + x̕₇⋅ω³ = y₃  →  { x̕₆ + x̕₇⋅ω³ = y₃     where ω⁴ = -1
   { x̕₆ + x̕₇⋅ω⁷ = y₇     { x̕₆ - x̕₇⋅ω³ = y₇

→  { y₃ + y₇ = 2⋅x̕₆     →  { (y₃ + y₇) / 2 = x̕₆
   { y₃ - y₇ = 2⋅x̕₇⋅ω³     { (y₃ - y₇) / ω³ / 2 = x̕₇

二、每回合除以2,改成最後一口氣除以N。節省時間。

三、除以ωⁿ,改成乘以ω⁻ⁿ。節省時間。

演算法(Gentleman-Sande Algorithm)

顛倒過程,亦可求得正向轉換:

一、正向轉換、逆向轉換,單位根ω互為倒數。

二、另一個差別:正向轉換不縮放,逆向轉換除以N。

三、逆向轉換,單位根倒數、乘以N,就是正向轉換。

四、正向轉換,顛倒過程、單位根倒數、乘以N,仍是正向轉換。

五、最後一口氣除以N、乘以N,兩者相消。

六、ωⁿ每回合開方,改成每回合翻倍。不必預先建立表格!

程式碼

Additive Convolution: Number Theoretic Transform

Number Theoretic Transform

「數論轉換」。輸入餘數數列,輸出餘數數列,長度相同。

ntt(x₀ x₁ x₂ x₃ x₄) ≡ (y₀ y₁ y₂ y₃ y₄) (mod m)
[ ω⁰  ω⁰  ω⁰  ω⁰  ω⁰  ] [ x₀ ]   [ y₀ ]
[ ω⁰  ω¹  ω²  ω³  ω⁴  ] [ x₁ ]   [ y₁ ]
[ ω⁰  ω²  ω⁴  ω⁶  ω⁸  ] [ x₂ ] ≡ [ y₂ ]   where ω⁵ ≡ 1 (mod m)
[ ω⁰  ω³  ω⁶  ω⁹  ω¹² ] [ x₃ ]   [ y₃ ]
[ ω⁰  ω⁴  ω⁸  ω¹² ω¹⁶ ] [ x₄ ]   [ y₄ ]

選擇模數

輸入輸出的每個數值都是餘數,大於等於零、小於模數。

當輸入輸出數值很大,那麼模數必須足夠大。

一、數列長度N,單位根ωᴺ ≡ 1,單位根數值皆異。

二、令次方值的模數是N的適當倍數,ord(m) = kN。

三、令模數是質數,根據費馬小定理,m = kN + 1。

四、令模數是質數,存在原根g。

五、令單位根ω ≡ gᵏ,使得ωᴺ ≡ gᵏᴺ ≡ 1。

程式碼

Dyadic Convolution: Walsh-Hadamard Transform(Under Construction!)

Walsh-Hadamard Transform

N=2的傅立葉轉換,推廣到高維度。

矩陣形式是Hadamard matrix或Walsh matrix。數值都是±1而且向量互相垂直(點積為零)。

因為沒有循環,所以沒有對偶?

http://iris.elf.stuba.sk/JEEEC/data/pdf/09-10_102-10.pdf
http://www.jjj.de/fxt/fxtbook.pdf   p457  p481 $23.8
http://en.wikipedia.org/wiki/Hadamard_matrix

Multiplicative Convolution

Multiplicative Convolution(Dirichlet Convolution)

「乘性卷積」、「狄利克雷卷積」。等價於狄利克雷級數乘法。

(a₀ a₁ a₂)         ⟷  a₀1ˣ + a₁2ˣ + a₂3ˣ
    ∗                          ×
(b₀ b₁ b₂)         ⟷  b₀1ˣ + b₁2ˣ + b₂3ˣ
    ‖                          ‖
(c₀ c₁ c₂ ... c₈)  ⟷  c₀0ˣ + c₁1ˣ + c₂2ˣ + ... + c₂9ˣ
(a∗b)(n) = sum a(i)b(j) = sum a(d)b(n/d) = sum a(n/d)b(d)
          i×j=n           d|n              d|n

Multiplicative Convolution: Abel's Summation Formula(Under Construction!)

Abel's Summation Formula

演算法是Divide and Conquer。計算一項的時間複雜度O(n¾)。乘性函數降至O(n)。

https://codeforces.com/blog/entry/54150

luogu P4213

Ramanujan's Sum(Under Construction!)

Fourier Ramanujan Transform

[totient multiplicative]
φ(mn) = φ(m)φ(n)
與m互質的數a  與n互質的數b 與mn互質的數x
{ x = a (mod m)     m個m個排整齊剛好是一堆column
{ x = b (mod n)     n個n個排整齊剛好是一堆column
窮舉a與b實施中國餘數定理得到x

[totient divisor sum]
sum φ(d) = n
d|n
1/20 2/20 .... 20/20
約分之後,分母只可能是因數
分母是20的有φ(20)個
分母是10的有φ(10)個
分母是5的有φ(5)個
分母是4 2 1 有 φ(4) φ(2) φ(1)個

[totient gcd]
 sum gcd(i,n) = sum φ(n) * n/d  原理同[euler divisor sum]  n/d是gcd的值
i=1~n           d|n
gcd數列前綴和   φ(n)數列只取因數項

[totient 排容]
φ(n) = sum μ(d) * n/d = n sum μ(d)/d
       d|n                d|n
μ(d) = n質因數分解  質因數各一個相乘  零個偶數個相乘+1 奇數個相乘-1  其他0
順帶一提,隨便一個poset都可以定義排容,μ(n)是定在因數關係的poset。

[euler product]
http://en.wikipedia.org/wiki/Proof_of_the_Euler_product_formula_for_the_Riemann_zeta_function

(=>)
質因數分解/算術基本定理的倒數版本
每一個質數的各種次方,利用多項式相乘,拼出所有數

  2的次方             3的次方              5的次方
  (1+1/2+1/4+1/8+...)*(1+1/3+1/9+1/27+...)*(1+1/5+1/25+1/125+...)*...
= 1+1/2+1/3+1/4+1/5+1/6+....
  全部的數

(<=)
篩法! wiki有證明,無窮級數運算
  1+1/2+1/3+1/4+1/5+1/6+....
= (1+1/2+1/4+1/8+...)*(1+1/3+1/9+1/27+...)*(1+1/5+1/25+1/125+...)*...

(推廣)
次方值一齊乘上相同倍率!

(不是倒數的版本)
次方值s代入-1
(1+2+4+8+...)*(1+3+9+27+...)*(1+5+25+125+...)*... = 1+2+3+4+5+6+... = -1/12

[Dirichlet convolution]
sum f(d) * g(n/d)
d|n

[數列變成多項式]
              c[6] = a[0]b[6] + a[1]b[5] + a[2]b[4] + ... + a[6][0]
正常卷積,基於加法分解,故補入x的各個次方 x^1  x^2  x^3 ...
                 可推廣 x^1a x^2a x^3a
              c[6] = a[1]b[6] + a[2]b[3] + a[3]b[2] + a[6]b[1]
滴哩卷積,基於乘法分解,故補入各個數倒數  1/1  1/2  1/3 ... (近似euler product)
                 可推廣 1^a  2^a  3^a ...
                        或者是想成開方值 1次單位根 2次單位根 3次單位根???
用倒數而不用自然數,可能是想分項,讓各項不會混在一起吧。
只有一個row而不是一個矩陣
複數系統(z-transform)  http://en.wikipedia.org/wiki/Dirichlet_series
餘數系統(數論轉換row)  http://en.wikipedia.org/wiki/Bell_series
還沒找到代入單位圓/循環波的相關式子 :I

不知道複數系統下有沒有矩陣滿足卷積/乘法對稱性?
Hibert matrix只有一點點點像,但不是。

[滴哩卷積的單位元素/脈衝函數]
delta(n) = sum μ(d) = { 1  when n = 1
           d|n        { 0  otherwise

μ(n)是 n做質因數分解  質因數各一個相乘  零個偶數個相乘+1 奇數個相乘-1  其他0

假設n有m個質因數,觀察n的每一個因數d
次方超過1的,μ(d)都是0,加起來沒影響
次方是0和1的,總共有2^m種
零個偶數個、奇數個,一樣多:
針對一種質因數的取捨方式,最後一個質因數從有切換成沒有,正負號就變了。

[滴哩卷積的倒數]
令c(n)=1,補到原本的等式裡面,弄成卷積的樣子
sum μ(d) * c(n/d) = { 1 when n = 1    [mobius divisor]
d|n                 { 0 otherwise
sum φ(d) * c(n/d) = n               [totient divisor]
d|n
sum μ(d) * n/d = φ(n)              [totient 排容]
d|n
https://math.stackexchange.com/questions/135351/

μ(n)和c(n)=1的卷積是脈衝函數。μ(n)的倒數是水平線c(n)=1。非常奇葩。
φ(n)和c(n)=1的卷積是45度斜線id(n)=n。也很奇葩。
μ(n)和id(n)=n的卷積是φ(n)。

 [1;31m簡單來說就是因數關係poset的排容原理,可以變成滴哩卷積! [m
 [1;31m凡是排容原理,也許都可以變成某種多項式乘法,變成某種卷積。 [m

[滴哩卷積的計算]
http://en.wikipedia.org/wiki/Multiplicative_function
  id   id
μ -> φ -> gcd前綴和

設函數id(n)=n
 sum gcd(i,n) = φ(n) 滴哩卷積 id(n)  [euler gcd]
i=1~n
 sum f(gcd(i,n)) = φ(n) 滴哩卷積 f(n)
i=1~n

 [1;31m要如何簡單表示id的倒數?這樣就可以很快用gcd前綴和來算出φ(n) [m
http://math.stackexchange.com/questions/423317/
http://oeis.org/A023900
φ^-1 = core * μ
φ    = id   * μ
delta = core * id * μ * μ
id^-1 = core * μ * μ       好醜

[Mobius inverse]
f convolution 1               = g
f convolution 1 conovlution μ = g convolution μ
f                             = g convolution μ
阿就這樣

[拉馬努金和]
http://en.wikipedia.org/wiki/Ramanujan%27s_sum
                 nk
eta_q(n) = sum  w      n個複數波,n是倍率,q是取樣數,k是索引值,通通加起來
          k=1~q  q      首先定參一下取樣數是多少

         = { 0  if q!|n 宛如餘數加法循環,各數(各等分角)剛好出現一遍,加起來是0
           { q  if q|n  取樣通通抽中0度角=1。

                nk
c_q(n) = sum   w        針對每一種波,取樣時取不可約分(互質)的那些點,φ(q)個
       (k,q)=1  q

         if q is prime  [totient φ(q)=q-1 if q is prime 的推廣]
       = { -1   if q!|n 全部的取樣和 - 可以約分的取樣和 = 不可約分的取樣和
         { q-1  if q|n      0或q         1(只有第n點)         -1或q-1

eta_q(n) = sum c_d(n)   [totient divisor sum]的推廣,從計數推廣成複數單位圓
           d|q          (需要仔細研究一下怎麼分配位置,因為循環所以可以位移吧)
sum eta_d(n) * μ(q/d) = c_q(n)  卷積倒數
d|q
sum d        * μ(q/d) = c_q(n)  照eta的公式,等於0的就刪掉,剩下有因數的
d|gcd(q,n)
sum μ(q/d)   * d      = c_q(n)  調整一下寫法
sum μ(d)     * q/d    = c_q(n)  調整一下寫法
                 [1;31m是因數的,滴哩卷積/排容後,得到互質的 [m

φ(q) = c_q(0)           平坦波,取樣必定都是0度角=1。照定義來,取小於n的互質數
μ(q) = c_q(1)           一倍波。這很奇葩,需要想一下。
c_q(0) = sum c_d(1) * q/d   [totient 排容]
         d|q
c_q(n) = sum c_d(1) * q/d   前面推導的式子 (n=0是特例,剛好是totient排容)
        d|gcd(q,n)          這可以DP:比較小的取樣結果,遞推比較大的取樣結果

[傅立葉轉換]
滴哩卷積下,μ(n)和φ(d)和g(n)=1和g(n)=n經過某種傅立葉轉換之後不知道是啥鬼。

循環卷積下,gcd數列dft的第一項(不是第零項)是φ(n)
      gcd數列與一倍頻率複數波的內積是φ(n)
不知道怎麼證

[傅立葉轉換,針對特例gcd]
http://www.kurims.kyoto-u.ac.jp/EMIS/journals/INTEGERS/papers/i50/i50.pdf
注意到沒有第零項
c_o(1)滴哩卷積id(1) ... c_o(n)滴哩卷積id(n) <===> gcd(1,n) ... gcd(n,n)
      φ(n)

正向傅立葉轉換

???

逆向傅立葉轉換
                     ij
gcd(i,n) = 1/n sum (w  * sum c_d(j) * n/d )  n是數量,i是索引值
              j=1~n      d|n

                 ij
         = sum (w  * sum c_d(j)/d )  外面1/n跟裡面n相消
          j=1~n      d|n

注意到j=1的那一項,sum c_d(1) * n/d = c_d(0) = φ(n)
                   d|n

還可以算複數系統的最大公因數???

[拉馬努金傅立葉轉換]
http://perso.ens-lyon.fr/patrick.flandrin/NewOrl2.pdf

把基底從複數波換成拉馬努金和
正規化的norm是φ(n)
逆轉換需要定義一下"無限長不規則數列"的內積

Fourier Ramanujan Transform

http://perso.ens-lyon.fr/patrick.flandrin/NewOrl2.pdf

複數波改成Ramanujan Sum。

norm從N改成φ(n)。垂直是互質。逆轉換需要定義無限長數列的內積。

Fourier Transform
       N-1
X(n) =  ∑ { x(k) ÷ e𝑖⋅2π⋅(k/N)⋅n }
       k=0
Fourier Ramanujan Transform
       N          q
X(n) =  ∑ { x(q) ⋅ ∑ e𝑖⋅2π⋅(k/q)⋅n }
       q=1        k=1,k⊥q

最大公因數的傅立葉轉換就是Euler's Totient Function

https://zhuanlan.zhihu.com/p/166530236

       N-1
φ(n) =  ∑ { gcd(n, k) ÷ e𝑖⋅2π⋅(k/N)⋅n }
       k=0