二次篩選法

二次篩選（Quadratic Sieve）演演算法是一個整數分解演演算法，在實際用途中為已知第二快的方法（目前第一快為普通數域篩選法）。但對於大約 100 位數以內的整數，它仍然是最快的演算法，而且比起普通數域篩選法來說簡潔得多。

這是一個通用的整數分解演演算法，意即其運算時間完全取決於欲分解的整數本身位數的大小，而不是在於特殊結構或特性。

二次篩選法是由卡爾·帕梅朗斯在1981年所發明，並作為理察·施羅佩爾的線性篩法之改良版。

## 基礎目標

此演演算法試圖去建立一個模 $n$ （$n$為欲分解的數）下的平方同餘，這往往即是$n$的因數分解。演演算法有兩個階段：「數據收集」，在此階段收集可能可以找到一個平方同餘的資料；以及「數據處理」，它把所有收集的數據放進一個矩陣裡，並解決、獲得一個平方同餘數。

數據收集的階段可以很輕易地使用多個處理器去平行化。但數據處理階段需要大量的記憶體，並且在多個運算節點之間有效地平行化相當困難，也可能每個運算節點的記憶體不夠足以儲存整個陣列。而Block Wiedemann演演算法可以使用在一些可以保存陣列的系統。

要找到一個平方同餘，一個較天然的方法便是隨機挑選數字，將其平方，並希望模$n$之後的非負餘數是一個完全平方數。例如：$80^2$ 模 $5959 = 441$，同時也是$21^2$。

這種方法對很大的$n$值而言，可以找到一個同餘的平方數的情況很罕見。但是當真的找到了一個時，在大多數情況下，同餘數為非平凡解而整數分解便完成了。這大致上即是費馬因式分解法（Fermat's factorization method）的核心。

而二次篩選法改良自狄克森因式分解法。

一般來說，二次篩選法的執行時間（去質數分解一個整數$n$時）為

$$e^{(1 + o(1))\sqrt{\ln n \ln\ln n}} =L_n\left[1/2,1\right]$$

參見 L-符號。

上式常數 e 為自然對數之底數。

## 解決方法

令 $x$ 模 $y$ 表示為 $x$ 除以 $y$ 之後所剩的餘數。為了分解整數，費馬因式分解法牽涉到需要尋找一個數字 $a$（$n^{1/2}<a<n-1$），使得$a^2\pmod{n}$是一個完全平方數。但這些$a$值相當難找到。二次篩選法包括對於好幾個$a$值去計算了$a^2\pmod{n}$，然後在$a$值與$a^2\pmod{n}$的集合中找到一個子集，當中的元素之乘積為完全平方數。而產生出一個平方同餘。

例如：$41^2$ 模 $1649 = 32$、$42^2$ 模 $1649=115$ 以及 $43^2$ 模 $1649$ 為 $200$。在這些數字（32、115、200）當中皆無完全平方數，但存在一乘積$32\times 200=6400=80^2$是一個平方數。模1649 之後，這個乘積$32\times 200=(41^2)\times(43^2)=(41\times 43)^2=114^2$（因為 $41\times 43$ 模 $1649 = 114$）。$32\times 200 = 80 ^2$ 的觀察因而給出了一個平方同餘：$114^2\equiv 80^2$ (模 $1649$)。

但是，如何將以下問題解決呢？「給予一組數字，找到一個子集使其乘積是平方數。」該解決方案使用了指數向量的概念。而指數向量，例如根據算術基本定理，504 可分解為 $2^3 3^2 5^0 7^1$。

這表示可以藉由指數向量 $(3,2,0,1)$，代表 $2,3,5,7$ 在因式分解的指數值。490 會同樣可分解為指數向量 $(1,0,1,2)$。將這些數字相乘相當餘把其指數向量的對應值一一相加：$504\times 490$ 得一向量 $4,2,1,3$。

有一些數字為平方數，其滿足每個在其指數向量的各個數字為偶數。例如，向量 $(3,0,0,1)$、$(1,2,0,1)$ 之和為 $(4,2,0,2)$，因此 56 乘以 126 是一個平方數。找尋一個平方數只需對於向量裡數字之的奇偶性之知識，所以有可能將整個向量簡化為模 2 的形式並作模 2 下的加法: $(1,0,0,1)+(1,0,0,1)=(0,0,0,0)$。在實作中，這相當地有效率，因其可以表示為一位元集（bitsets）且模 2 之加法將變為位元運算互斥或（XOR）。

