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