cubic_my.m

来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 33 行

function [nreal,root] = cubic_my(p,q,r) 
      
pi=4.*atan(1.);
a=(3.*q-p*p)/3.;
b=(2.*p^3-9.*p*q+27.*r)/27.;
d=b*b/4.+a^3/27.;
if d > 0.
	nreal=1;
else
	nreal=3;
end

if nreal ~= 3
	arg1=-b/2.+sqrt(d);
	sarg1=arg1/abs(arg1);
	arg2=-b/2.-sqrt(d);
	sarg2=arg2/abs(arg2);
	aa=sarg1*((abs(arg1))^(1./3.));
	bb=sarg2*((abs(arg2))^(1./3.));
	root(1)=aa+bb-p/3.;
else
	c=2.*sqrt(-a/3.);
	z=(3.*b)/(a*c);
	alpha=atan((sqrt(1.-z*z))/abs(z));
	if z < 0
		alpha=pi-alpha;
	end
	y=(1./3.)*alpha;
	root(1)=c*cos(y)-p/3.;
	root(2)=c*cos(y+2.*pi/3.)-p/3.;
	root(3)=c*cos(y+4.*pi/3.)-p/3.;
end