於是此問題變化為：「給予一個 0, 1向量構成的集合，找到一個子集，其中所有向量之和為模 2 的零向量。」而這是一個線性代數的問題；且該解答為線性相依的。

線性相依是線性代數中的一個定理：當有比每個向量中含有的元素還要多的向量時，這種相依關係必然存在。而它可以被高效率地找到，例如：把所有向量一列一列地排在一個矩陣裡，然後使用高斯消去法。比起實數來說，此方法尤其容易套用到模 2 後的整數上。而此演演算法所需的平方數即是那些向量所對應的數字之積。

然而，純粹地只去將一堆隨機數字平方並模 $n$ 會產生很大量的、不同的質因數，也因此會產生出很長的向量以及一個非常大的矩陣。解決方法是去找到一些特別的數字 $a$，使得 $a^2$ 模 $n$ 之值只由很小的質因數組成（它們都是光滑數）。此種數字很難被找到，但是若僅使用光滑數將可以保持向量和矩陣之尺寸更小、更容易處理。

而二次篩選法使用一種之後會提及名為「篩法」（sieving）的技巧去找尋光滑數，也就是此演演算法的命名由來。

## 演演算法

總地來說，二次篩演算法基本有以下主要的步驟：

- 選擇一個光滑數之上界 $B$。以質數計數函數 $\pi(B)$ 表示小於 $B$ 的質數之數量，其將控制之後向量的長度以及需要的向量之數量。
- 使用篩法找到 $\pi(B)+1$ 個數字 $a_i$ 使得 這樣 $b_i=(a_i^2$ 模 $n)$ 為一個 $B$-光滑數。
- 將 $b_i$ 作質因數分解生成一個指數向量（每個數字都要模 $2$）。
- 套用線性代數的概念找到一個子集，其中的每個向量之和為一零向量。把這些向量所對應的 $a_i$ 相乘、$b_i$ 相乘並模 $n$：便得到一個 $B$-光滑數 $b^2$.
- 現在，我們得到方程式 $a^2=b^2 ($ 模 $n)$因為從步驟 4 得到 $(a^2$ 模 $n)$ 的兩個平方根，一個即是對整數 $b^2$ 取平方根，也就是 $b$；另一個則是步驟 4 得到的 $a$ 本身。
- 因此現在，我們掌握了所需的恆等式：$(a+b)(a-b)\equiv0 \pmod n$ $(a+b)(a-b)\equiv0 \pmod n$ $(a+b)(a-b)\equiv0\pmod{n}$。計算 $n$ 與 $a-b$ （或是 $a+b$）的 GCD （最大公因數，Greatest Common Divisor）。此將產生 $n$ 的其中一個因數，儘管有可能是一個平凡（trivial）因數 (即為 $n$ 或是 1)。如果此因數是平凡的，使用不同的線性相依或是不同的 $a$ 值去再次嘗試。

本文的剩餘部分將解釋這個基本演演算法的細節和延伸。

## 二次篩法（QS）如何最佳化找尋同餘

二次篩選法試圖找到一整數對 $x$ 和 $y(x)$ （其中 $y(x)$ 為 $x$ 的函數）其滿足比 $x^2\equiv y ($模 $n)$ 還要弱得多的條件。它選擇一些質數作為一集合作為「因數基底」，並試圖找到 $x$，使得 $y(x)=x^2$ 模 $n$ 之值的質因數只會在此因數基底。此時可稱 $y$ 值：對於此因數基底是光滑的。

$y(x)$ 的其中一個值之因式分解（為因數基底的一部分），跟 $x$ 一起，被稱為「關係」（relation）。二次篩選法藉由採取接近 $n$ 的平方根之 $x$ 值，以加速尋找這類「關係」的過程。這將確保 $y(x)$ 會較小，因而具有更大的可能性是光滑的。

$$y(x)=\left(\left\lceil\sqrt{n}\right\rceil+x\right)^2-n\hbox{ (where }x\hbox{ is a small integer)}$$

$$y(x)\approx 2x\left\lceil\sqrt{n}\right\rceil$$

