1//<-- CLI SHELL MODE -->
2// =============================================================================
3// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4// Copyright (C) ????-2008 - INRIA
5//
6//  This file is distributed under the same license as the Scilab package.
7// =============================================================================
8//
9
10deff('[ok]=cmpr(h1,h2,eps)',['h1=h1-h2;';
11         'if norm(coeff(h1(2)))>eps then ok=0,else ok=1,end'])
12
13s = poly(0,'s');
14//
15//rationals
16//
17num = 1;
18den = 1 + s;
19eps = 5000 * %eps;
20assert_checkequal(cmpr(num / den, rlist(num, den, []), eps), 1);
21assert_checkequal(cmpr(den \ num, rlist(num, den, []), eps), 1);
22assert_checkequal(cmpr(num./den, rlist(num, den, []), eps), 1);
23assert_checkequal(cmpr(den.\num, rlist(num, den, []), eps), 1);
24
25num = 1.5 + s ** 3;
26assert_checkequal(cmpr(num / den, rlist(num, den, []), eps), 1);
27assert_checkequal(cmpr(den \ num, rlist(num, den, []), eps), 1);
28assert_checkequal(cmpr(num./den, rlist(num, den, []), eps), 1);
29assert_checkequal(cmpr(den.\num, rlist(num, den, []), eps), 1);
30
31h1 = num / den;
32x = 1.5;
33assert_checkequal(cmpr(h1 + x, rlist(num + x * den, den, []), eps), 1);
34assert_checkequal(cmpr(x + h1, rlist(num + x * den, den, []), eps), 1);
35assert_checkequal(cmpr(h1 - x, rlist(num - x * den, den, []), eps), 1);
36assert_checkequal(cmpr(x - h1, rlist(-num + x * den, den, []), eps), 1);
37
38x = 1.5 + 3 * s;
39assert_checkequal(cmpr(h1 + x, rlist(num + x * den, den, []), eps), 1);
40assert_checkequal(cmpr(x + h1, rlist(num + x * den, den, []), eps), 1);
41assert_checkequal(cmpr(h1 - x, rlist(num - x * den, den, []), eps), 1);
42assert_checkequal(cmpr(x - h1, rlist(-num + x * den, den, []), eps), 1);
43
44y= s ** 3;
45h2 = x / y;
46assert_checkequal(cmpr(h1 + h2, rlist(num * y + x * den, den * y, []), eps), 1);
47assert_checkequal(cmpr(h1 - h2, rlist(num * y - x * den, den * y, []), eps), 1);
48
49//
50// concatenations
51//
52h1 = num / den;
53x = 1.5;
54assert_checkequal(cmpr([h1, x], rlist([num, x], [den, 1], []), eps), 1);
55assert_checkequal(cmpr([x, h1], rlist([x, num], [1, den], []), eps), 1);
56assert_checkequal(cmpr([h1; x], rlist([num; x], [den; 1], []), eps), 1);
57assert_checkequal(cmpr([x; h1], rlist([x; num], [1; den], []), eps), 1);
58
59x = 1.5 + 3 * s;
60assert_checkequal(cmpr([h1, x], rlist([num, x], [den, 1], []), eps), 1);
61assert_checkequal(cmpr([x, h1], rlist([x, num], [1, den], []), eps), 1);
62assert_checkequal(cmpr([h1; x], rlist([num; x], [den; 1], []), eps), 1);
63assert_checkequal(cmpr([x; h1], rlist([x; num], [1; den], []), eps), 1);
64
65y = -0.5 + s ** 3;
66h2 = x / y;
67assert_checkequal(cmpr([h1, h2], rlist([num, x], [den, y], []), eps), 1);
68assert_checkequal(cmpr([h1; h2], rlist([num; x], [den; y], []), eps), 1);
69
70h1 = [num / den, den / num];
71x = [0.3 1.5];
72assert_checkequal(cmpr([h1, x], rlist([h1(2), x], [h1(3), ones(x)], []), eps), 1);
73assert_checkequal(cmpr([x, h1], rlist([x, h1(2)], [ones(x), h1(3)], []), eps), 1);
74assert_checkequal(cmpr([h1; x], rlist([h1(2); x], [h1(3); ones(x)], []), eps), 1);
75assert_checkequal(cmpr([x; h1], rlist([x; h1(2)], [ones(x); h1(3)], []), eps), 1);
76
77h1 = [num / den; den / num];
78x = [0.3; -1.5];
79assert_checkequal(cmpr([h1, x], rlist([h1(2), x], [h1(3), ones(x)], []), eps), 1);
80assert_checkequal(cmpr([x, h1], rlist([x, h1(2)], [ones(x), h1(3)], []), eps), 1);
81assert_checkequal(cmpr([h1; x], rlist([h1(2); x], [h1(3); ones(x)], []), eps), 1);
82assert_checkequal(cmpr([x; h1], rlist([x; h1(2)], [ones(x); h1(3)], []), eps), 1);
83
84x = [1.5 + 3 * s; -1 + s ** 3];
85assert_checkequal(cmpr([h1, x], rlist([h1(2), x], [h1(3), ones(x)], []), eps), 1);
86assert_checkequal(cmpr([x, h1], rlist([x, h1(2)], [ones(x), h1(3)], []), eps), 1);
87assert_checkequal(cmpr([h1; x], rlist([h1(2); x], [h1(3); ones(x)], []), eps), 1);
88assert_checkequal(cmpr([x; h1], rlist([x; h1(2)], [ones(x); h1(3)], []), eps), 1);
89
90h1 = [num / den, den / num];
91x = [1.5 + 3 * s, -1 + s ** 2];
92assert_checkequal(cmpr([h1, x], rlist([h1(2), x], [h1(3), ones(x)], []), eps), 1);
93assert_checkequal(cmpr([x, h1], rlist([x, h1(2)], [ones(x), h1(3)], []), eps), 1);
94assert_checkequal(cmpr([h1; x], rlist([h1(2); x], [h1(3); ones(x)], []), eps), 1);
95assert_checkequal(cmpr([x; h1], rlist([x; h1(2)], [ones(x); h1(3)], []), eps), 1);
96
97h1 = [num / den; den / num];
98y = -0.5 + s ** 3;
99h2 = [num / y; y * y / (y + 1)];
100assert_checkequal(cmpr([h1, h2], rlist([h1(2), h2(2)], [h1(3), h2(3)], []), eps), 1);
101assert_checkequal(cmpr([h1; h2], rlist([h1(2); h2(2)], [h1(3); h2(3)], []), eps), 1);
102
103h1 = [num / den, den / num];
104y = -0.5 + s ** 3;
105h2 = [num / y, y * y / (y + 1)];
106assert_checkequal(cmpr([h1, h2], rlist([h1(2), h2(2)], [h1(3), h2(3)], []), eps), 1);
107assert_checkequal(cmpr([h1; h2], rlist([h1(2); h2(2)], [h1(3); h2(3)], []), eps), 1);
108
109//
110// extraction
111//
112h1 = [num / den, den / num];
113assert_checkequal(cmpr(h1(1, 1), num / den, eps), 1);
114assert_checkequal(cmpr(h1(1, [2 1]), [den / num, num / den], eps), 1);
115
116h1 = [num / den; den / num];
117assert_checkequal(cmpr(h1(2, 1), den / num, eps), 1);
118assert_checkequal(cmpr(h1(1:2, 1), h1, eps), 1);
119assert_checkequal(cmpr(h1([2 1], 1), [den / num ; num / den], eps), 1);
120
121y = -0.5 + s ** 3;
122h1 = [num / den, den / num];
123h2 = [num / den, den / num; num / y, y * y / (y + 1)];
124assert_checkequal(cmpr(h2(2, 1), num / y, eps), 1);
125assert_checkequal(cmpr(h2(1:2, 1:2), h2, eps), 1);
126assert_checkequal(cmpr(h2([2 1], 1), [num / y; num / den], eps), 1);
127assert_checkequal(cmpr(h2(:, 1), [num / den ; num / y], eps), 1);
128assert_checkequal(cmpr(h2(2, :), [num / y, y * y / (y + 1)], eps), 1);
129assert_checkequal(cmpr(h2(:, :), h2, eps), 1);
130
131//
132// insertions
133//
134h1 = [num / den, den / num];
135x = 0.33;
136hh = h1;
137hh(1, 1) = x;
138assert_checkequal(cmpr(hh, [x, h1(1, 2)], eps), 1);
139
140x = [-2.67 0.8];
141hh = h1;
142hh(1, 1:2) = x;
143assert_checkequal(cmpr(hh, x, eps), 1);
144
145hh = h1;
146hh(1, [2 1]) = x;
147assert_checkequal(cmpr(hh, x([2, 1]), eps), 1);
148
149hh = h1;
150hh(1, [2; 1]) = x;
151assert_checkequal(cmpr(hh, x([2, 1]), eps), 1);
152
153hh = h1;
154hh(1, :) = x;
155assert_checkequal(cmpr(hh, x, eps), 1);
156
157hh = h1;
158hh(:,:) = x;
159assert_checkequal(cmpr(hh, x, eps), 1);
160
161h1 = [num / den; den / num];
162x = 0.33;
163hh = h1;
164hh(1,1) = x;
165assert_checkequal(cmpr(hh, [x; h1(2,1)], eps), 1);
166
167x = [-2.67; 0.8];
168hh = h1;
169hh(1:2,1) = x;
170assert_checkequal(cmpr(hh, x, eps), 1);
171
172hh = h1;
173hh(:,1) = x;
174assert_checkequal(cmpr(hh, x, eps), 1);
175
176hh = h1;
177hh(:,:) = x;
178assert_checkequal(cmpr(hh, x, eps), 1);
179
180hh = h1;
181hh([2 1],1) = x;
182assert_checkequal(cmpr(hh, x([2 1]), eps), 1);
183
184hh = h1;
185hh([2;1], 1) = x;
186assert_checkequal(cmpr(hh, x([2 1]), eps), 1);
187
188h1 = [num / den, den / num];
189x = 0.33 * s + 1;
190hh = h1;
191hh(1,1) = x;
192assert_checkequal(cmpr(hh, [x, h1(1, 2)], eps), 1);
193
194x = [-2.67 0.8 + s ** 3];
195hh = h1;
196hh(1,1:2) = x;
197assert_checkequal(cmpr(hh, x, eps), 1);
198
199hh = h1;
200hh(1,[2 1]) = x;
201assert_checkequal(cmpr(hh, x([2 1]), eps), 1);
202
203hh = h1;
204hh(1,[2;1]) = x;
205assert_checkequal(cmpr(hh, x([2 1]), eps), 1);
206
207hh = h1;
208hh(1,:) = x;
209assert_checkequal(cmpr(hh, x, eps), 1);
210
211hh = h1;
212hh(:,:) = x;
213assert_checkequal(cmpr(hh, x, eps), 1);
214
215h1 = [num / den; den / num];
216x = -0.33 + 38 * s;
217hh = h1;
218hh(1,1) = x;
219assert_checkequal(cmpr(hh, [x; h1(2,1)], eps), 1);
220
221x = [-2.67 - s * 8; 0.8];
222hh = h1;
223hh(1:2,1) = x;
224assert_checkequal(cmpr(hh, x, eps), 1);
225
226hh = h1;
227hh(:,1) = x;
228assert_checkequal(cmpr(hh, x, eps), 1);
229
230hh = h1;
231hh(:,:) = x;
232assert_checkequal(cmpr(hh, x, eps), 1);
233
234hh = h1;
235hh([2 1], 1) = x;
236assert_checkequal(cmpr(hh, x([2 1]), eps), 1);
237
238hh = h1;
239hh([2;1],1) = x;
240assert_checkequal(cmpr(hh, x([2 1]), eps), 1);
241
242h1 = [num / den, den / num];
243y = 0.33 * s + 1;
244x = y * y / (y + 1);
245hh = h1;
246hh(1,1) = x;
247assert_checkequal(cmpr(hh, [x, h1(1,2)], eps), 1);
248
249x = [num / y, y * y / (y + 1)];
250hh = h1;
251hh(1,1:2) = x;
252assert_checkequal(cmpr(hh, x, eps), 1);
253
254hh = h1;
255hh(1,[2 1]) = x;
256assert_checkequal(cmpr(hh, x(1,[2,1]), eps), 1);
257
258hh = h1;
259hh(1,[2;1]) = x;
260assert_checkequal(cmpr(hh, x(1,[2,1]), eps), 1);
261
262hh = h1;
263hh(1,:) = x;
264assert_checkequal(cmpr(hh, x, eps), 1);
265
266hh = h1;
267hh(:,:) = x;
268assert_checkequal(cmpr(hh, x, eps), 1);
269
270h1 = [num / den; den / num];
271x = [num / y; y *y / (y + 1)];
272hh = h1;
273hh(1:2,1) = x;
274assert_checkequal(cmpr(hh, x, eps), 1);
275
276hh = h1;
277hh(:,1) = x;
278assert_checkequal(cmpr(hh, x, eps), 1);
279
280hh = h1;
281hh(:,:) = x;
282assert_checkequal(cmpr(hh, x, eps), 1);
283
284hh = h1;
285hh([2 1],1) = x;
286assert_checkequal(cmpr(hh, x([2 1], 1), eps), 1);
287
288hh = h1;
289hh([2;1],1) = x;
290assert_checkequal(cmpr(hh, x([2 1], 1), eps), 1);
291
292//
293// matrix operations
294//
295h1 = [num / den, den / num];
296x = [0.3 1.5];
297assert_checkequal(cmpr(-h1, (-1) * h1, eps), 1);
298assert_checkequal(cmpr(h1 + x, [h1(1,1) + x(1,1) h1(1,2) + x(1,2)], eps), 1);
299assert_checkequal(cmpr(x + h1, [h1(1,1) + x(1,1) h1(1,2) + x(1,2)], eps), 1);
300assert_checkequal(cmpr(h1 - x, [h1(1,1) - x(1,1) h1(1,2) - x(1,2)], eps), 1);
301assert_checkequal(cmpr(x - h1, [-h1(1,1) + x(1,1) -h1(1,2) + x(1,2)], eps), 1);
302
303h1 = [num / den; den / num];
304x = [0.3; 1.5];
305assert_checkequal(cmpr(-h1, (-1) * h1, eps), 1);
306assert_checkequal(cmpr(h1 + x, [h1(1,1) + x(1,1); h1(2,1) + x(2,1)], eps), 1);
307assert_checkequal(cmpr(x + h1, [h1(1,1) + x(1,1); h1(2,1) + x(2,1)], eps), 1);
308assert_checkequal(cmpr(h1 - x, [h1(1,1) - x(1,1); h1(2,1) - x(2,1)], eps), 1);
309assert_checkequal(cmpr(x - h1, [-h1(1,1) + x(1,1); -h1(2,1) + x(2,1)], eps), 1);
310
311h1 = [num / den, den / num];
312x=[0.3 + s, 1.5];
313assert_checkequal(cmpr(-h1, (-1) * h1, eps), 1);
314assert_checkequal(cmpr(h1 + x, [h1(1,1) + x(1,1) h1(1,2) + x(1,2)], eps), 1);
315assert_checkequal(cmpr(x + h1,[h1(1,1) + x(1,1) h1(1,2) + x(1,2)], eps), 1);
316assert_checkequal(cmpr(h1 - x,[h1(1,1) - x(1,1) h1(1,2) - x(1,2)], eps), 1);
317assert_checkequal(cmpr(x - h1,[-h1(1,1) + x(1,1) -h1(1,2) + x(1,2)], eps), 1);
318
319h1 = [num / den; den / num];
320x = [0.3; 1.5 - 3 * s];
321assert_checkequal(cmpr(-h1, (-1) * h1, eps), 1);
322assert_checkequal(cmpr(h1 + x, [h1(1,1) + x(1,1); h1(2,1) + x(2,1)], eps), 1);
323assert_checkequal(cmpr(x + h1, [h1(1,1) + x(1,1); h1(2,1) + x(2,1)], eps), 1);
324assert_checkequal(cmpr(h1 - x, [h1(1,1) - x(1,1); h1(2,1) - x(2,1)], eps), 1);
325assert_checkequal(cmpr(x - h1, [-h1(1,1) + x(1,1); -h1(2,1) + x(2,1)], eps), 1);
326
327//
328h1 = [num / den, den / num];
329assert_checkequal(cmpr([num, den] / den, [num / den, 1], eps), 1);
330assert_checkequal(cmpr(den \ [num, den], [num / den, 1], eps), 1);
331assert_checkequal(cmpr([num, den]./[den, num], h1, eps), 1);
332assert_checkequal(cmpr([den, num].\[num, den], h1, eps), 1);
333
334h1 = [num / den, den / num];
335x = [0.3 1.5];
336assert_checkequal(cmpr(h1 / x(1), [h1(1,1) / x(1), h1(1,2) / x(1)], eps), 1);
337assert_checkequal(cmpr(x(1) \ h1, [h1(1,1) / x(1), h1(1,2) / x(1)], eps), 1);
338assert_checkequal(cmpr(h1./ x, [h1(1,1) / x(1), h1(1,2) / x(2)], eps), 1);
339assert_checkequal(cmpr(x.\ h1, [h1(1,1) / x(1), h1(1,2) / x(2)], eps), 1);
340assert_checkequal(cmpr(h1 * x(1), [h1(1,1) * x(1), h1(1,2) * x(1)], eps), 1);
341assert_checkequal(cmpr(h1 .* x, [h1(1,1) * x(1), h1(1,2) * x(2)], eps), 1);
342
343h1 = [num / den; den / num];
344assert_checkequal(cmpr([num; den] / den, [num / den; 1], eps), 1);
345assert_checkequal(cmpr(den \ [num; den],[num / den; 1], eps), 1);
346assert_checkequal(cmpr([num; den]./[den; num], h1, eps), 1);
347assert_checkequal(cmpr([den; num].\[num; den], h1, eps), 1);
348
349x = [0.3; 1.5];
350assert_checkequal(cmpr(h1 / x(1), [h1(1,1) / x(1); h1(2,1) / x(1)], eps), 1);
351assert_checkequal(cmpr(x(1) \ h1, [h1(1,1) / x(1) ; h1(2,1) / x(1)], eps), 1);
352assert_checkequal(cmpr(h1 ./ x, [h1(1,1) / x(1); h1(2,1) / x(2)], eps), 1);
353assert_checkequal(cmpr(x .\ h1, [h1(1,1) / x(1); h1(2,1) / x(2)], eps), 1);
354assert_checkequal(cmpr(h1 * x(1), [h1(1,1) * x(1); h1(2,1) * x(1)], eps), 1);
355assert_checkequal(cmpr(h1 .* x, [h1(1,1) * x(1); h1(2,1) * x(2)], eps), 1);
356assert_checkequal(cmpr(h1' * x, h1(1,1) * x(1) + h1(2,1) * x(2), eps), 1);
357assert_checkequal(cmpr(h1 * x', [h1(1,1) * x(1), h1(1,1) * x(2); h1(2,1) * x(1), h1(2,1) * x(2)], eps), 1);
358
359h1 = [num / den, den / num];
360x = [0.3 + 3 * s, 1.5];
361assert_checkequal(cmpr(h1 / x(1), [h1(1,1) / x(1), h1(1,2) / x(1)], eps), 1);
362assert_checkequal(cmpr(x(1) \ h1, [h1(1,1) / x(1), h1(1,2) / x(1)], eps), 1);
363assert_checkequal(cmpr(h1 ./ x, [h1(1,1) / x(1), h1(1,2) / x(2)], eps), 1);
364assert_checkequal(cmpr(x .\ h1, [h1(1,1) / x(1), h1(1,2) / x(2)], eps), 1);
365assert_checkequal(cmpr(h1 * x(1), [h1(1,1) * x(1), h1(1,2) * x(1)], eps), 1);
366assert_checkequal(cmpr(h1 .* x, [h1(1,1) * x(1), h1(1,2) * x(2)], eps), 1);
367assert_checkequal(cmpr(h1 * x', h1(1,1) * x(1) + h1(1,2) * x(2), eps), 1);
368assert_checkequal(cmpr(h1' * x, [h1(1,1) * x(1), h1(1,1) * x(2); h1(1,2) * x(1), h1(1,2) * x(2)], eps), 1);
369
370h1=[num/den;den/num];x=[0.3+3*s; 1.5];
371assert_checkequal(cmpr(h1 / x(1), [h1(1,1) / x(1); h1(2,1) / x(1)], eps), 1);
372assert_checkequal(cmpr(x(1) \ h1, [h1(1,1) / x(1) ; h1(2,1) / x(1)], eps), 1);
373assert_checkequal(cmpr(h1 ./ x, [h1(1,1) / x(1); h1(2,1) / x(2)], eps), 1);
374assert_checkequal(cmpr(x .\ h1, [h1(1,1) / x(1); h1(2,1) / x(2)], eps), 1);
375assert_checkequal(cmpr(h1 * x(1), [h1(1,1) * x(1); h1(2,1) * x(1)], eps), 1);
376assert_checkequal(cmpr(h1 .* x, [h1(1,1) * x(1); h1(2,1) * x(2)], eps), 1);
377assert_checkequal(cmpr(h1' * x, h1(1,1) * x(1) + h1(2,1) * x(2), eps), 1);
378assert_checkequal(cmpr(h1 * x', [h1(1,1) * x(1), h1(1,1) * x(2); h1(2,1) * x(1), h1(2,1) * x(2)], eps), 1);
379
380h1 = [num / den, den / num];
381x = [0.3 / s, 1.5 - s**2 / (1 + s**2)];
382assert_checkequal(cmpr(h1 / x(1,1), [h1(1,1) / x(1,1), h1(1,2) / x(1,1)], eps), 1);
383assert_checkequal(cmpr(x(1,1) \ h1, [h1(1,1) / x(1,1), h1(1,2) / x(1,1)], eps), 1);
384assert_checkequal(cmpr(h1 ./ x, [h1(1,1) / x(1,1), h1(1,2) / x(1,2)], eps), 1);
385assert_checkequal(cmpr(x .\ h1, [h1(1,1) / x(1,1), h1(1,2) / x(1,2)], eps), 1);
386assert_checkequal(cmpr(h1 * x(1,1), [h1(1,1) * x(1,1), h1(1,2) * x(1,1)], eps), 1);
387assert_checkequal(cmpr(h1 .* x, [h1(1,1) * x(1,1), h1(1,2) * x(1,2)], eps), 1);
388
389h1 = [num /den; den / num];
390x = [0.3 / s; 1.5 - s**2 / (1 + s**2)];
391assert_checkequal(cmpr(h1 / x(1,1), [h1(1,1) / x(1,1); h1(2,1) /  x(1,1)], eps), 1);
392assert_checkequal(cmpr(x(1,1) \ h1, [h1(1,1) / x(1,1); h1(2,1) / x(1,1)], eps), 1);
393assert_checkequal(cmpr(h1 ./ x, [h1(1,1) / x(1,1); h1(2,1) / x(2,1)], eps), 1);
394assert_checkequal(cmpr(x .\ h1, [h1(1,1) / x(1,1); h1(2,1) / x(2,1)], eps), 1);
395assert_checkequal(cmpr(h1 * x(1,1), [h1(1,1) * x(1,1); h1(2,1) * x(1,1)], eps), 1);
396assert_checkequal(cmpr(h1 .* x, [h1(1,1) * x(1,1); h1(2,1) * x(2,1)], eps), 1);
397assert_checkequal(cmpr(h1' * x, h1(1,1) * x(1,1) + h1(2,1) * x(2,1), eps), 1);
398assert_checkequal(cmpr(h1 * x', [h1(1,1) * x(1,1), h1(1,1) * x(2,1); h1(2,1) * x(1,1), h1(2,1) * x(2,1)], eps), 1);
399
400h1 = h1 * x';
401x = [1 2; 3 4];
402assert_checkequal(cmpr(h1 / x, h1 * inv(x), eps), 1);
403assert_checkequal(cmpr(x \ h1, inv(x) * h1, eps), 1);
404
405x = [s, s * s + s - 1; 1, s + 1];
406xi = [1 + s, 1 - s - s * s; -1, s];
407if norm(coeff(invr(x)-xi))>eps then pause,end
408assert_checkequal(cmpr(h1 / x, h1 * xi, eps), 1);
409assert_checkequal(cmpr(x \ h1, xi * h1, eps), 1);
410assert_checkequal(cmpr(x ** (-1), xi, eps), 1);
411
412x = [1 / (1 + s), 1 / (s * (s +1 )) - 1; 1, 1/s];
413xi = [1/s, -1 / (s * (s + 1)) + 1; -1, 1 / (1 + s)];
414assert_checkequal(cmpr(invr(x), xi, eps), 1);
415assert_checkequal(cmpr(h1 / x, h1 * xi, eps), 1);
416assert_checkequal(cmpr(x \ h1, xi * h1, eps), 1);
417assert_checkequal(cmpr((1/(1 + s)) ** 3, 1 / ((1 + s) ** 3), eps), 1);
418
419x = [1 / (1 + s), 1; 0, 1 / s];
420x3 = [x(1,1) ** 3, x(1,1) ** 2 + x(2,2) * (x(2,2) + x(1,1)); 0, x(2,2) ** 3];
421assert_checkequal(cmpr(x ** 3, x3, eps), 1);
422
423//Bezout
424
425//mode(5)
426//test
427un = poly(1, 's', 'c');
428zer = 0 * un;
429s = poly(0, 's');
430[p, q] = bezout(un, s);
431assert_checkalmostequal(norm(coeff([un s] * q - [p 0])), 0, [], 10 * %eps);
432
433[p, q] = bezout(s, un);
434assert_checkalmostequal(norm(coeff([s un] * q - [p 0])), 0, [], 10 * %eps);
435
436[p, q] = bezout(un, un);
437assert_checkalmostequal(norm(coeff([un un] * q - [p 0])), 0, [], 10 * %eps);
438
439[p, q] = bezout(zer, s);
440assert_checkalmostequal(norm(coeff([zer s] * q - [p 0])), 0, [], 10 * %eps);
441
442[p, q] = bezout(s, zer);
443assert_checkalmostequal(norm(coeff([s zer] * q - [p 0])), 0, [], 10 * %eps);
444
445[p, q] = bezout(zer, zer);
446assert_checkalmostequal(norm(coeff([zer zer] * q - [p 0])), 0, [], 10 * %eps);
447
448[p, q] = bezout(zer, un);
449assert_checkalmostequal(norm(coeff([zer un] * q - [p 0])), 0, [], 10 * %eps);
450
451[p,q] = bezout(un, zer);
452assert_checkalmostequal(norm(coeff([un zer] * q - [p 0])), 0, [], 10 * %eps);
453
454//simple
455a = poly([1 2 3], 'z');
456b = poly([4 1], 'z');
457
458[p q] = bezout(a, b);
459dt = q(1,1) * q(2,2) - q(1,2) * q(2,1);
460dt0 = coeff(dt, 0);
461assert_checkalmostequal(norm(coeff(dt / dt0 - 1)), 0, [], 10 * %eps);
462
463qi = [q(2,2) -q(1,2); -q(2,1) q(1,1)] / dt0;
464assert_checkalmostequal(norm(coeff([p 0] * qi - [a b])), 0, [], 100 * %eps);
465
466//more difficult
467b1 = poly([4 1 + 1000 * %eps], 'z');
468del = 10 * norm(coeff(b1 - b));
469[p, q] = bezout(a, b1);
470dt = q(1,1) * q(2,2) - q(1,2) * q(2,1);
471dt0 = coeff(dt, 0);
472assert_checkalmostequal(norm(coeff(dt / dt0 - 1)), 0, [], del);
473
474qi = [q(2,2) -q(1,2); -q(2,1) q(1,1)] / dt0;
475assert_checkalmostequal(norm(coeff([p 0] * qi - [a b1])), 0, [], del);
476
477b1 = poly([4 1 + .001], 'z');
478del = 10 * norm(coeff(b1 - b));
479[p, q]=bezout(a, b1);
480dt = q(1,1) * q(2,2) - q(1,2) * q(2,1);
481dt0 = coeff(dt, 0);
482assert_checkalmostequal(norm(coeff(dt / dt0 - 1)), 0, [], 100000*%eps);
483
484qi = [q(2,2) -q(1,2); -q(2,1) q(1,1)] / dt0;
485assert_checkalmostequal(norm(coeff([p 0] * qi - [a b1])), 0, [], 100000*%eps);
486
487b1 = poly([4 1 + 10000 * %eps], 'z');
488del = 10 * norm(coeff(b1 - b));
489[p, q] = bezout(a, b1);
490dt = q(1,1) * q(2,2) - q(1,2) * q(2,1);
491dt0 = coeff(dt, 0);
492assert_checkalmostequal(norm(coeff(dt / dt0 - 1)), 0, [], del);
493
494qi = [q(2,2) -q(1,2); -q(2,1) q(1,1)] / dt0;
495assert_checkalmostequal(norm(coeff([p 0] * qi - [a b1])), 0, [], del);
496
497z = poly(0, 'z');
498num =   0.99999999999999922 + z *(-4.24619123578530730 + ...
499        z * (10.0503215227460350 + z * (-14.6836461849931740 + ...
500        z * (13.924822877119892 + z * (-5.63165998008533460 + ...
501        z * (-5.63165998008530710 + z * (13.9248228771198730 + ...
502        z * (-14.683646184993167 + z * (10.0503215227460330 + ...
503        z * (-4.24619123578530910 + z * (0.99999999999999989)))))))))));
504den =   -0.17323463717888873 + z * (1.91435457459735380 + ...
505        z * (-9.90126732768255560 + z * (31.6286096652930410 + ...
506        z * (-69.3385546880304280 + z * (109.586435800377690+ ...
507        z * (-127.516160100808290 + z * (109.388684898145950+ ...
508        z * (-67.92645394857864+z*(29.1602681026148110+ ...
509        z * (-7.8212498781094952+z*(0.99999999999999989)))))))))));
510
511[p, q] = bezout(num, den);
512del = 1.d-4;
513dt = q(1,1) * q(2,2) - q(1,2) * q(2,1);
514dt0 = coeff(dt,0);
515assert_checkalmostequal(norm(coeff(dt / dt0 - 1)), 0, [], del);
516
517qi = [q(2,2) -q(1,2); -q(2,1) q(1,1)] / dt0;
518// JPC
519del = 3 * del;
520assert_checkalmostequal(norm(coeff([p 0] * qi - [num den])), 0, [], del);
521assert_checkalmostequal(degree(p), 0, [], del);
522
523// une autre "solution" telle que l'erreur directe est petite mais l'erreur
524// inverse grande
525q1 = [];
526q1(1,1) = poly([0.011533588674 , 3.570224012502 , -28.78817740957 ,...
527            102.9479419903, -210.8258579715 , 266.2028963639 ,...
528           -207.427018299 , 92.74478640319, -18.5259652457],'z','c');
529q1(1,2) = poly([0.000270220664 , -0.002465565223 , 0.010039706635 ,...
530            -0.023913827007, 0.036387144032 , -0.036175176058 ,...
531            0.022954475171 , -0.008514798968,  0.001417382492],'z','c');
532q1(2,1) = poly([0.000639302018 , 20.502472606 , -26.66621972106 ,...
533            39.74001534015,  -5.945830824342 , -7.973226036298 ,...
534            39.84118622788 , -26.51337424414, 18.5259652457],'z','c');
535q1(2,2) = poly( [0.001562930385 , -0.003589174974 , 0.005076237479 ,...
536            -0.003483568682,  -0.000135940266 , 0.003651509121 ,...
537            -0.005059502869 , 0.003447573440, -0.001417382492],'z','c');
538p1 = poly( [0.011422839421 , -0.029264363070 , 0.030070175223 ,...
539         -0.012596066108],'z','c');
540
541//simplification
542num = [0.03398330733500143, -0.20716057008572933, 0.64660689206696986,...
543     -1.97665462134021790, 3.38751027286891300, -3.58940006392108120,...
544      5.09956159043231680, 5.2514918861834694, 1.00000000000000020];
545den = [0.03398330733500360, -0.20716057008816283, 0.64660689204312,...
546      -1.97665462079896660, 3.38751027286891300, -3.58940006392108350,...
547       5.09956159043232040, 5.2514918861834712, 1];
548num = poly(num, 'z', 'c');
549den = poly(den, 'z', 'c');
550[p,q] = bezout(num, den);
551del = 1.d-8;
552dt = q(1,1) * q(2,2) - q(1,2) * q(2,1);
553dt0 = coeff(dt, 0);
554assert_checkalmostequal(norm(coeff(dt / dt0 - 1)), 0, [], del);
555
556qi = [q(2,2) -q(1,2); -q(2,1) q(1,1)] / dt0;
557assert_checkalmostequal(norm(coeff([p 0] * qi - [num den])), 0, [], del);
558assert_checkalmostequal(degree(p), 8, [], 8);
559
560
561