getf SCI/util/testexamples.sci reinit_for_test() %U=mopen('SCI/tests/automatic_tests/sfact_data.ref','rb'); //Simple polynomial example z = poly(0, 'z'); p = (z - 1/2) * (2 - z); if load_ref('p') then pause,end, w = sfact(p); %ans = w * numer(horner(w, 1/z)); if load_ref('%ans') then pause,end, //matrix example F1 = [z - 1/2,z + 1/2,z^2 + 2;1,z,-z;z^3 + 2 * z,z,1/2 - z]; P = F1 * gtild(F1, 'd');//P is symmetric F = sfact(P); if load_ref('F') then pause,end, %ans = roots(det(P)); if load_ref('%ans') then pause,end, %ans = roots(det(gtild(F, 'd'))); if load_ref('%ans') then pause,end, //The stable roots %ans = roots(det(F)); if load_ref('%ans') then pause,end, //The antistable roots %ans = clean(P - F * gtild(F, 'd')); if load_ref('%ans') then pause,end, //Example of continuous time use s = poly(0, 's'); p = -3 * (s + 1 + %i) * (s + 1 - %i) * (s + 0.5) * (s - 0.5) * (s - (1 + %i)) * (s - (1 - %i));p = real(p); //p(s) = polynomial in s^2 , looks for stable f such that p=f(s)*f(-s) w = horner(p, (1 - s)/(1 + s));// bilinear transform w=p((1-s)/(1+s)) wn = numer(w);//take the numerator fn = sfact(wn);f = numer(horner(fn, (1 - s)/(s + 1)));//Factor and back transform f = f/sqrt(horner(f * gtild(f, 'c'), 0));f = f * sqrt(horner(p, 0));//normalization %ans = roots(f); if load_ref('%ans') then pause,end, //f is stable %ans = clean(f * gtild(f, 'c') - p); if load_ref('%ans') then pause,end, //f(s)*f(-s) is p(s) xdel_run(winsid()); mclose(%U);