001/* 002 * $Id$ 003 */ 004 005package edu.jas.root; 006 007 008import java.util.ArrayList; 009import java.util.Collections; 010import java.util.List; 011 012import edu.jas.arith.BigDecimal; 013import edu.jas.arith.BigRational; 014import edu.jas.arith.Roots; 015import edu.jas.poly.GenPolynomial; 016import edu.jas.poly.GenPolynomialRing; 017import edu.jas.poly.TermOrder; 018import edu.jas.structure.Power; 019 020import junit.framework.Test; 021import junit.framework.TestCase; 022import junit.framework.TestSuite; 023 024 025/** 026 * RealRoot tests with JUnit. 027 * @author Heinz Kredel 028 */ 029 030public class RealRootTest extends TestCase { 031 032 033 /** 034 * main. 035 */ 036 public static void main(String[] args) { 037 junit.textui.TestRunner.run(suite()); 038 } 039 040 041 /** 042 * Constructs a <CODE>RealRootTest</CODE> object. 043 * @param name String. 044 */ 045 public RealRootTest(String name) { 046 super(name); 047 } 048 049 050 /** 051 */ 052 public static Test suite() { 053 TestSuite suite = new TestSuite(RealRootTest.class); 054 return suite; 055 } 056 057 058 TermOrder to = new TermOrder(TermOrder.INVLEX); 059 060 061 GenPolynomialRing<BigRational> dfac; 062 063 064 BigRational ai, bi, ci, di, ei, eps; 065 066 067 GenPolynomial<BigRational> a, b, c, d, e; 068 069 070 int rl = 1; 071 072 073 int kl = 5; 074 075 076 int ll = 7; 077 078 079 int el = 7; 080 081 082 float q = 0.7f; 083 084 085 @Override 086 protected void setUp() { 087 a = b = c = d = e = null; 088 ai = bi = ci = di = ei = null; 089 String[] vars = new String[] { "x" }; 090 dfac = new GenPolynomialRing<BigRational>(new BigRational(1), rl, to, vars); 091 // eps = new BigRational(1L,1000000L*1000000L*1000000L); 092 eps = Power.positivePower(new BigRational(1L, 10L), BigDecimal.DEFAULT_PRECISION); 093 } 094 095 096 @Override 097 protected void tearDown() { 098 a = b = c = d = e = null; 099 ai = bi = ci = di = ei = null; 100 dfac = null; 101 eps = null; 102 } 103 104 105 /** 106 * Test Sturm sequence. 107 */ 108 public void testSturmSequence() { 109 a = dfac.random(kl, ll, el, q); 110 //System.out.println("a = " + a); 111 112 RealRootsSturm<BigRational> rrs = new RealRootsSturm<BigRational>(); 113 114 List<GenPolynomial<BigRational>> S = rrs.sturmSequence(a); 115 //System.out.println("S = " + S); 116 117 try { 118 b = a.remainder(S.get(0)); 119 } catch (Exception e) { 120 fail("not S(0)|f " + e); 121 } 122 assertTrue("a mod S(0) == 0 ", b.isZERO()); 123 124 assertTrue("S(-1) == 1 ", S.get(S.size() - 1).isConstant()); 125 } 126 127 128 /** 129 * Test root bound. 130 */ 131 public void testRootBound() { 132 a = dfac.random(kl, ll, el, q); 133 //System.out.println("a = " + a); 134 135 RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>(); 136 137 // used root bound 138 BigRational M = rr.realRootBound(a); 139 //System.out.println("M = " + M); 140 assertTrue("M >= 1 ", M.compareTo(BigRational.ONE) >= 0); 141 Interval<BigRational> v1 = new Interval<BigRational>(M.negate(), M); 142 long r1 = rr.realRootCount(v1, a); 143 //System.out.println("v1 = " + v1 + ", r1 = " + r1); 144 145 a = a.monic(); 146 //System.out.println("a = " + a); 147 BigDecimal ar = M.getDecimal(); 148 //System.out.println("ar = " + ar); 149 assertTrue("ar >= 1 ", ar.compareTo(BigDecimal.ONE) >= 0); 150 151 // maxNorm root bound 152 BigRational mr = a.maxNorm().getRational().sum(BigRational.ONE); 153 BigDecimal dr = mr.getDecimal(); 154 //System.out.println("dr = " + dr); 155 //assertTrue("ar >= maxNorm(a): " + (ar.subtract(dr)), ar.compareTo(dr) >= 0); 156 Interval<BigRational> v2 = new Interval<BigRational>(mr.negate(), mr); 157 long r2 = rr.realRootCount(v2, a); 158 //System.out.println("v2 = " + v2 + ", r2 = " + r2); 159 assertTrue("r1 == r2: " + (r2 - r1), r1 == r2); 160 161 // squareNorm root bound 162 BigRational qr = a.squareNorm().getRational(); 163 BigDecimal ir = Roots.sqrt(qr.getDecimal()); 164 //qr = Roots.sqrt(qr); 165 //System.out.println("ir = " + ir); 166 //assertTrue("ar >= squareNorm(a): " + (ar.subtract(ir)), ar.compareTo(ir) >= 0); 167 Interval<BigRational> v3 = new Interval<BigRational>(qr.negate(), qr); 168 long r3 = rr.realRootCount(v3, a); 169 //System.out.println("v3 = " + v3 + ", r3 = " + r3); 170 assertTrue("r1 == r3: " + (r3 - r1), r1 == r3); 171 172 // sumNorm root bound 173 BigRational pr = a.sumNorm().getRational(); 174 BigDecimal sr = pr.getDecimal(); 175 //System.out.println("sr = " + sr); 176 //assertTrue("ar >= squareNorm(a): " + (ar.subtract(sr)), ar.compareTo(sr) >= 0); 177 Interval<BigRational> v4 = new Interval<BigRational>(pr.negate(), pr); 178 long r4 = rr.realRootCount(v4, a); 179 //System.out.println("v4 = " + v4 + ", r4 = " + r4); 180 assertTrue("r1 == r4: " + (r4 - r1), r1 == r4); 181 182 // minimal root bound 183 BigDecimal dri = dr.sum(BigDecimal.ONE).inverse(); 184 //System.out.println("dri = " + dri + ", sign(dri) = " + dri.signum()); 185 assertTrue("minimal root > 0: " + dri, dri.signum() > 0); 186 BigDecimal mri = rr.realMinimalRootBound(a).getDecimal(); 187 //System.out.println("mri = " + mri + ", sign(mri) = " + mri.signum()); 188 BigDecimal s = dri.subtract(mri).abs(); 189 eps = eps.multiply(BigRational.ONE.fromInteger(10)); 190 //System.out.println("s = " + s + ", eps = " + eps.getDecimal()); 191 assertTrue("minimal root: " + dri, s.compareTo(eps.getDecimal()) < 0); 192 193 // minimal root separation 194 long n = a.degree(); 195 if (n > 0) { 196 BigDecimal sep = sr.sum(BigDecimal.ONE).power(2 * n).multiply(sr.fromInteger(n).power(n + 1)) 197 .inverse(); 198 //System.out.println("sep = " + sep + ", sign(sep) = " + sep.signum()); 199 assertTrue("separation(a) > 0: " + sep, sep.signum() > 0); 200 BigDecimal sri = rr.realMinimalRootSeparation(a).getDecimal(); 201 BigDecimal ss = sep.subtract(sri).abs(); 202 assertTrue("minimal separation: " + sep, ss.compareTo(eps.getDecimal()) < 0); 203 } 204 } 205 206 207 /** 208 * Test real root isolation. 209 */ 210 public void testRealRootIsolation() { 211 a = dfac.random(kl, ll * 2, el * 2, q); 212 //a = a.multiply( dfac.univariate(0) ); 213 //System.out.println("a = " + a); 214 215 RealRoots<BigRational> rr = new RealRootsSturm<BigRational>(); 216 217 List<Interval<BigRational>> R = rr.realRoots(a); 218 //System.out.println("R = " + R); 219 //assertTrue("#roots >= 0 ", R.size() >= 0); 220 assertTrue("#roots >= 0 ", R != null); 221 } 222 223 224 /** 225 * Test Thom lemma real root sign sequence. 226 */ 227 public void testRealRootSignSequence() { 228 a = dfac.random(kl, ll * 2, el * 2, q); 229 if (a.degree() % 2 == 0) { 230 a = a.multiply(dfac.univariate(0).subtract(dfac.getONE())); 231 } 232 if (a.trailingExpVector().degree() > 0) { 233 a = a.subtract(dfac.getONE()); // exclude root 0 234 } 235 //System.out.println("a = " + a); 236 RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>(); 237 238 List<Interval<BigRational>> R = rr.realRoots(a); 239 //System.out.println("R = " + R); 240 //assertTrue("#roots >= 0 ", R.size() >= 0); 241 assertTrue("#roots >= 0 ", R != null); 242 243 int l = R.size(); 244 Interval<BigRational> v = R.get(l - 1); 245 Interval<BigRational> u = R.get(0); 246 if (u.left.isZERO() && u.right.isZERO()) { 247 Interval<BigRational> w = v; 248 v = u; 249 u = w; 250 } 251 Interval<BigRational> vm = new Interval<BigRational>(u.left, v.right); 252 //System.out.println("v = " + v); 253 //System.out.println("u = " + u); 254 //System.out.println("vm = " + vm); 255 long rc = rr.realRootCount(vm, a); 256 //System.out.println("rc = " + rc); 257 assertTrue("root number: R = " + R + ", a = " + a + ", rc = " + rc, rc == l); 258 long rn = rr.realRootNumber(a, vm); 259 assertTrue("root number == " + rn, rn == l); 260 261 long d = a.degree(); 262 List<GenPolynomial<BigRational>> fs = rr.fourierSequence(a); 263 //System.out.println("fs = " + fs); 264 assertTrue("len(fs) == " + (d + 1 - fs.size()), fs.size() == (d + 1)); 265 266 //List<Integer> ss = rr.signSequence(a, v); 267 //System.out.println("ss = " + ss); 268 //assertTrue("len(ss) == " + (d-ss.size()), ss.size() == d); 269 for (Interval<BigRational> t : R) { 270 List<Integer> ss = rr.signSequence(a, t); 271 //System.out.println("ss = " + ss); 272 assertTrue("len(ss) == " + (d - ss.size()), ss.size() == d); 273 } 274 } 275 276 277 /** 278 * Test real root isolation Wilkinson polynomials. p = 279 * (x-0)*(x-1)*(x-2)*(x-3)*...*(x-n) 280 */ 281 public void testRealRootIsolationWilkinson() { 282 final int N = 10; 283 d = dfac.getONE(); 284 e = dfac.univariate(0); 285 286 List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N); 287 a = d; 288 for (int i = 0; i < N; i++) { 289 c = dfac.fromInteger(i); 290 Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient())); 291 b = e.subtract(c); 292 a = a.multiply(b); 293 } 294 //System.out.println("a = " + a); 295 296 RealRoots<BigRational> rr = new RealRootsSturm<BigRational>(); 297 298 List<Interval<BigRational>> R = rr.realRoots(a); 299 //System.out.println("R = " + R); 300 301 assertTrue("#roots = " + N + " ", R.size() == N); 302 303 eps = eps.multiply(new BigRational("1/10")); 304 //System.out.println("eps = " + eps); 305 306 R = rr.refineIntervals(R, a, eps); 307 //System.out.println("R = " + R); 308 int i = 0; 309 for (Interval<BigRational> v : R) { 310 BigDecimal dd = v.toDecimal(); 311 BigDecimal di = Rn.get(i++).toDecimal(); 312 //System.out.println("v = " + dd); 313 //System.out.println("vi = " + di); 314 //System.out.println("|dd - di| < eps: " + dd.compareTo(di)); 315 assertTrue("|dd - di| < eps ", dd.compareTo(di) == 0); 316 } 317 } 318 319 320 /** 321 * Test real root isolation Wilkinson polynomials inverse. p = 322 * (x-1)*(x-1/2)*(x-1/3)*...*(x-1/n) 323 */ 324 public void testRealRootIsolationWilkinsonInverse() { 325 final int N = 9; 326 d = dfac.getONE(); 327 e = dfac.univariate(0); 328 329 List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N); 330 a = d; 331 for (int i = 1; i < N; i++) { // use only for i > 0, since reverse 332 c = dfac.fromInteger(i); 333 if (i != 0) { 334 c = d.divide(c); 335 } 336 Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient())); 337 b = e.subtract(c); 338 a = a.multiply(b); 339 } 340 //System.out.println("a = " + a); 341 //System.out.println("Rn = " + Rn); 342 Collections.reverse(Rn); 343 //System.out.println("Rn = " + Rn); 344 345 RealRoots<BigRational> rr = new RealRootsSturm<BigRational>(); 346 347 List<Interval<BigRational>> R = rr.realRoots(a); 348 //System.out.println("R = " + R); 349 350 assertTrue("#roots = " + (N - 1) + " ", R.size() == (N - 1)); 351 352 eps = eps.multiply(new BigRational("1/100")); 353 //System.out.println("eps = " + eps); 354 355 R = rr.refineIntervals(R, a, eps); 356 //System.out.println("R = " + R); 357 int i = 0; 358 for (Interval<BigRational> v : R) { 359 BigDecimal dd = v.toDecimal(); //.sum(eps1); 360 BigDecimal di = Rn.get(i++).toDecimal(); 361 //System.out.println("v = " + dd); 362 //System.out.println("vi = " + di); 363 //System.out.println("|dd - di| < eps: " + dd.compareTo(di)); 364 assertTrue("|dd - di| < eps ", dd.compareTo(di) == 0); 365 } 366 } 367 368 369 /** 370 * Test real algebraic number sign. 371 */ 372 public void testRealAlgebraicNumberSign() { 373 d = dfac.fromInteger(2); 374 e = dfac.univariate(0); 375 376 a = e.multiply(e); 377 // a = a.multiply(e).multiply(e).multiply(e); 378 a = a.subtract(d); // x^2 -2 379 //System.out.println("a = " + a); 380 381 RealRoots<BigRational> rr = new RealRootsSturm<BigRational>(); 382 383 ai = new BigRational(1); 384 bi = new BigRational(2); 385 Interval<BigRational> iv = new Interval<BigRational>(ai, bi); 386 //System.out.println("iv = " + iv); 387 assertTrue("sign change", rr.signChange(iv, a)); 388 389 b = dfac.random(kl, (int) a.degree() + 1, (int) a.degree(), 1.0f); 390 //b = dfac.getZERO(); 391 //b = dfac.random(kl,ll,el,q); 392 //b = b.multiply(b); 393 //b = b.abs().negate(); 394 //System.out.println("b = " + b); 395 if (b.isZERO()) { 396 int s = rr.realSign(iv, a, b); 397 assertTrue("algebraic sign", s == 0); 398 return; 399 } 400 401 int as = rr.realSign(iv, a, b); 402 //System.out.println("as = " + as); 403 // how to test? 404 int asn = rr.realSign(iv, a, b.negate()); 405 //System.out.println("asn = " + asn); 406 assertTrue("algebraic sign", as != asn); 407 408 iv = new Interval<BigRational>(bi.negate(), ai.negate()); 409 //System.out.println("iv = " + iv); 410 assertTrue("sign change", rr.signChange(iv, a)); 411 412 int as1 = rr.realSign(iv, a, b); 413 //System.out.println("as1 = " + as1); 414 // how to test? 415 int asn1 = rr.realSign(iv, a, b.negate()); 416 //System.out.println("asn1 = " + asn1); 417 assertTrue("algebraic sign", as1 != asn1); 418 419 assertTrue("algebraic sign", as * as1 == asn * asn1); 420 } 421 422 423 /** 424 * Test real root isolation and decimal refinement of Wilkinson polynomials. 425 * p = (x-0)*(x-1)*(x-2)*(x-3)*...*(x-n) 426 */ 427 public void testRealRootIsolationDecimalWilkinson() { 428 final int N = 10; 429 d = dfac.getONE(); 430 e = dfac.univariate(0); 431 432 List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N); 433 a = d; 434 for (int i = 0; i < N; i++) { 435 c = dfac.fromInteger(i); 436 Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient())); 437 b = e.subtract(c); 438 a = a.multiply(b); 439 } 440 //System.out.println("a = " + a); 441 442 RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>(); 443 444 List<Interval<BigRational>> R = rr.realRoots(a); 445 //System.out.println("R = " + R); 446 447 assertTrue("#roots = " + N + " ", R.size() == N); 448 449 eps = eps.multiply(new BigRational(100000)); 450 //System.out.println("eps = " + eps); 451 BigDecimal eps1 = new BigDecimal(eps); 452 BigDecimal eps2 = eps1.multiply(new BigDecimal("100")); 453 //System.out.println("eps1 = " + eps1); 454 //System.out.println("eps2 = " + eps2); 455 456 try { 457 int i = 0; 458 for (Interval<BigRational> v : R) { 459 //System.out.println("v = " + v); 460 BigDecimal dd = rr.approximateRoot(v, a, eps); 461 BigDecimal di = Rn.get(i++).toDecimal(); 462 //System.out.println("di = " + di); 463 //System.out.println("dd = " + dd); 464 assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0); 465 } 466 } catch (NoConvergenceException e) { 467 fail(e.toString()); 468 } 469 } 470 471 472 /** 473 * Test real root isolation and decimal refinement of Wilkinson polynomials, 474 * inverse roots. p = (x-1)*(x-1/2)*(x-1/3)*...*(x-1/n) 475 */ 476 public void testRealRootIsolationDecimalWilkinsonInverse() { 477 final int N = 10; 478 d = dfac.getONE(); 479 e = dfac.univariate(0); 480 481 List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N); 482 a = d; 483 for (int i = 1; i < N; i++) { // use only for i > 0, since reverse 484 c = dfac.fromInteger(i); 485 if (i != 0) { 486 c = d.divide(c); 487 } 488 Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient())); 489 b = e.subtract(c); 490 a = a.multiply(b); 491 } 492 //System.out.println("a = " + a); 493 //System.out.println("Rn = " + Rn); 494 Collections.reverse(Rn); 495 //System.out.println("Rn = " + Rn); 496 497 RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>(); 498 499 List<Interval<BigRational>> R = rr.realRoots(a); 500 //System.out.println("R = " + R); 501 502 assertTrue("#roots = " + (N - 1) + " ", R.size() == (N - 1)); 503 504 eps = eps.multiply(new BigRational(1000000)); 505 //System.out.println("eps = " + eps); 506 BigDecimal eps1 = new BigDecimal(eps); 507 BigDecimal eps2 = eps1.multiply(new BigDecimal("10")); 508 //System.out.println("eps1 = " + eps1); 509 //System.out.println("eps2 = " + eps2); 510 511 try { 512 int i = 0; 513 for (Interval<BigRational> v : R) { 514 //System.out.println("v = " + v); 515 BigDecimal dd = rr.approximateRoot(v, a, eps); 516 BigDecimal di = Rn.get(i++).toDecimal(); 517 //System.out.println("di = " + di); 518 //System.out.println("dd = " + dd); 519 assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0); 520 } 521 } catch (NoConvergenceException e) { 522 fail(e.toString()); 523 } 524 } 525 526 527 /** 528 * Test real root isolation and decimal refinement of Wilkinson polynomials, 529 * all roots. p = (x-0)*(x-1)*(x-2)*(x-3)*...*(x-n) 530 */ 531 public void testRealRootIsolationDecimalWilkinsonAll() { 532 final int N = 10; 533 d = dfac.getONE(); 534 e = dfac.univariate(0); 535 536 List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N); 537 a = d; 538 for (int i = 0; i < N; i++) { 539 c = dfac.fromInteger(i); 540 Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient())); 541 b = e.subtract(c); 542 a = a.multiply(b); 543 } 544 //System.out.println("a = " + a); 545 546 RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>(); 547 548 eps = eps.multiply(new BigRational(10000)); 549 //System.out.println("eps = " + eps); 550 BigDecimal eps1 = new BigDecimal(eps); 551 BigDecimal eps2 = eps1.multiply(new BigDecimal("100")); 552 //System.out.println("eps1 = " + eps1); 553 //System.out.println("eps2 = " + eps2); 554 555 List<BigDecimal> R = null; 556 R = rr.approximateRoots(a, eps); 557 //System.out.println("R = " + R); 558 assertTrue("#roots = " + N + " ", R.size() == N); 559 560 int i = 0; 561 for (BigDecimal dd : R) { 562 //System.out.println("dd = " + dd); 563 BigDecimal di = Rn.get(i++).toDecimal(); 564 //System.out.println("di = " + di); 565 assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0); 566 } 567 boolean t = rr.isApproximateRoot(R, a, eps); 568 assertTrue("some |a(dd)| < eps ", t); 569 } 570 571 572 /** 573 * Test real root isolation and decimal refinement of Wilkinson polynomials, 574 * inverse roots, all roots. p = (x-1)*(x-1/2)*(x-1/3)*...*(x-1/n) 575 */ 576 public void testRealRootIsolationDecimalWilkinsonInverseAll() { 577 final int N = 10; 578 d = dfac.getONE(); 579 e = dfac.univariate(0); 580 581 List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N); 582 a = d; 583 for (int i = 1; i < N; i++) { // use only for i > 0, since reverse 584 c = dfac.fromInteger(i); 585 if (i != 0) { 586 c = d.divide(c); 587 } 588 Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient())); 589 b = e.subtract(c); 590 a = a.multiply(b); 591 } 592 //System.out.println("a = " + a); 593 //System.out.println("Rn = " + Rn); 594 Collections.reverse(Rn); 595 //System.out.println("Rn = " + Rn); 596 597 RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>(); 598 599 eps = eps.multiply(new BigRational(1000000)); 600 //System.out.println("eps = " + eps); 601 BigDecimal eps1 = new BigDecimal(eps); 602 BigDecimal eps2 = eps1.multiply(new BigDecimal("10")); 603 //System.out.println("eps1 = " + eps1); 604 //System.out.println("eps2 = " + eps2); 605 606 List<BigDecimal> R = null; 607 R = rr.approximateRoots(a, eps); 608 //System.out.println("R = " + R); 609 assertTrue("#roots = " + (N - 1) + " ", R.size() == (N - 1)); 610 611 int i = 0; 612 for (BigDecimal dd : R) { 613 //System.out.println("dd = " + dd); 614 BigDecimal di = Rn.get(i++).toDecimal(); 615 //System.out.println("di = " + di); 616 assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0); 617 } 618 boolean t = rr.isApproximateRoot(R, a, eps); 619 assertTrue("some |a(dd)| < eps ", t); 620 } 621 622}