001/* 002 * $Id: ElementaryIntegration.java 4965 2014-10-17 20:07:51Z kredel $ 003 */ 004 005package edu.jas.integrate; 006 007 008import java.util.ArrayList; 009import java.util.List; 010import java.util.SortedMap; 011 012import org.apache.log4j.Logger; 013 014import edu.jas.poly.AlgebraicNumber; 015import edu.jas.poly.AlgebraicNumberRing; 016import edu.jas.poly.GenPolynomial; 017import edu.jas.poly.GenPolynomialRing; 018import edu.jas.poly.PolyUtil; 019import edu.jas.structure.GcdRingElem; 020import edu.jas.structure.Power; 021import edu.jas.structure.RingFactory; 022import edu.jas.ufd.FactorAbstract; 023import edu.jas.ufd.FactorFactory; 024import edu.jas.ufd.GCDFactory; 025import edu.jas.ufd.GreatestCommonDivisorAbstract; 026import edu.jas.ufd.GreatestCommonDivisorSubres; 027import edu.jas.ufd.PolyUfdUtil; 028import edu.jas.ufd.Quotient; 029import edu.jas.ufd.QuotientRing; 030import edu.jas.ufd.SquarefreeAbstract; 031import edu.jas.ufd.SquarefreeFactory; 032 033 034/** 035 * Methods related to elementary integration. In particular there are methods 036 * for Hermite reduction and Rothstein-Trager integration of the logarithmic 037 * part. 038 * 039 * @author Axel Kramer 040 * @author Heinz Kredel 041 * @param <C> coefficient type 042 */ 043 044public class ElementaryIntegration<C extends GcdRingElem<C>> { 045 046 047 private static final Logger logger = Logger.getLogger(ElementaryIntegration.class); 048 049 050 private final boolean debug = logger.isDebugEnabled(); 051 052 053 /** 054 * Engine for factorization. 055 */ 056 public final FactorAbstract<C> irr; 057 058 059 /** 060 * Engine for squarefree decomposition. 061 */ 062 public final SquarefreeAbstract<C> sqf; 063 064 065 /** 066 * Engine for greatest common divisors. 067 */ 068 public final GreatestCommonDivisorAbstract<C> ufd; 069 070 071 /** 072 * Constructor. 073 */ 074 public ElementaryIntegration(RingFactory<C> br) { 075 ufd = GCDFactory.<C> getProxy(br); 076 sqf = SquarefreeFactory.<C> getImplementation(br); 077 irr = /*(FactorAbsolute<C>)*/FactorFactory.<C> getImplementation(br); 078 } 079 080 081 /** 082 * Integration of a rational function. 083 * 084 * @param r rational function 085 * @return Integral container, such that integrate(r) = sum_i(g_i) + sum_j( 086 * an_j log(hd_j) ) 087 */ 088 public QuotIntegral<C> integrate(Quotient<C> r) { 089 Integral<C> integral = integrate(r.num, r.den); 090 return new QuotIntegral<C>(r.ring, integral); 091 } 092 093 094 /** 095 * Integration of a rational function. 096 * 097 * @param a numerator 098 * @param d denominator 099 * @return Integral container, such that integrate(a/d) = sum_i(gn_i/gd_i) + 100 * integrate(h0) + sum_j( an_j log(hd_j) ) 101 */ 102 public Integral<C> integrate(GenPolynomial<C> a, GenPolynomial<C> d) { 103 if (d == null || a == null || d.isZERO()) { 104 throw new IllegalArgumentException("zero or null not allowed"); 105 } 106 if (a.isZERO()) { 107 return new Integral<C>(a, d, a); 108 } 109 if (d.isONE()) { 110 GenPolynomial<C> pi = PolyUtil.<C> baseIntegral(a); 111 return new Integral<C>(a, d, pi); 112 } 113 GenPolynomialRing<C> pfac = d.ring; 114 if (pfac.nvar > 1) { 115 throw new IllegalArgumentException("only for univariate polynomials " + pfac); 116 } 117 if (!pfac.coFac.isField()) { 118 throw new IllegalArgumentException("only for field coefficients " + pfac); 119 } 120 121 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(a, d); 122 GenPolynomial<C> p = qr[0]; 123 GenPolynomial<C> r = qr[1]; 124 125 GenPolynomial<C> c = ufd.gcd(r, d); 126 if (!c.isONE()) { 127 r = PolyUtil.<C> basePseudoQuotientRemainder(r, c)[0]; 128 d = PolyUtil.<C> basePseudoQuotientRemainder(d, c)[0]; 129 } 130 List<GenPolynomial<C>>[] ih = integrateHermite(r, d); 131 List<GenPolynomial<C>> rat = ih[0]; 132 List<GenPolynomial<C>> log = ih[1]; 133 134 GenPolynomial<C> pp = log.remove(0); 135 p = p.sum(pp); 136 GenPolynomial<C> pi = PolyUtil.<C> baseIntegral(p); 137 138 if (debug) { 139 logger.debug("pi = " + pi); 140 logger.debug("rat = " + rat); 141 logger.debug("log = " + log); 142 } 143 if (log.size() == 0) { 144 return new Integral<C>(a, d, pi, rat); 145 } 146 147 List<LogIntegral<C>> logi = new ArrayList<LogIntegral<C>>(log.size() / 2); 148 for (int i = 0; i < log.size(); i++) { 149 GenPolynomial<C> ln = log.get(i++); 150 GenPolynomial<C> ld = log.get(i); 151 LogIntegral<C> pf = integrateLogPart(ln, ld); 152 logi.add(pf); 153 } 154 if (debug) { 155 logger.debug("logi = " + logi); 156 } 157 return new Integral<C>(a, d, pi, rat, logi); 158 } 159 160 161 /** 162 * Integration of the rational part, Hermite reduction step. 163 * 164 * @param a numerator 165 * @param d denominator, gcd(a,d) == 1 166 * @return [ [ gn_i, gd_i ], [ h0, hn_j, hd_j ] ] such that integrate(a/d) = 167 * sum_i(gn_i/gd_i) + integrate(h0) + sum_j( integrate(hn_j/hd_j) ) 168 */ 169 @SuppressWarnings("cast") 170 public List<GenPolynomial<C>>[] integrateHermite(GenPolynomial<C> a, GenPolynomial<C> d) { 171 if (d == null || d.isZERO()) { 172 throw new IllegalArgumentException("d == null or d == 0"); 173 } 174 if (a == null || a.isZERO()) { 175 throw new IllegalArgumentException("a == null or a == 0"); 176 } 177 178 // get squarefree decomposition 179 SortedMap<GenPolynomial<C>, Long> sfactors = sqf.squarefreeFactors(d); 180 181 List<GenPolynomial<C>> D = new ArrayList<GenPolynomial<C>>(sfactors.keySet()); 182 List<GenPolynomial<C>> DP = new ArrayList<GenPolynomial<C>>(); 183 for (GenPolynomial<C> f : D) { 184 long e = sfactors.get(f); 185 GenPolynomial<C> dp = Power.<GenPolynomial<C>> positivePower(f, e); 186 DP.add(dp); 187 } 188 //System.out.println("D: " + D); 189 //System.out.println("DP: " + DP); 190 191 // get partial fraction decompostion 192 List<GenPolynomial<C>> Ai = ufd.basePartialFraction(a, DP); 193 //System.out.println("Ai: " + Ai); 194 195 List<GenPolynomial<C>> G = new ArrayList<GenPolynomial<C>>(); 196 List<GenPolynomial<C>> H = new ArrayList<GenPolynomial<C>>(); 197 H.add(Ai.remove(0)); // P 198 199 GenPolynomialRing<C> fac = d.ring; 200 int i = 0; 201 for (GenPolynomial<C> v : D) { 202 //System.out.println("V:" + v.toString()); 203 GenPolynomial<C> Ak = Ai.get(i++); 204 //System.out.println("Ak: " + Ak.toString()); 205 int k = sfactors.get(v).intValue(); // assert low power 206 for (int j = k - 1; j >= 1; j--) { 207 //System.out.println("Step(" + k + "," + j + ")"); 208 GenPolynomial<C> DV_dx = PolyUtil.<C> baseDeriviative(v); 209 GenPolynomial<C> Aik = Ak.divide(fac.fromInteger(-j)); 210 GenPolynomial<C>[] BC = ufd.baseGcdDiophant(DV_dx, v, Aik); 211 GenPolynomial<C> b = BC[0]; 212 GenPolynomial<C> c = BC[1]; 213 GenPolynomial<C> vj = Power.<GenPolynomial<C>> positivePower(v, j); 214 G.add(b); // B 215 G.add(vj); // v^j 216 Ak = fac.fromInteger(-j).multiply(c).subtract(PolyUtil.<C> baseDeriviative(b)); 217 //System.out.println("B: " + b.toString()); 218 //System.out.println("C: " + c.toString()); 219 } 220 //System.out.println("V:" + v.toString()); 221 //System.out.println("Ak: " + Ak.toString()); 222 if (!Ak.isZERO()) { 223 H.add(Ak); // A_k 224 H.add(v); // v 225 } 226 } 227 List<GenPolynomial<C>>[] ret = (List<GenPolynomial<C>>[]) new List[2]; 228 ret[0] = G; 229 ret[1] = H; 230 return ret; 231 } 232 233 234 /** 235 * Univariate GenPolynomial integration of the logaritmic part, 236 * Rothstein-Trager algorithm. 237 * @param A univariate GenPolynomial, deg(A) < deg(P). 238 * @param P univariate squarefree GenPolynomial, gcd(A,P) == 1. 239 * @return logarithmic part container. 240 */ 241 public LogIntegral<C> integrateLogPart(GenPolynomial<C> A, GenPolynomial<C> P) { 242 if (P == null || P.isZERO()) { 243 throw new IllegalArgumentException(" P == null or P == 0"); 244 } 245 if (A == null || A.isZERO()) { 246 throw new IllegalArgumentException(" A == null or A == 0"); 247 } 248 //System.out.println("\nP_base_algeb_part = " + P); 249 GenPolynomialRing<C> pfac = P.ring; // K[x] 250 if (pfac.nvar > 1) { 251 throw new IllegalArgumentException("only for univariate polynomials " + pfac); 252 } 253 if (!pfac.coFac.isField()) { 254 throw new IllegalArgumentException("only for field coefficients " + pfac); 255 } 256 List<C> cfactors = new ArrayList<C>(); 257 List<GenPolynomial<C>> cdenom = new ArrayList<GenPolynomial<C>>(); 258 List<AlgebraicNumber<C>> afactors = new ArrayList<AlgebraicNumber<C>>(); 259 List<GenPolynomial<AlgebraicNumber<C>>> adenom = new ArrayList<GenPolynomial<AlgebraicNumber<C>>>(); 260 261 // P linear 262 if (P.degree(0) <= 1) { 263 cfactors.add(A.leadingBaseCoefficient()); 264 cdenom.add(P); 265 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom); 266 } 267 List<GenPolynomial<C>> Pfac = irr.baseFactorsSquarefree(P); 268 //System.out.println("\nPfac = " + Pfac); 269 270 List<GenPolynomial<C>> Afac = ufd.basePartialFraction(A, Pfac); 271 272 GenPolynomial<C> A0 = Afac.remove(0); 273 if (!A0.isZERO()) { 274 throw new RuntimeException(" A0 != 0: deg(A)>= deg(P)"); 275 } 276 277 // algebraic and linear factors 278 int i = 0; 279 for (GenPolynomial<C> pi : Pfac) { 280 GenPolynomial<C> ai = Afac.get(i++); 281 if (pi.degree(0) <= 1) { 282 cfactors.add(ai.leadingBaseCoefficient()); 283 cdenom.add(pi); 284 continue; 285 } 286 LogIntegral<C> pf = integrateLogPartIrreducible(ai, pi); 287 cfactors.addAll(pf.cfactors); 288 cdenom.addAll(pf.cdenom); 289 afactors.addAll(pf.afactors); 290 adenom.addAll(pf.adenom); 291 } 292 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom); 293 } 294 295 296 /** 297 * Univariate GenPolynomial integration of the logaritmic part, 298 * Rothstein-Trager algorithm. 299 * @param A univariate GenPolynomial, deg(A) < deg(P). 300 * @param P univariate irreducible GenPolynomial. // gcd(A,P) == 1 automatic 301 * @return logarithmic part container. 302 */ 303 public LogIntegral<C> integrateLogPartIrreducible(GenPolynomial<C> A, GenPolynomial<C> P) { 304 if (P == null || P.isZERO()) { 305 throw new IllegalArgumentException("P == null or P == 0"); 306 } 307 //System.out.println("\nP_base_algeb_part = " + P); 308 GenPolynomialRing<C> pfac = P.ring; // K[x] 309 if (pfac.nvar > 1) { 310 throw new IllegalArgumentException("only for univariate polynomials " + pfac); 311 } 312 if (!pfac.coFac.isField()) { 313 throw new IllegalArgumentException("only for field coefficients " + pfac); 314 } 315 List<C> cfactors = new ArrayList<C>(); 316 List<GenPolynomial<C>> cdenom = new ArrayList<GenPolynomial<C>>(); 317 List<AlgebraicNumber<C>> afactors = new ArrayList<AlgebraicNumber<C>>(); 318 List<GenPolynomial<AlgebraicNumber<C>>> adenom = new ArrayList<GenPolynomial<AlgebraicNumber<C>>>(); 319 320 // P linear 321 if (P.degree(0) <= 1) { 322 cfactors.add(A.leadingBaseCoefficient()); 323 cdenom.add(P); 324 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom); 325 } 326 327 // deriviative 328 GenPolynomial<C> Pp = PolyUtil.<C> baseDeriviative(P); 329 //no: Pp = Pp.monic(); 330 //System.out.println("\nP = " + P); 331 //System.out.println("Pp = " + Pp); 332 333 // Q[t] 334 String[] vars = new String[] { "t" }; 335 GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(pfac.coFac, 1, pfac.tord, vars); 336 GenPolynomial<C> t = cfac.univariate(0); 337 //System.out.println("t = " + t); 338 339 // Q[x][t] 340 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(pfac, cfac); // sic 341 //System.out.println("rfac = " + rfac.toScript()); 342 343 // transform polynomials to bi-variate polynomial 344 GenPolynomial<GenPolynomial<C>> Ac = PolyUfdUtil.<C> introduceLowerVariable(rfac, A); 345 //System.out.println("Ac = " + Ac); 346 GenPolynomial<GenPolynomial<C>> Pc = PolyUfdUtil.<C> introduceLowerVariable(rfac, P); 347 //System.out.println("Pc = " + Pc); 348 GenPolynomial<GenPolynomial<C>> Pcp = PolyUfdUtil.<C> introduceLowerVariable(rfac, Pp); 349 //System.out.println("Pcp = " + Pcp); 350 351 // Q[t][x] 352 GenPolynomialRing<GenPolynomial<C>> rfac1 = Pc.ring; 353 //System.out.println("rfac1 = " + rfac1.toScript()); 354 355 // A - t P' 356 GenPolynomial<GenPolynomial<C>> tc = rfac1.getONE().multiply(t); 357 //System.out.println("tc = " + tc); 358 GenPolynomial<GenPolynomial<C>> At = Ac.subtract(tc.multiply(Pcp)); 359 //System.out.println("At = " + At); 360 361 GreatestCommonDivisorSubres<C> engine = new GreatestCommonDivisorSubres<C>(); 362 // = GCDFactory.<C>getImplementation( cfac.coFac ); 363 GreatestCommonDivisorAbstract<AlgebraicNumber<C>> aengine = null; 364 365 GenPolynomial<GenPolynomial<C>> Rc = engine.recursiveUnivariateResultant(Pc, At); 366 //System.out.println("Rc = " + Rc); 367 GenPolynomial<C> res = Rc.leadingBaseCoefficient(); 368 //no: res = res.monic(); 369 //System.out.println("\nres = " + res); 370 371 SortedMap<GenPolynomial<C>, Long> resfac = irr.baseFactors(res); 372 //System.out.println("resfac = " + resfac + "\n"); 373 374 for (GenPolynomial<C> r : resfac.keySet()) { 375 //System.out.println("\nr(t) = " + r); 376 if (r.isConstant()) { 377 continue; 378 } 379 //vars = new String[] { "z_" + Math.abs(r.hashCode() % 1000) }; 380 vars = pfac.newVars("z_"); 381 pfac = pfac.copy(); 382 @SuppressWarnings("unused") 383 String[] unused = pfac.setVars(vars); 384 r = pfac.copy(r); // hack to exchange the variables 385 //System.out.println("r(z_) = " + r); 386 AlgebraicNumberRing<C> afac = new AlgebraicNumberRing<C>(r, true); // since irreducible 387 logger.debug("afac = " + afac.toScript()); 388 AlgebraicNumber<C> a = afac.getGenerator(); 389 //no: a = a.negate(); 390 //System.out.println("a = " + a); 391 392 // K(alpha)[x] 393 GenPolynomialRing<AlgebraicNumber<C>> pafac = new GenPolynomialRing<AlgebraicNumber<C>>(afac, 394 Pc.ring); 395 //System.out.println("pafac = " + pafac.toScript()); 396 397 // convert to K(alpha)[x] 398 GenPolynomial<AlgebraicNumber<C>> Pa = PolyUtil.<C> convertToAlgebraicCoefficients(pafac, P); 399 //System.out.println("Pa = " + Pa); 400 GenPolynomial<AlgebraicNumber<C>> Pap = PolyUtil.<C> convertToAlgebraicCoefficients(pafac, Pp); 401 //System.out.println("Pap = " + Pap); 402 GenPolynomial<AlgebraicNumber<C>> Aa = PolyUtil.<C> convertToAlgebraicCoefficients(pafac, A); 403 //System.out.println("Aa = " + Aa); 404 405 // A - a P' 406 GenPolynomial<AlgebraicNumber<C>> Ap = Aa.subtract(Pap.multiply(a)); 407 //System.out.println("Ap = " + Ap); 408 409 if (aengine == null) { 410 aengine = GCDFactory.<AlgebraicNumber<C>> getImplementation(afac); 411 } 412 GenPolynomial<AlgebraicNumber<C>> Ga = aengine.baseGcd(Pa, Ap); 413 //System.out.println("Ga = " + Ga); 414 if (Ga.isConstant()) { 415 //System.out.println("warning constant gcd ignored"); 416 continue; 417 } 418 afactors.add(a); 419 adenom.add(Ga); 420 // special quadratic case 421 if (P.degree(0) == 2 && Ga.degree(0) == 1) { 422 GenPolynomial<AlgebraicNumber<C>>[] qra = PolyUtil 423 .<AlgebraicNumber<C>> basePseudoQuotientRemainder(Pa, Ga); 424 GenPolynomial<AlgebraicNumber<C>> Qa = qra[0]; 425 if (!qra[1].isZERO()) { 426 throw new ArithmeticException("remainder not zero"); 427 } 428 //System.out.println("Qa = " + Qa); 429 afactors.add(a.negate()); 430 adenom.add(Qa); 431 } 432 // todo: eventually implement special cases deg = 3, 4 433 } 434 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom); 435 } 436 437 438 /** 439 * Derivation of a univariate rational function. 440 * 441 * @param r rational function 442 * @return dr/dx 443 */ 444 public Quotient<C> deriviative(Quotient<C> r) { 445 GenPolynomial<C> num = r.num; 446 GenPolynomial<C> den = r.den; 447 GenPolynomial<C> nump = PolyUtil.<C> baseDeriviative(num); 448 if (den.isONE()) { 449 return new Quotient<C>(r.ring, nump, den); 450 } 451 GenPolynomial<C> denp = PolyUtil.<C> baseDeriviative(den); 452 453 GenPolynomial<C> n = den.multiply(nump).subtract(num.multiply(denp)); 454 GenPolynomial<C> d = den.multiply(den); 455 456 Quotient<C> der = new Quotient<C>(r.ring, n, d); 457 return der; 458 } 459 460 461 /** 462 * Test of integration of a rational function. 463 * 464 * @param ri integral 465 * @return true, if ri is an integral, else false. 466 */ 467 public boolean isIntegral(QuotIntegral<C> ri) { 468 Quotient<C> r = ri.quot; 469 QuotientRing<C> qr = r.ring; 470 Quotient<C> i = r.ring.getZERO(); 471 for (Quotient<C> q : ri.rational) { 472 Quotient<C> qd = deriviative(q); 473 i = i.sum(qd); 474 } 475 if (ri.logarithm.size() == 0) { 476 return r.equals(i); 477 } 478 for (LogIntegral<C> li : ri.logarithm) { 479 Quotient<C> q = new Quotient<C>(qr, li.num, li.den); 480 i = i.sum(q); 481 } 482 boolean t = r.equals(i); 483 if (!t) { 484 return false; 485 } 486 for (LogIntegral<C> li : ri.logarithm) { 487 t = isIntegral(li); 488 if (!t) { 489 return false; 490 } 491 } 492 return true; 493 } 494 495 496 /** 497 * Test of integration of the logarithmic part of a rational function. 498 * 499 * @param rl logarithmic part of an integral 500 * @return true, if rl is an integral, else false. 501 */ 502 public boolean isIntegral(LogIntegral<C> rl) { 503 QuotientRing<C> qr = new QuotientRing<C>(rl.den.ring); 504 Quotient<C> r = new Quotient<C>(qr, rl.num, rl.den); 505 506 Quotient<C> i = qr.getZERO(); 507 int j = 0; 508 for (GenPolynomial<C> d : rl.cdenom) { 509 GenPolynomial<C> dp = PolyUtil.<C> baseDeriviative(d); 510 dp = dp.multiply(rl.cfactors.get(j++)); 511 Quotient<C> f = new Quotient<C>(qr, dp, d); 512 i = i.sum(f); 513 } 514 if (rl.afactors.size() == 0) { 515 return r.equals(i); 516 } 517 r = r.subtract(i); 518 QuotientRing<AlgebraicNumber<C>> aqr = new QuotientRing<AlgebraicNumber<C>>(rl.adenom.get(0).ring); 519 Quotient<AlgebraicNumber<C>> ai = aqr.getZERO(); 520 521 GenPolynomial<AlgebraicNumber<C>> aqn = PolyUtil.<C> convertToAlgebraicCoefficients(aqr.ring, r.num); 522 GenPolynomial<AlgebraicNumber<C>> aqd = PolyUtil.<C> convertToAlgebraicCoefficients(aqr.ring, r.den); 523 Quotient<AlgebraicNumber<C>> ar = new Quotient<AlgebraicNumber<C>>(aqr, aqn, aqd); 524 525 j = 0; 526 for (GenPolynomial<AlgebraicNumber<C>> d : rl.adenom) { 527 GenPolynomial<AlgebraicNumber<C>> dp = PolyUtil.<AlgebraicNumber<C>> baseDeriviative(d); 528 dp = dp.multiply(rl.afactors.get(j++)); 529 Quotient<AlgebraicNumber<C>> f = new Quotient<AlgebraicNumber<C>>(aqr, dp, d); 530 ai = ai.sum(f); 531 } 532 boolean t = ar.equals(ai); 533 if (t) { 534 return true; 535 } 536 logger.warn("log integral not verified"); 537 //System.out.println("r = " + r); 538 //System.out.println("afactors = " + rl.afactors); 539 //System.out.println("adenom = " + rl.adenom); 540 //System.out.println("ar = " + ar); 541 //System.out.println("ai = " + ai); 542 return true; 543 } 544 545}