有理函數系列展開的最佳算法

我需要用C++編寫函數來有效地找到給定有理函數（P（x）/ Q（x））的泰勒級數的係數。有理函數系列展開的最佳算法

函數參數將是多項式的冪（在分母和分母中相等），兩個具有多項式係數和展開項數的數組。

我的想法如下。考慮身份

P(x)/Q(x) = R(x) + ...

哪裏R(x)是項等於我需要找到係數的數數的多項式。然後，我可以用Q(x)乘雙方並獲得

P(x) = R(x) * Q(x) 

R(x) * Q(x) - P(x) = 0

因此，所有的係數應該是零。這是具有O（n^3）算法來解決的方程組。 O（n^3）並不像我想要的那麼快。

有沒有更快的算法？

我知道系列的係數是滿足線性遞推關係的。這讓我覺得O（n）算法是可能的。

來源

2014-04-15 Somnium

@ DavidEisenstat謝謝，這簡化了任務。但是我怎樣才能找到長分裂的'1/Q（x）'？我認爲與分工我只找到商和餘數。 – Somnium

@DavidEisenstat謝謝！避免分離也會加速事情，因爲乘法將是* O（n * m）*，其中n - 項的數量，m - 度。順便說一句，我不會理解爲什麼我的答案被標記爲過於寬泛，它需要一個特定的算法。 – Somnium

@DavidEisenstat你去了:)也許我們應該啓動一個關於meta的[算法]狀態的話題。對於是否將這些問題提交給CS –

我將要描述的算法在數學上由formal power series來證明。泰勒系列的每個功能都有正式的功率系列。反過來說並不真實，但如果我們用泰勒級數對函數進行算術運算，並用泰勒級數得到函數，那麼我們可以用正式冪級數做相同的算術並得到相同的答案。

正式冪級數的長除法算法就像您可能在學校學到的long division算法。我將在示例(1 + 2 x)/(1 - x - x^2)上演示它，其係數等於Lucas numbers。

分母必須有一個非零常數項。我們首先寫出分子，這是第一個殘差。

   -------- 
1 - x - x^2) 1 + 2 x

[ 我們除以分母的常數項（1）殘留的最低階項（1），並把該商數往上頂。

   1 
      -------- 
1 - x - x^2) 1 + 2 x

現在我們通過1乘1 - x - x^2，並從當前剩餘減去它。

   1 
      -------- 
1 - x - x^2) 1 + 2 x 
       1 - x - x^2 
       ------------- 
        3 x + x^2

再次做。

   1 + 3 x 
      -------- 
1 - x - x^2) 1 + 2 x 
       1 - x - x^2 
       --------------- 
        3 x + x^2 
        3 x - 3 x^2 - 3 x^3 
        ------------------- 
         4 x^2 + 3 x^3

再次。

   1 + 3 x + 4 x^2 
      ---------------- 
1 - x - x^2) 1 + 2 x 
       1 - x - x^2 
       --------------- 
        3 x + x^2 
        3 x - 3 x^2 - 3 x^3 
        ------------------- 
         4 x^2 + 3 x^3 
         4 x^2 - 4 x^3 - 4 x^4 
         --------------------- 
           7 x^3 + 4 x^4

再次。

   1 + 3 x + 4 x^2 + 7 x^3 
      ------------------------ 
1 - x - x^2) 1 + 2 x 
       1 - x - x^2 
       --------------- 
        3 x + x^2 
        3 x - 3 x^2 - 3 x^3 
        ------------------- 
         4 x^2 + 3 x^3 
         4 x^2 - 4 x^3 - 4 x^4 
         --------------------- 
           7 x^3 + 4 x^4 
           7 x^3 - 7 x^4 - 7 x^4 
           --------------------- 
             11 x^4 + 7 x^5

各個部門都有點無聊，因爲我使用的除數與一家領先的1，但如果我用了，說，2 - 2 x - 2 x^2，那麼所有的商係數將由2分。

來源

2014-04-15 20:24:46

