免費論壇 繁體 | 簡體
Sclub交友聊天~加入聊天室當版主
分享
返回列表 发帖

Matlab 多项式 GCD

本帖最后由 hbghlyj 于 2021-11-20 08:04 编辑

Matlab没有内置的求多项式的GCD的函数.这里有人做了一个函数
  1. function g = poly_gcd(p,q)
  2.      if conv(p,q) == 0,  g = 1;   return, end;
  3.         g1 = p(min(find(p)):end)/p(min(find(p)));
  4.         g2 = q(min(find(q)):end)/q(min(find(q)));
  5.     for k = 1:length(g1)+length(g2)
  6.         t12 = length(g1)-length(g2); t21 = -t12;
  7.         g3 = [g1,zeros(1,t21)]-[g2,zeros(1,t12)];
  8.         g3 = [g3(min(find(abs(g3)>0.0001)):end)];
  9.      if conv(g3,g2) == 0, g = g2; return; end;
  10.      if t12 >= 0, g1 = g2; end;   g2 = g3/g3(1);
  11.     end;
复制代码
求整系数多项式的精确GCD(来自这里)
  1. function [g,P] = intgcd(p,q)
  2. %
  3. %    *** Exact GCD of a pair of integer polynomials ***
  4. %        By F C Chang     10/28/08
  5. %
  6. %    The polynomial GCD of two given polynomials can be found   
  7. %    exactly if the polynomial coefficients are all integers.
  8. %
  9. %    The process involves only integer arithematics operations.
  10. %    The floating point errors may be completely eliminated.
  11. %
  12. %    The expected result is obtained by applying recurrently
  13. %    "Elemination of leading coefficients of two equal-degree
  14. %    polynomials." It is modified from "Monic polynomials      
  15. %    subtraction," presented previously.
  16. %
  17. %    One of the main application is to factorize a integer
  18. %    polynomial with multiple roots.
  19. %
  20. %    References:
  21. %      "GCD of Polynomials," 26 Jul 2008 (Updated 05 Sep 2008).
  22. %      "Factoring a multiple-root polynomial." 18 Mar 2008.
  23. %
  24. %
  25. % -------------------------------------------------------------
  26. %
  27. %%  Sample run on MATLAB:
  28. %
  29. %  >> %  Create a pair of test integer polymonial p(x) and its
  30. %  >> %  derivative q(x) to find the exact GCD of p(x) and q(x).
  31. %  >> %  Also printout "Polynomial remainder sequence" (PRS).
  32. %  >> %  
  33. %
  34. %  >> Format long g
  35. %  >> p = poly([-1 -1 -1 -1 -2 -2 -2 -3 -3 -4]), q = polyder(p),
  36. %
  37. %       p =
  38. %                1         20        175        882       2835
  39. %             6072       8777       8458       5204       1848
  40. %              288     
  41. %       q =
  42. %               10        180       1400       6174      17010
  43. %            30360      35108      25374      10408       1848
  44. %
  45. %  >> [g,P] = intgcd(p,q);  celldisp(P), g,
  46. %
  47. %       P{1} =
  48. %                1         20        175        882       2835
  49. %             6072       8777       8458       5204       1848
  50. %              288      
  51. %       P{2} =
  52. %               10        180       1400       6174      17010
  53. %            30360      35108      25374      10408       1848
  54. %       P{3} =
  55. %               10        175       1323       5670      15180
  56. %            26331      29603      20816       8316       1440     
  57. %       P{4} =
  58. %                5         77        504       1830       4029
  59. %             5505       4558       2092        408     
  60. %       P{5} =
  61. %                7        105        670       2374       5107
  62. %             6829       5544       2500        480     
  63. %       P{6} =
  64. %                7         89        470       1334       2195
  65. %             2093       1072        228     
  66. %       P{7} =
  67. %                2         25        130        364        592
  68. %              559        284         60      
  69. %       P{8} =
  70. %                1         10         40         82         91
  71. %               52         12      
  72. %       P{9} =
  73. %                1         10         40         82         91
  74. %               52         12
  75. %
  76. %       g =
  77. %                1         10         40         82         91
  78. %               52         12      
  79. %
  80. %  >> % Complete PRS may even be obtained by hand calculation! %
  81. %
  82. % ----------------------------------------------------------------
  83. %
  84. %
  85. %
  86. %function [g,P] = intgcd(p,q)
  87. %
  88. % *** Exact GCD of a pair of integer polynomials ***
  89. %       by F C Chang    10/28/08

  90.   if length(p) < 2,   g = 1;  P{1} = 1;    return, end;
  91.       n = length(p)-1;  nc = max(find(abs(p)>0))-1;
  92.       m = length(q)-1;  mc = max(find(abs(q)>0))-1;
  93.       p2 = p(1:nc+1);
  94.       p3 = q(1:mc+1);
  95. for k = 1:nc+nc+2;
  96.       p1 =  recast(p2);
  97.       p2 = [recast(p3),zeros(1,length(p2)-length(p3))];
  98.       p3 = p1*p2(1)-p2*p1(1);
  99.       p3 = p3(min(find(abs(p3)>0)):end);
  100.       P{k} = p1(1:max(find(abs(p1)>0)));
  101.   if isempty(find(abs(p3)>0)), P{k+1} = P{k}; break; end;   
  102. end;
  103.       gc = P{k};
  104.       g = [gc,zeros(1,min(n-nc,m-mc))];
  105.       celldisp(P); g,
  106. function [q,v] = recast(p)
  107. %
  108. %  *** Recast a vector by dividing it with its own GCD ***
  109. %       by F C Chang     09/28/08
  110.     if abs(p) == 0,  q = abs(p);   v = 1;   return;  end;
  111.        v = p(min(find(abs(p)>0)));   p = p*sign(v);
  112.    for j = 1:length(p),  
  113.        v = gcd(v,p(j));   if v==1,   break,  end;
  114.    end;
  115.        q = p/v;
复制代码

返回列表 回复 发帖