繁體
|
簡體
Sclub交友聊天~加入聊天室當版主
(檢舉)
分享
新浪微博
QQ空间
人人网
腾讯微博
Facebook
Google+
Plurk
Twitter
Line
标题:
Matlab 多项式 GCD
[打印本页]
作者:
hbghlyj
时间:
2021-11-20 07:57
标题:
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;
复制代码
欢迎光临 悠闲数学娱乐论坛(第2版) (http://kuing.orzweb.net/)
Powered by Discuz! 7.2