謝謝，這有助於很多！我想找到的所有擴展都是整數係數，這可能保證領先係數爲1. – Somnium

如果你仔細看看你的計劃中得到的系統，你會發現它已經是對角線的，並且不需要O（n^3）來解決。它簡單地退化成一個線性遞歸（P []，Q []和R []是相應的多項式的係數）：

R[0] = P[0]/Q[0] 
R[n] = (P[n] - sum{0..n-1}(R[i] * Q[n-i]))/Q[0]

由於Q是一個多項式，總和具有不超過deg(Q)術語（從而花費不變的時間來計算），使整體複雜度漸近線性。您也可以查看遞歸的矩陣表示（可能）更好的漸近。

來源

2014-04-15 22:21:19 user58697

它對於一個R（i）值是線性的，但對於所有值R（0），...，R（n）它變成* O（n^2）* – Somnium

不會。如果Q具有無限數量的係數，它將變成O（n^2）。總和沒有比deg（Q）更多的項，所以需要一段時間來計算。我會編輯我的答案 - 目前尚不清楚。謝謝。 – user58697

然後它變成O（n * q），其中Q - 的度數。看起來和接受的帖子一樣。 – Somnium

這可以在O(n log n)的時間內完成對任意P和n程度Q。更確切地說，這可以在M(n)中完成，其中M(n)是多項式乘法的複雜度，其本身可以在O(n log n)中完成。

首先，第一個n系列擴展的術語可以簡單地看作是n-1的多項式。

假設您有興趣在P(x)/Q(x)系列擴展的第n條款。存在如上定義的計算M(n)時間內的Q的倒數的算法。

逆T(x)的Q(x)滿足T(x) * Q(x) = 1 + O(x^N)。即T(x) * Q(x)恰恰是1加上一些誤差項，其係數都在我們感興趣的第一項術語後面，所以我們可以放棄它們。

現在P(x)/Q(x)只是P(x) * T(x)，這只是另一個多項式乘法。

您可以在我的開源庫Altruct中找到一個計算上述反轉的實現。請參閱series.h文件。假設你已經有了一個計算兩個多項式乘積的方法，那麼計算倒數的代碼約爲10行（分而治之）。

實際算法如下：假設Q(x) = 1 + a1*x + a2*x^2 + ...。如果a0不是1，則可以簡單地將Q(x)及其後的逆T(x)與a0分開。 Asume，在每一步你有L條款的逆使Q(x) * T_L(x) = 1 + x^L * E_L(x)一些錯誤E_L(x)。最初爲T_1(X) = 1。如果你在上面插入這個，你將得到Q(x) * T_1(x) = Q(x) = 1 + x^1 * E_1(x)對於一些E_1(x)這意味着這適用於L=1。現在讓我們在每一步加倍L。您可以從上一步獲得E_L(x)作爲E_L(x) = (Q(x) * T_L(x) - 1)/x^L，或執行方式，只需刪除產品的第一個L係數。然後，您可以將上一步中的T_2L(x)計算爲T_2L(x) = T_L(x) - x^L * E_L(x) * T_L(x)。錯誤將是E_2L(x) = - E_L(x)^2。現在讓我們來檢查誘導步驟是否成立。

Q(x) * T_2L(x) 
= Q(x) * (T_L(x) - x^L * E_L(x) * T_L(x)) 
= Q(x) * T_L(x) * (1 - x^L * E_L(x)) 
= (1 + x^L * E_L(x)) * (1 - x^L * E_L(x)) 
= 1^2 - (x^L * E_L(x))^2 
= 1 + x^2L * E_2L(x)

Q.E.D.

我敢肯定這是不可能計算多項式除法比乘法更高效，你可以在下表中看到，該算法比一次乘法慢了3次：

n  mul  inv  factor 
10^4  24 ms  80 ms 3,33x 
10^5  318 ms  950 ms 2,99x 
10^6 4.162 ms 12.258 ms 2,95x 
10^7 101.119 ms 294.894 ms 2,92x

來源

2016-10-24 23:43:52 plamenko

有理函數系列展開的最佳算法

回答

相關問題