Matlab 多项式 GCD

hbghlyj

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

复制代码