這意味著， $y$ 在 $2x\sqrt{n}$ 的數量級上.。然而，這也意味著 $y$ 的增長幅度與 $x$ 乘以 ($n$的平方根) 成正比。

另一個可以增加光滑的可能性是，即是單純地增大因數基底的大小。然而，比起因數基底的質數數量，至少找到一個光滑的「關係」還是必要許多，其確保存在一個線性相依。

### 「部分關係」以及循環

即使對於某些「關係」來說， $y(x)$ 並非光滑的。但如果兩個 $y$ 剛好是由因數基底以外的相同質數之乘積，也可能可以合併這兩個部分「關係」，以形成一個完整的「關係」。註：此形同於因數基底的擴展。

例如：如果因數基底為$\{2,3,5,7\}$和 $n=01$，存在「部分關係」（partial relations）：

$${21^2\equiv 7^1\cdot 11\pmod{91}}$$

$${29^2\equiv 2^1\cdot 11\pmod{91}}$$

將上面兩式乘在一起：

$${(21\cdot 29)^2\equiv2^1\cdot7^1\cdot11^2\pmod{91}}$$

並將等號兩邊皆乘上 $\left(11^{-1}\right)^2$ 模 $91$。而 $11^{-1}$ 對 $91$ 取模為 $58$，所以：

$$(58\cdot 21\cdot 29)^2\equiv 2^1\cdot7^1\pmod{91}$$

$$14^2\equiv 2^1\cdot7^1\pmod{91}$$

即產生了一個完整的「關係」。這樣的一個完整的「關係」（藉由結合「部分關係」所獲得的）稱為循環。有時候，從兩個「部份關係」形成的循環，可以直接導向一個平方同餘，但是此情況非常罕見。

### 藉由篩選來檢查光滑度

有好幾種方法可以 $y$ 值們的光滑度。最直覺的是藉由試除法，儘管這樣會增加數據收集階段的運行時間。

另一個方法較能被接受的方法是橢圓曲線因式分解（ECM）。

而在實作中，稱為篩選的方法比較會被經常使用。設 $f(x)$ 為多項式 $f(x)=x^2-n$ 我們得：

$$f(x)=x^2-n$$

$$f(x+kp)=(x+kp)^2-n$$

$$f(x+kp)=x^2+2xkp+(kp)^2-n$$

$$f(x+kp)=f(x)+2xkp+(kp)^2\equiv f(x)\pmod{p}$$

