001/* 002 * $Id: ComplexRootTest.java 4073 2012-07-28 10:04:29Z kredel $ 003 */ 004 005package edu.jas.root; 006 007 008import java.util.ArrayList; 009import java.util.Collections; 010import java.util.List; 011 012import junit.framework.Test; 013import junit.framework.TestCase; 014import junit.framework.TestSuite; 015 016import edu.jas.arith.BigDecimal; 017import edu.jas.arith.BigRational; 018import edu.jas.kern.ComputerThreads; 019import edu.jas.poly.Complex; 020import edu.jas.poly.ComplexRing; 021import edu.jas.poly.GenPolynomial; 022import edu.jas.poly.GenPolynomialRing; 023import edu.jas.poly.PolyUtil; 024import edu.jas.poly.TermOrder; 025import edu.jas.structure.Power; 026import edu.jas.ufd.Squarefree; 027import edu.jas.ufd.SquarefreeFactory; 028 029 030/** 031 * RootUtil tests with JUnit. 032 * @author Heinz Kredel. 033 */ 034 035public class ComplexRootTest extends TestCase { 036 037 038 /** 039 * main. 040 */ 041 public static void main(String[] args) { 042 //BasicConfigurator.configure(); 043 junit.textui.TestRunner.run(suite()); 044 ComputerThreads.terminate(); 045 } 046 047 048 /** 049 * Constructs a <CODE>ComplexRootTest</CODE> object. 050 * @param name String. 051 */ 052 public ComplexRootTest(String name) { 053 super(name); 054 } 055 056 057 /** 058 */ 059 public static Test suite() { 060 TestSuite suite = new TestSuite(ComplexRootTest.class); 061 return suite; 062 } 063 064 065 TermOrder to = new TermOrder(TermOrder.INVLEX); 066 067 068 GenPolynomialRing<Complex<BigRational>> dfac; 069 070 071 ComplexRing<BigRational> cfac; 072 073 074 BigRational eps; 075 076 077 Complex<BigRational> ceps; 078 079 080 GenPolynomial<Complex<BigRational>> a; 081 082 083 GenPolynomial<Complex<BigRational>> b; 084 085 086 GenPolynomial<Complex<BigRational>> c; 087 088 089 GenPolynomial<Complex<BigRational>> d; 090 091 092 GenPolynomial<Complex<BigRational>> e; 093 094 095 int rl = 1; 096 097 098 int kl = 3; 099 100 101 int ll = 3; 102 103 104 int el = 5; 105 106 107 float q = 0.7f; 108 109 110 @Override 111 protected void setUp() { 112 a = b = c = d = e = null; 113 cfac = new ComplexRing<BigRational>(new BigRational(1)); 114 String[] vars = new String[] { "x" }; 115 dfac = new GenPolynomialRing<Complex<BigRational>>(cfac, rl, to, vars); 116 eps = Power.positivePower(new BigRational(1L, 10L), BigDecimal.DEFAULT_PRECISION); 117 ceps = new Complex<BigRational>(cfac, eps); 118 } 119 120 121 @Override 122 protected void tearDown() { 123 a = b = c = d = e = null; 124 dfac = null; 125 cfac = null; 126 eps = null; 127 } 128 129 130 /** 131 * Test root bound. 132 * 133 */ 134 public void testRootBound() { 135 //a = dfac.random(kl, ll, el, q); 136 a = dfac.univariate(0, 2L).sum(dfac.getONE()); // x^2 + 1 137 //System.out.println("a = " + a); 138 139 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 140 141 Complex<BigRational> Mb = cr.rootBound(a); 142 BigRational M = Mb.getRe(); 143 144 //System.out.println("M = " + M); 145 assertTrue("M >= 1 ", M.compareTo(BigRational.ONE) >= 0); 146 147 //a = a.monic(); 148 a = a.multiply(dfac.fromInteger(5)); 149 //System.out.println("a = " + a); 150 M = cr.rootBound(a).getRe(); 151 152 //System.out.println("M = " + M); 153 assertTrue("M >= 1 ", M.compareTo(BigRational.ONE) >= 0); 154 } 155 156 157 /** 158 * Test Cauchy index. 159 * 160 */ 161 public void testCauchyIndex() { 162 a = dfac.random(kl, ll, el, q); 163 b = dfac.random(kl, ll, el, q); 164 //a = dfac.univariate(0,2L).sum(dfac.getONE()); // x^2 + 1 165 //System.out.println("a = " + a); 166 //System.out.println("b = " + b); 167 168 BigRational l = new BigRational(0); 169 BigRational r = new BigRational(1); 170 GenPolynomialRing<BigRational> fac = new GenPolynomialRing<BigRational>(l, dfac); 171 172 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 173 174 GenPolynomial<BigRational> f = PolyUtil.<BigRational> realPartFromComplex(fac, a); 175 GenPolynomial<BigRational> g = PolyUtil.<BigRational> imaginaryPartFromComplex(fac, b); 176 //System.out.println("re(a) = " + f); 177 //System.out.println("im(b) = " + g); 178 179 long ci = cr.indexOfCauchy(l, r, g, f); 180 //System.out.println("ci = " + ci); 181 182 assertTrue("ci >= 0 ", ci >= -a.degree(0)); 183 } 184 185 186 /** 187 * Test Routh. 188 * 189 */ 190 public void testRouth() { 191 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 192 193 //a = dfac.random(kl, ll, el, q); 194 //b = dfac.random(kl, ll, el, q); 195 Complex<BigRational> I = cfac.getIMAG(); 196 GenPolynomial<Complex<BigRational>> X = dfac.univariate(0); 197 //System.out.println("I = " + I); 198 //System.out.println("X = " + X); 199 200 //a = dfac.univariate(0,2L).sum( dfac.getONE().multiply(I) ); // x^2 + i 201 202 //b = X.subtract( dfac.getONE().multiply( I ) ); // x - i 203 b = X.subtract(dfac.getONE().multiply(I.negate())); // x + i 204 c = X.subtract(dfac.getONE().multiply(I.multiply(cfac.fromInteger(3)))); // x - 3i 205 d = X.subtract(dfac.getONE().multiply(I.multiply(cfac.fromInteger(4)))); // x - 4i 206 e = X.subtract(dfac.getONE().multiply(I.multiply(cfac.fromInteger(5)))); // x - 5i 207 208 a = b.multiply(c).multiply(d).multiply(e); 209 //System.out.println("a = " + a.toScript()); 210 //System.out.println("i = " + cfac.getIMAG()); 211 212 Complex<BigRational> Mb = cr.rootBound(a); 213 BigRational M = Mb.getRe(); 214 //System.out.println("M = " + M); 215 216 BigRational minf = M.negate(); // - infinity 217 BigRational pinf = M; // + infinity 218 GenPolynomialRing<BigRational> fac = new GenPolynomialRing<BigRational>(pinf, dfac); 219 220 GenPolynomial<BigRational> f = PolyUtil.<BigRational> realPartFromComplex(fac, a); 221 GenPolynomial<BigRational> g = PolyUtil.<BigRational> imaginaryPartFromComplex(fac, a); 222 //System.out.println("re(a) = " + f.toScript()); 223 //System.out.println("im(a) = " + g.toScript()); 224 225 long[] ri = cr.indexOfRouth(minf, pinf, f, g); 226 //System.out.println("ri = [" + ri[0] + ", " + ri[1] + " ]"); 227 long deg = ri[0] + ri[1]; 228 assertTrue("sum(ri) == deg(a) ", deg >= a.degree(0)); 229 } 230 231 232 /** 233 * Test winding number. 234 * 235 */ 236 @SuppressWarnings("unchecked") 237 public void testWindingNumber() { 238 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 239 //Complex<BigRational> I = cfac.getIMAG(); 240 241 a = dfac.univariate(0, 2L).sum(cfac.fromInteger(1)); // x^2 + 1 242 //a = dfac.random(kl, ll, el, q); 243 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1 244 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I 245 //a = dfac.univariate(0,1L); // x 246 //System.out.println("a = " + a); 247 248 Complex<BigRational> Mb = cr.rootBound(a); 249 BigRational M = Mb.getRe(); 250 //System.out.println("M = " + M); 251 BigRational eps = new BigRational(1, 1000); 252 //System.out.println("eps = " + eps); 253 254 Complex<BigRational>[] corner = new Complex[4]; 255 256 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw 257 corner[1] = new Complex<BigRational>(cfac, M.negate(), M.negate()); // sw 258 corner[2] = new Complex<BigRational>(cfac, M, M.negate()); // se 259 corner[3] = new Complex<BigRational>(cfac, M, M); // ne 260 261 long v = 0; 262 try { 263 v = cr.windingNumber(new Rectangle<BigRational>(corner), a); 264 } catch (InvalidBoundaryException e) { 265 fail("" + e); 266 } 267 //System.out.println("winding number = " + v); 268 assertTrue("wind(rect,a) == 2 ", v == 2); 269 270 //if ( true ) return; 271 272 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw 273 corner[1] = new Complex<BigRational>(cfac, M.negate(), eps); // sw 274 corner[2] = new Complex<BigRational>(cfac, M, eps); // se 275 corner[3] = new Complex<BigRational>(cfac, M, M); // ne 276 277 try { 278 v = cr.windingNumber(new Rectangle<BigRational>(corner), a); 279 } catch (InvalidBoundaryException e) { 280 fail("" + e); 281 } 282 //System.out.println("winding number = " + v); 283 assertTrue("wind(rect,a) == 1 ", v == 1); 284 285 corner[0] = new Complex<BigRational>(cfac, eps.negate(), eps); // nw 286 corner[1] = new Complex<BigRational>(cfac, eps.negate(), eps.negate()); // sw 287 corner[2] = new Complex<BigRational>(cfac, eps, eps.negate()); // se 288 corner[3] = new Complex<BigRational>(cfac, eps, eps); // ne 289 290 try { 291 v = cr.windingNumber(new Rectangle<BigRational>(corner), a); 292 } catch (InvalidBoundaryException e) { 293 fail("" + e); 294 } 295 //System.out.println("winding number = " + v); 296 assertTrue("wind(rect,a) == 0 ", v == 0); 297 } 298 299 300 /** 301 * Test complex roots, sqrt(-1). 302 * 303 */ 304 @SuppressWarnings("unchecked") 305 public void testComplexRootsImag() { 306 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 307 //Complex<BigRational> I = cfac.getIMAG(); 308 309 a = dfac.univariate(0, 2L).sum(cfac.fromInteger(1)); // x^2 + 1 310 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1 311 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I 312 //a = dfac.univariate(0,1L); // x 313 //System.out.println("a = " + a); 314 315 Complex<BigRational> Mb = cr.rootBound(a); 316 BigRational M = Mb.getRe(); 317 //System.out.println("M = " + M); 318 319 Complex<BigRational>[] corner = new Complex[4]; 320 321 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw 322 corner[1] = new Complex<BigRational>(cfac, M.negate(), M.negate()); // sw 323 corner[2] = new Complex<BigRational>(cfac, M, M.negate()); // se 324 corner[3] = new Complex<BigRational>(cfac, M, M); // ne 325 326 Rectangle<BigRational> rect = new Rectangle<BigRational>(corner); 327 328 List<Rectangle<BigRational>> roots = null; 329 try { 330 roots = cr.complexRoots(rect, a); 331 } catch (InvalidBoundaryException e) { 332 fail("" + e); 333 } 334 //System.out.println("roots = " + roots); 335 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 336 } 337 338 339 /** 340 * Test complex roots. 341 */ 342 @SuppressWarnings("unchecked") 343 public void testComplexRootsRand() { 344 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 345 //Complex<BigRational> I = cfac.getIMAG(); 346 347 a = dfac.random(kl, ll, el, q); 348 Squarefree<Complex<BigRational>> engine = SquarefreeFactory 349 .<Complex<BigRational>> getImplementation(cfac); 350 a = engine.squarefreePart(a); 351 352 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1 353 //a = dfac.univariate(0,2L).sum(cfac.fromInteger(1)); // x^2 + 1 354 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I 355 //a = dfac.univariate(0,1L); // x 356 //System.out.println("a = " + a); 357 358 Complex<BigRational> Mb = cr.rootBound(a); 359 BigRational M = Mb.getRe(); 360 //System.out.println("M = " + M); 361 362 Complex<BigRational>[] corner = new Complex[4]; 363 364 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw 365 corner[1] = new Complex<BigRational>(cfac, M.negate(), M.negate()); // sw 366 corner[2] = new Complex<BigRational>(cfac, M, M.negate()); // se 367 corner[3] = new Complex<BigRational>(cfac, M, M); // ne 368 369 Rectangle<BigRational> rect = new Rectangle<BigRational>(corner); 370 371 372 List<Rectangle<BigRational>> roots = null; 373 try { 374 roots = cr.complexRoots(rect, a); 375 } catch (InvalidBoundaryException e) { 376 fail("" + e); 377 } 378 //System.out.println("a = " + a); 379 //System.out.println("roots = " + roots); 380 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 381 } 382 383 384 /** 385 * Test complex roots. 386 */ 387 public void testComplexRoots() { 388 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 389 390 a = dfac.random(kl, ll, el + 1, q); 391 //System.out.println("a = " + a); 392 393 List<Rectangle<BigRational>> roots = cr.complexRoots(a); 394 //System.out.println("a = " + a); 395 //System.out.println("roots = " + roots); 396 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 397 } 398 399 400 /** 401 * Test complex root refinement. 402 */ 403 public void testComplexRootRefinement() { 404 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 405 406 a = dfac.random(kl, ll, el - 1, q); 407 //a = dfac.parse("( (x-1)^3 )"); 408 Squarefree<Complex<BigRational>> engine = SquarefreeFactory 409 .<Complex<BigRational>> getImplementation(cfac); 410 //System.out.println("a = " + a); 411 a = engine.squarefreePart(a); 412 //System.out.println("a = " + a); 413 414 List<Rectangle<BigRational>> roots = cr.complexRoots(a); 415 //System.out.println("a = " + a); 416 //System.out.println("roots = " + roots); 417 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 418 419 BigRational len = new BigRational(1, 1000); 420 //System.out.println("len = " + len); 421 422 for (Rectangle<BigRational> root : roots) { 423 try { 424 Rectangle<BigRational> refine = cr.complexRootRefinement(root, a, len); 425 //System.out.println("refine = " + refine); 426 assertFalse("refine != null", refine == null); 427 } catch (InvalidBoundaryException e) { 428 fail("" + e); 429 } 430 } 431 } 432 433 434 /** 435 * Test complex root refinement full. 436 */ 437 public void testComplexRootRefinementFull() { 438 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 439 440 a = dfac.random(kl, ll, el - 1, q); 441 //a = dfac.parse("( (x-1)^3 )"); 442 //a = dfac.parse("( x^4-2 )"); 443 //System.out.println("a = " + a); 444 445 BigRational len = new BigRational(1, 1000); 446 //System.out.println("len = " + len); 447 448 List<Rectangle<BigRational>> refine = cr.complexRoots(a, len); 449 //System.out.println("refine = " + refine); 450 assertTrue("#roots == deg(a) ", refine.size() == a.degree(0)); 451 } 452 453 454 /** 455 * Test winding number with wrong precondition. 456 */ 457 @SuppressWarnings("unchecked") 458 public void testWindingNumberWrong() { 459 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 460 //Complex<BigRational> I = cfac.getIMAG(); 461 462 a = dfac.univariate(0, 2L).sum(cfac.fromInteger(1)); // x^2 + 1 463 //a = dfac.random(kl, ll, el, q); 464 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1 465 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I 466 //a = dfac.univariate(0,1L); // x 467 //System.out.println("a = " + a); 468 469 Complex<BigRational> Mb = cfac.fromInteger(1); //.divide(cfac.fromInteger(2)); //cr.rootBound(a); 470 BigRational M = Mb.getRe(); 471 //System.out.println("M = " + M); 472 //BigRational eps = new BigRational(1, 1000); 473 //System.out.println("eps = " + eps); 474 BigRational zero = new BigRational(); 475 //BigRational one = new BigRational(1); 476 477 Complex<BigRational>[] corner = new Complex[4]; 478 479 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw 480 corner[1] = new Complex<BigRational>(cfac, M.negate(), zero); // sw 481 corner[2] = new Complex<BigRational>(cfac, M, zero); // se 482 corner[3] = new Complex<BigRational>(cfac, M, M); // ne 483 484 Rectangle<BigRational> rect = new Rectangle<BigRational>(corner); 485 //System.out.println("rect = " + rect); 486 487 try { 488 long v = cr.windingNumber(rect, a); 489 System.out.println("winding number = " + v); 490 fail("wind(rect,a) must throw an exception"); 491 } catch (InvalidBoundaryException e) { 492 } 493 } 494 495 496 /** 497 * Test complex root approximation. 498 */ 499 public void testComplexRootApproximation() { 500 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 501 502 //a = dfac.random(kl, ll, el-1, q); 503 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*(x - { 0i1 })*(x-5) )*( x^4-2 )"); 504 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*( x^4-2 ) )"); 505 //a = dfac.parse("( (x-2)*( x^4-2 ) )"); 506 a = dfac.parse("( ( x^4-2 ) )"); 507 b = dfac.parse("( (x-1)*(x-2)*(x-3) )"); 508 c = dfac.parse("( x^4-2 )"); 509 d = dfac.parse("( (x - { 0i1 })*(x-5) )"); 510 //a = c; // b; //.multiply(c); //.multiply(d); 511 //System.out.println("a = " + a); 512 //System.out.println("b = " + b); 513 //System.out.println("c = " + c); 514 //System.out.println("d = " + d); 515 //a = b.multiply(c).multiply(d); 516 //System.out.println("a = " + a); 517 Squarefree<Complex<BigRational>> engine = SquarefreeFactory 518 .<Complex<BigRational>> getImplementation(cfac); 519 a = engine.squarefreePart(a); 520 //System.out.println("a = " + a); 521 522 eps = eps.multiply(new BigRational(1000000)); 523 //System.out.println("eps = " + eps); 524 //BigDecimal eps1 = new BigDecimal(eps); 525 //BigDecimal eps2 = eps1.multiply(new BigDecimal("10")); 526 //System.out.println("eps1 = " + eps1); 527 //System.out.println("eps2 = " + eps2); 528 529 List<Rectangle<BigRational>> roots = cr.complexRoots(a); 530 //System.out.println("a = " + a); 531 //System.out.println("roots = " + roots); 532 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 533 534 for (Rectangle<BigRational> root : roots) { 535 try { 536 Complex<BigDecimal> cd = cr.approximateRoot(root, a, eps); 537 //System.out.println("cd = " + cd); 538 assertFalse("cd != 0", cd.isZERO()); 539 } catch (NoConvergenceException e) { 540 //fail("" + e); 541 } 542 } 543 } 544 545 546 /** 547 * Test complex root approximation full algorithm. 548 */ 549 public void testComplexRootApproximationFull() { 550 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 551 552 //a = dfac.random(kl, ll, el-1, q); 553 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*(x - { 0i1 })*(x-5) )*( x^4-2 )"); 554 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*( x^4-2 ) )"); 555 a = dfac.parse("( (x-2)*( x^4-2 ) )"); 556 //a = dfac.parse("( ( x^4-2 ) )"); 557 b = dfac.parse("( (x-1)*(x-2)*(x-3) )"); 558 c = dfac.parse("( x^4-2 )"); 559 d = dfac.parse("( (x - { 0i1 })*(x-5) )"); 560 //a = c; // b; //.multiply(c); //.multiply(d); 561 //System.out.println("a = " + a); 562 //System.out.println("b = " + b); 563 //System.out.println("c = " + c); 564 //System.out.println("d = " + d); 565 //a = b.multiply(c).multiply(d); 566 //System.out.println("a = " + a); 567 568 eps = eps.multiply(new BigRational(1000000)); 569 //System.out.println("eps = " + eps); 570 //BigDecimal eps1 = new BigDecimal(eps); 571 //BigDecimal eps2 = eps1.multiply(new BigDecimal("10")); 572 //System.out.println("eps1 = " + eps1); 573 //System.out.println("eps2 = " + eps2); 574 575 List<Complex<BigDecimal>> roots = cr.approximateRoots(a, eps); 576 //System.out.println("a = " + a); 577 //System.out.println("roots = " + roots); 578 //now always true: 579 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 580 } 581 582 583 /** 584 * Test complex root approximation full algorithm with Wilkinson 585 * polynomials. p = (x-i0)*(x-i1)*(x-i2)*(x-i3*...*(x-in) 586 */ 587 public void testComplexRootApproximationWilkinsonFull() { 588 final int N = 4; 589 d = dfac.getONE(); 590 e = dfac.univariate(0); 591 592 BigDecimal br = new BigDecimal(); 593 ComplexRing<BigDecimal> cf = new ComplexRing<BigDecimal>(br); 594 Complex<BigDecimal> I = cf.getIMAG(); 595 Complex<BigDecimal> cc = null; 596 Complex<BigRational> Ir = cfac.getIMAG(); 597 598 List<Complex<BigDecimal>> Rn = new ArrayList<Complex<BigDecimal>>(N); 599 a = d; 600 for (int i = 0; i < N; i++) { 601 cc = cf.fromInteger(i).multiply(I); 602 Rn.add(cc); 603 c = dfac.fromInteger(i).multiply(Ir); 604 b = e.subtract(c); 605 a = a.multiply(b); 606 } 607 //System.out.println("a = " + a); 608 Collections.reverse(Rn); 609 //System.out.println("Rn = " + Rn); 610 611 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 612 613 eps = eps.multiply(new BigRational(100000)); 614 //System.out.println("eps = " + eps); 615 BigDecimal eps1 = new BigDecimal(eps); 616 BigDecimal eps2 = eps1.multiply(new BigDecimal("10")); 617 //System.out.println("eps1 = " + eps1); 618 //System.out.println("eps2 = " + eps2); 619 620 List<Complex<BigDecimal>> roots = cr.approximateRoots(a, eps); 621 //System.out.println("a = " + a); 622 //System.out.println("roots = " + roots); 623 //now always true: 624 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 625 626 int i = 0; 627 for (Complex<BigDecimal> dd : roots) { 628 Complex<BigDecimal> di = Rn.get(i++); 629 //System.out.print("di = " + di + ", "); 630 //System.out.println("dd = " + dd); 631 assertTrue("|dd - di| < eps ", dd.subtract(di).norm().getRe().compareTo(eps2) <= 0); 632 } 633 } 634 635 636 /** 637 * Test complex root approximation full algorithm with Wilkinson 638 * polynomials, inverse roots. p = (x-1/i1)*(x-1/i2)*(x-1/i3*...*(x-1/in) 639 */ 640 public void testComplexRootApproximationWilkinsonInverseFull() { 641 final int N = 5; 642 d = dfac.getONE(); 643 e = dfac.univariate(0); 644 645 BigDecimal br = new BigDecimal(); 646 ComplexRing<BigDecimal> cf = new ComplexRing<BigDecimal>(br); 647 Complex<BigDecimal> I = cf.getIMAG(); 648 Complex<BigDecimal> cc = null; 649 Complex<BigRational> Ir = cfac.getIMAG(); 650 651 List<Complex<BigDecimal>> Rn = new ArrayList<Complex<BigDecimal>>(N); 652 a = d; 653 for (int i = 1; i < N; i++) { 654 cc = cf.fromInteger(i).multiply(I); 655 cc = cc.inverse(); 656 Rn.add(cc); 657 c = dfac.fromInteger(i).multiply(Ir); 658 c = d.divide(c); 659 b = e.subtract(c); 660 a = a.multiply(b); 661 } 662 //System.out.println("a = " + a); 663 //Collections.sort(Rn); 664 //System.out.println("Rn = " + Rn); 665 666 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 667 668 eps = eps.multiply(new BigRational(100000)); 669 //System.out.println("eps = " + eps); 670 BigDecimal eps1 = new BigDecimal(eps); 671 BigDecimal eps2 = eps1.multiply(new BigDecimal("10")); 672 //System.out.println("eps1 = " + eps1); 673 //System.out.println("eps2 = " + eps2); 674 675 List<Complex<BigDecimal>> roots = cr.approximateRoots(a, eps); 676 //System.out.println("a = " + a); 677 //System.out.println("roots = " + roots); 678 //now always true: 679 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 680 //Collections.sort(roots); 681 //System.out.println("roots = " + roots); 682 683 for (Complex<BigDecimal> dd : roots) { 684 //System.out.println("dd = " + dd); 685 boolean t = false; 686 for (Complex<BigDecimal> di : Rn) { 687 //System.out.println("di = " + di); 688 t = dd.subtract(di).norm().getRe().compareTo(eps2) <= 0; 689 if (t) { 690 break; 691 } 692 } 693 if (!t) { 694 //assertTrue("|dd - di| < eps ", dd.subtract(di).norm().getRe().compareTo(eps2) <= 0); 695 fail("|dd - di| < eps "); 696 } 697 } 698 } 699 700 701 /** 702 * Test complex root invariant rectangle. 703 */ 704 public void testComplexRootInvariant() { 705 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 706 707 a = dfac.random(kl, ll, el - 1, q); 708 b = dfac.random(kl, ll, 2, q); 709 //a = dfac.parse("( (x-1)^3 )"); 710 //a = dfac.parse("( x^4-2 )"); 711 if (a.degree() == 0) { 712 return; 713 } 714 Squarefree<Complex<BigRational>> engine = SquarefreeFactory 715 .<Complex<BigRational>> getImplementation(cfac); 716 a = engine.squarefreePart(a); 717 b = engine.squarefreePart(b); 718 //System.out.println("a = " + a); 719 //System.out.println("b = " + b); 720 721 List<Rectangle<BigRational>> roots = cr.complexRoots(a); 722 //System.out.println("roots = " + roots); 723 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 724 725 Rectangle<BigRational> rect = roots.get(0); 726 //System.out.println("rect = " + rect); 727 728 try { 729 Rectangle<BigRational> ref = cr.invariantRectangle(rect, a, b); 730 //System.out.println("rect = " + rect); 731 //System.out.println("ref = " + ref); 732 assertTrue("rect subseteq ref ", rect.contains(ref)); 733 } catch (InvalidBoundaryException e) { 734 e.printStackTrace(); 735 fail("bad boundary"); 736 } 737 } 738 739 740 /** 741 * Test complex root invariant magnitude rectangle. 742 */ 743 public void testComplexRootInvariantMagnitude() { 744 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac); 745 746 a = dfac.random(kl, ll, el - 1, q); 747 b = dfac.random(kl, ll, 3, q); 748 //a = dfac.parse("( x^2 + 1 )"); 749 //b = dfac.parse("( x - 0i1 )"); 750 if (a.degree() == 0) { 751 return; 752 } 753 Squarefree<Complex<BigRational>> engine = SquarefreeFactory 754 .<Complex<BigRational>> getImplementation(cfac); 755 a = engine.squarefreePart(a); 756 b = engine.squarefreePart(b); 757 //System.out.println("a = " + a); 758 //System.out.println("b = " + b); 759 760 List<Rectangle<BigRational>> roots = cr.complexRoots(a); 761 //System.out.println("roots = " + roots); 762 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0)); 763 764 Rectangle<BigRational> rect = roots.get(0); 765 //System.out.println("rect = " + rect); 766 767 try { 768 Rectangle<BigRational> ref = cr.invariantMagnitudeRectangle(rect, a, b, eps); 769 //System.out.println("ref = " + ref); 770 assertTrue("rect subseteq ref ", rect.contains(ref)); 771 Complex<BigRational> mag = cr.complexRectangleMagnitude(ref, a, b); 772 //System.out.println("mag = " + mag); 773 Complex<BigRational> cmag = cr.complexMagnitude(ref, a, b, eps); 774 //System.out.println("cmag = " + cmag); 775 assertEquals("mag == cmag: " + cmag, mag, cmag); 776 //BigRational rmag = cmag.getRe(); 777 //System.out.println("rmag = " + new BigDecimal(cmag.getRe()) + " i " + new BigDecimal(cmag.getIm())); 778 } catch (InvalidBoundaryException e) { 779 e.printStackTrace(); 780 fail("bad boundary"); 781 } 782 } 783 784}