Sclub交友聊天~加入聊天室當版主 # Matlab 多项式 GCD

 本帖最后由 hbghlyj 于 2021-11-20 08:04 编辑 Matlab没有内置的求多项式的GCD的函数.这里有人做了一个函数function g = poly_gcd(p,q)      if conv(p,q) == 0,  g = 1;   return, end;         g1 = p(min(find(p)):end)/p(min(find(p)));         g2 = q(min(find(q)):end)/q(min(find(q)));     for k = 1:length(g1)+length(g2)         t12 = length(g1)-length(g2); t21 = -t12;         g3 = [g1,zeros(1,t21)]-[g2,zeros(1,t12)];         g3 = [g3(min(find(abs(g3)>0.0001)):end)];      if conv(g3,g2) == 0, g = g2; return; end;      if t12 >= 0, g1 = g2; end;   g2 = g3/g3(1);     end; 复制代码求整系数多项式的精确GCD(来自这里)function [g,P] = intgcd(p,q) % %    *** Exact GCD of a pair of integer polynomials *** %        By F C Chang     10/28/08 % %    The polynomial GCD of two given polynomials can be found    %    exactly if the polynomial coefficients are all integers. % %    The process involves only integer arithematics operations. %    The floating point errors may be completely eliminated. % %    The expected result is obtained by applying recurrently %    "Elemination of leading coefficients of two equal-degree %    polynomials." It is modified from "Monic polynomials       %    subtraction," presented previously. % %    One of the main application is to factorize a integer %    polynomial with multiple roots. % %    References: %      "GCD of Polynomials," 26 Jul 2008 (Updated 05 Sep 2008). %      "Factoring a multiple-root polynomial." 18 Mar 2008. % % % ------------------------------------------------------------- % %%  Sample run on MATLAB: % %  >> %  Create a pair of test integer polymonial p(x) and its %  >> %  derivative q(x) to find the exact GCD of p(x) and q(x). %  >> %  Also printout "Polynomial remainder sequence" (PRS). %  >> %   % %  >> Format long g %  >> p = poly([-1 -1 -1 -1 -2 -2 -2 -3 -3 -4]), q = polyder(p), % %       p = %                1         20        175        882       2835 %             6072       8777       8458       5204       1848 %              288      %       q = %               10        180       1400       6174      17010 %            30360      35108      25374      10408       1848 % %  >> [g,P] = intgcd(p,q);  celldisp(P), g, % %       P{1} = %                1         20        175        882       2835 %             6072       8777       8458       5204       1848 %              288       %       P{2} = %               10        180       1400       6174      17010 %            30360      35108      25374      10408       1848 %       P{3} = %               10        175       1323       5670      15180 %            26331      29603      20816       8316       1440      %       P{4} = %                5         77        504       1830       4029 %             5505       4558       2092        408      %       P{5} = %                7        105        670       2374       5107 %             6829       5544       2500        480      %       P{6} = %                7         89        470       1334       2195 %             2093       1072        228      %       P{7} = %                2         25        130        364        592 %              559        284         60       %       P{8} = %                1         10         40         82         91 %               52         12       %       P{9} = %                1         10         40         82         91 %               52         12 % %       g = %                1         10         40         82         91 %               52         12       % %  >> % Complete PRS may even be obtained by hand calculation! % % % ---------------------------------------------------------------- % % % %function [g,P] = intgcd(p,q) % % *** Exact GCD of a pair of integer polynomials *** %       by F C Chang    10/28/08   if length(p) < 2,   g = 1;  P{1} = 1;    return, end;       n = length(p)-1;  nc = max(find(abs(p)>0))-1;       m = length(q)-1;  mc = max(find(abs(q)>0))-1;       p2 = p(1:nc+1);       p3 = q(1:mc+1); for k = 1:nc+nc+2;       p1 =  recast(p2);       p2 = [recast(p3),zeros(1,length(p2)-length(p3))];       p3 = p1*p2(1)-p2*p1(1);       p3 = p3(min(find(abs(p3)>0)):end);       P{k} = p1(1:max(find(abs(p1)>0)));   if isempty(find(abs(p3)>0)), P{k+1} = P{k}; break; end;    end;       gc = P{k};       g = [gc,zeros(1,min(n-nc,m-mc))];       celldisp(P); g, function [q,v] = recast(p) % %  *** Recast a vector by dividing it with its own GCD *** %       by F C Chang     09/28/08     if abs(p) == 0,  q = abs(p);   v = 1;   return;  end;        v = p(min(find(abs(p)>0)));   p = p*sign(v);    for j = 1:length(p),          v = gcd(v,p(j));   if v==1,   break,  end;    end;        q = p/v;复制代码

 返回列表 回复 发帖