因此解決出 $f(x)\equiv 0 ($模 $p)$對於某個 $x$ 值，將產生出一整個序列，當中的每個數值 $y~(y=f(x)$ 皆可被 $p$ 整除。此問題便是對某個質數取模下找到一個平方根，對其存在著高效率的演演算法，例如謝克斯–托內里演演算法的。（這便是二次篩法的名稱來由： $y$ 是一個 $x$ 的二次多項式且篩選過程中的運算類似埃拉托斯特尼篩法。）

篩選一開始將一個大陣列 $A[]$ 每個「元」（entry）的每個位元組設為零。對於每一個 $p$，去解決模$p$ 下的二次方程式並得到兩個根 $\alpha$ 和 $\beta$，然後在每個 $y(x)=0$ 模 $p$「元」之中加入一個近似於 $\log(p)$ 之值……也就是 $A[kp+\alpha]$ 和 $A[kp+\beta]$。為了辨識數字是否可被因數基底中的質數之平方所整除，解決幾個模 ( $p$ 的小次方) 下的二次方程式也是必要的。

在因數基底的尾端，任何 $A[]$ 有包含一個值超過大約為 $\log(x^2-n)$ 的臨界值，將會對應到一個 $y(x)$ 值，其由因數基底的部分組成。那些包含了確定 $y(x)$ 可以被哪些質數整除的資訊已經遺失掉了，但是因為其只包含一些小的因數，而且已知有很多優良的演演算法可以去分解那些已知只有小因數的數字。例如小質數的試除法、SQUFOF、波拉德 ρ，以及ECM，以上是經常作為一起使用的方法。

基本上很多 $y(x)$ 值都會是可行的，因此因式分解過程的尾聲不需要是完全可信的；通常此過程大約有 5% 的輸入會出現異常，此時需要做少量的額外篩選。

## 基本篩選的例子

以下例子將演示沒有作對數優化或是質數次方的標準的二次篩法。令要分解的數為 $N=15347$，因此平方根 $N$ 無條件進位為124。由於 $N$ 很小，因此基本的多項式即足夠了： $y(x)=(x+124)^2-15347$。

### 數據收集

因為 $N$ 為小數字，所以只需 4 個質數。滿足在模 $p$ 下 15347 有一平方根的前 4 個質數 $p$ 為 2、17、23 以及 29(換句話說，對這些質數來說，15347是一個模這些數字的二次剩餘)。這些質數將是篩選的基礎。

現在我們要建造出我們的篩選 $V_X$從 $Y(X) = (X + \lceil\sqrt{N}\rceil)^2 - N = (X+124)^2-15347$並開始對基底裡每個質數進行篩選，以下選擇篩出 $0\leq X<100$ 的那些 $Y(X)$：

$$\begin{align}V &= \begin{bmatrix} Y(0) & Y(1) & Y(2) & Y(3) & Y(4) & Y(5) & \cdots & Y(99) \end{bmatrix} \\ <pre> & =\begin{bmatrix} 29 & 278 & 529 & 782 & 1037 & 1294 & \cdots & 34382 \end{bmatrix}\end{align}$$

下一步即是去作篩選的動作。對於我們的質數基底 $\lbrace 2, 17, 23, 29\rbrace$中的每一個質數 $p$ 值去解決以下方程式：

$$Y(X) \equiv (X + \left\lceil\sqrt{N}\right\rceil)^2 - N \equiv 0 \pmod{p}$$

找到陣列 $V$ 之中可被 $p$ 所整除的那些「元」。

對於 $p=2$ 解出 $(X + 124)^2 - 15347 \equiv 0 \pmod{2}$ 得到了 $X \equiv \sqrt{15347}-124 \equiv 1 \pmod{2}$。

所以，從 $X=1$ 開始每次 $+2$，每個「元」可被 2 整除。把那些元除以 2 之後得到：

$$V = \begin{bmatrix} 29 & 139 & 529 & 391 & 1037 & 647 & \cdots & 17191 \end{bmatrix}$$

同理，對於剩下的質數 $p$ $\lbrace 17, 23, 29\rbrace$ 方程式$X \equiv \sqrt{15347} - 124 \pmod{p}$ 也解決了。

值得注意的是，對於每一個 $p>2$，因為有兩個模平方根，因此得到 2 個線性方程式。

$$\begin{align} <pre> X & \equiv \sqrt{15347} - 124 & \equiv 8 - 124 & \equiv 3\pmod{17} \\ & & \equiv 9 - 124 & \equiv 4\pmod{17} \\ X & \equiv \sqrt{15347} - 124 & \equiv 11 - 124 & \equiv 2\pmod{23} \\ & & \equiv 12 - 124 & \equiv 3\pmod{23} \\ X & \equiv \sqrt{15347} - 124 & \equiv 8 - 124 & \equiv 0\pmod{29} \\ & & \equiv 21 - 124 & \equiv 13\pmod{29} \\ </pre> \end{align}$$

每個方程式 $X \equiv a \pmod{p}$ 導致 $V_x$ 從 $x=a$ 和之後每一次遞增一個 $p$ 值的那些項次皆可被 $p$ 整除。把 $V$ 中的 $a$、$a+p$、$a+2p$、$a+3p$等等的位置除以 $p$ , 如此對於每個在基底中的質數可以找到為相異質數的乘積（一次方）之光滑數。

$$V = \begin{bmatrix} 1 & 139 & 23 & 1 & 61 & 647 & \cdots & 17191 \end{bmatrix}$$

在 $V$ 之中的值等於一的那些「元」皆對應到一個光滑數。因為 $V_0$, $V_3$， $V_{71}$ 等於一，因此對應到：

$X+124$｜$Y$｜因數

$124$｜$29$｜$2^0\cdot17^0\cdot23^0\cdot29^1$

$127$｜$782$｜$2^1\cdot17^1\cdot23^1\cdot29^0$

$195$｜$22678$｜$2^1\cdot17^1\cdot23^1\cdot29^1$

### 矩陣處理

由於根據 $Y \equiv Z^2 \pmod{N}$的性質我們已經找到平滑數 $Y$，而演演算法接著的剩餘部分等同於狄克森因式分解法中的任何變體。

將方程式中的一個子集裡的指數乘積

$$\begin{align} 29 &= 2^0 \cdot 17^0 \cdot 23^0 \cdot 29^1 \\ 782 &= 2^1 \cdot 17^1 \cdot 23^1 \cdot 29^0 \\ 22678 &= 2^1 \cdot 17^1 \cdot 23^1 \cdot 29^1 \\ \end{align}$$

轉為一個矩陣形式 （在模 2 下）得到以下方程式：

$$S \cdot \begin{bmatrix} 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 0 \\ 1 & 1 & 1 & 1 \end{bmatrix} \equiv \begin{bmatrix} 0 & 0 & 0 & 0 \end{bmatrix} \pmod{2}$$

此方程式可由零空間（null space）的概念所給出一個解，為：

$$S = \begin{bmatrix}1 & 1 & 1 \end{bmatrix}$$

因此三個方程式的乘積產生了一個平方數（模 $N$ 之下）：

$$29 \cdot 782 \cdot 22678 = 22678^2$$

以及

$$124^2 \cdot 127^2 \cdot 195^2 = 3070860^2$$

所以此演演算法找到了

$$22678^2 \equiv 3070860^2 \pmod{15347}$$

測試其結果得到 $\gcd(3070860 - 22678, 15347)=103$，為 15347 的一個非平凡因數，而另一個為149。

而以上恰好顯示出，二次篩法只適用於 $n$ 值較大時。對於例如像 15347 這類的小數字，此演算法顯得過猶不及。試除法或是波拉德 ρ都可以在少量許多的計算之下找到一個因數。

## 倍數多項式

在實際用途上，有許多相異的多項式用在 $y$ 上，因為僅僅一個多項式通常不足以產生出對於因數基底的光滑數對 $(x,y)$。使用的多項式使用必須要有一個特別形式，因為它們需要為模 $n$. 下的平方數。多項式必定會與原始的 $y(x)=x^2-n$ 有類似的形式：

$$y(x)=(Ax+B)^2-n \qquad A,B\in\mathbb{Z}$$

假設 $B^2-n$ 是 $A$ 的一個倍數，則 $B^2-n = AC$ 且多項式 $y(x)$ 可以被寫作 $y(x) = A\cdot(Ax^2+2Bx+C)$。而如果 $A$ 為一個完全平方數，則只需考慮 $(Ax^2+2Bx+C)$的部分。

此方式（稱為 MPQS，倍數多項式二次篩選法（Multiple Polynomial Quadratic Sieve））非常適合平行運算，因為每一個處理因式分解的處理器可以單純的給入 $n$、因數基底以及多項式的集合，且直到運算完多項式之前都不須跟中央處理器作任何傳輸。

## 大質數

### 單一的大質數

如果在除以所有小於 $A$ 的因數之後，剩餘的數字（餘因子）小於 $A^2$，那麼這個餘因子必為質數。實際上，藉由對於餘因子去排序「關係」表，則它可以添加進因數基底裡。如果 $y(a)=7\times 11\times 23\times 137$ 且 $y(b)=3\times 5\times 7\times 137$, 則 $y(a)\times y(b)==3\times 5\times 11\times 23\times 7^2\times 137^2$。此可以降低以上完整執行因式分解的篩選陣列之「元」的臨界值。

### 更多的大質數

甚至更進一步去降低臨界值，並且使用一個高效處理將 $y(x)$ 之值分解為一些更大的質數之積（ECM 適合處理這樣子的東西）可以找到因數大多在因數基底，但有兩個甚至三個大質數的「關係」。循環的尋找過程因此允許一個共享好幾個質數的「關係」集合，合併成為單一的「關係」。

## 實際例子下的參數

為了展示在一個有包含多個多項式以及大質數優化下的實作方式去跑實際例子，會有的典型參數選取，將一個 267 位元的半質數輸入進 msieve 中，產出了以下的參數：

- 試除因數分解截止於：27 位元
- 篩選區間（對於每個多項式）：393216（12 個大小為 32768 的區塊）
- 光滑數之上界：1300967 (共 50294 個質數)
- 對於多項式 A 的係數之因數數量：10 (見上面倍數多項式條目)
- 大質數之上界：128795733 (26 位元) (見上面大質數條目)
- 光滑數的發現數：有 25952 為直接篩出，另外的 24462 為藉由合併那些有大質數的數字所得出
- 最終矩陣的大小：50294 × 50414，藉由過濾法減少到 35750 × 35862
- 非平凡的線性相依之發現數：15
- 總執行時間 (在 1.6 GHz UltraSparc III 上)：35 分 39 秒
- 最大記憶體使用量：8 MB

## 整數分解的紀錄

直到發現普通數體篩選法（number field sieve, 簡稱 NFS）之前，二次篩法（QS）曾是已知漸近最快的通用整數分解演演算法。現在，倫斯特拉橢圓曲線因式分解具有跟 QS 有相同的漸近運行時間（在 $n$ 由兩個相同大小級別的質數相乘所得的情況下），但在實際情況中，QS 速度更快，因為它採用的是單精度浮點數操作而不是橢圓曲線所使用的高精度計算。

在 1994 年的 4 月，RSA-129 的因數分解藉由 QS 完成了。其為一個由兩個大質數相乘的129 位數數字一個因數為 64 位長而另一個為 65 位。此因數分解的因數基底包含了 524339 個質數。數據收集階段花了 5000 個 MIPS 年，其完成於網際網路上的分散式計算。數據收集總量為 2 GB。數據處理花了45個小時在 Bellcore ( 現為 Telcordia 科技公司) 的 MasPar （大規模的平行化）超級電腦。這曾是最大的、藉由通用演演算法的公開分解，直至 NFS 被用於分解 RSA-130，於 1996 年 4 月 10 日 完成。所有自此以後分解的 RSA號碼 皆使用 NFS。

目前 QS 的紀錄是 $2^{803}-2^{402}+1$ 的一個 135 位數長之餘因子，其為 $2^{1606}+1$的一個 Aurifeuillian因數，在 2001 年分解為 66 位以及 69 位數長的質因數。

## 實作

- PPMPQS and PPSIQS
- mpqs
- SIMPQS 是由 William Hart 編寫的自我初始化（self-initializing）的倍數多項式二次篩選法的快速實作。其提供大質數的優化變體，並使用 Jason Papadopoulos' block Lanczos 程式碼於線性代數階段.。SIMPQS 可以使用 qsieve 指令在 SageMath 電腦代數套件上存取，或是原始來源裡下載。SIMPQS 被優化用於 Athlon 和 Opteron 機器上，但仍可在最常見的 32、 64 位元的結構上運行。而其完全是由 C 語言編寫而成的。
- 由 Dario Alpern 所提供的 factoring applet，其在特定狀況之下會使用二次篩選法。
- PARI/GP 電腦代數套件包含著自我初始化（self-initializing）的倍數多項式二次篩選法的一個實作並有著大質數的優化變體。其源自於 Thomas Papanikolaou 以及 Xavier Roblot 的一個編寫給 LiDIA 計畫的篩選法。自我初始化的方法是基於一個 Thomas Sosnowski 的一篇論文上的一個點子。
- 一個二次篩選法的變體開放於 MAGMA 電腦代數套件。其基於 1995 年 Arjen Lenstra 一個使用於他自己的「透過電子郵件分解整數」計畫的一次實作。
- msieve，一個支援單個或雙個大質數的倍數多項式二次篩選法之實作，由 Jason Papadopoulos 所編寫。原始碼以及 Windows 的二進位檔案皆是公開的。
- YAFU，由 Ben Buhrow 所編寫，與 msieve 相似但是對於現今大多的處理器來說更快。其使用 Jason Papadopoulos' block Lanczos 程式碼。原始碼以及 Windows、Linus 的二進位檔案皆是公開的。
- Ariel，一個用於教學用途的二次篩選法 Java 簡易實作。
- java-math-library 包含著也許是編寫於 Java 最快的二次篩選法（PSIQS 4.0 的後繼者）。
- Java QS，一個開源的 Java 計畫包含著 QS 的基本實作。由 Ilya Gazman 於 2016 年 2 月 4 日所釋出。
- C Quadratic Sieve – 使用C語言編寫的二次篩選法演算法的公共領域實現。於2022年發布，支持對最多330位的數進行批量分解，並以JSON或CSV格式輸出結果。

## 外部連結

- Reference paper "The Quadratic Sieve Factoring Algorithm" by Eric Landquist