By the way, I forgot to mention that I think it would be better to just support negative exponents directly in Poly rather than creating a LaurentPoly class. At the very least, Poly shouldn't pretend it is doing Laurent polynomials (like with Poly(x + 1/x)) when it really isn't.
Aaron Meurer On Jun 16, 2011, at 6:19 PM, Aaron S. Meurer wrote: > Well, the thing that bugs me the most is that you can do the same things with > what would appear to be rational functions: > > In [7]: Poly(1/x) > Out[7]: Poly(1/x, 1/x, domain='ZZ') > > In [8]: Poly(1/x)*x > Out[8]: Poly(x*1/x, 1/x, domain='ZZ[x]') > > In [9]: Poly(1/x)*x - 1 > Out[9]: Poly(x*1/x - 1, 1/x, domain='ZZ[x]') > > In [10]: (Poly(1/x)*x - 1).is_zero > Out[10]: False > > This was the original problem in issue 2032. > > On Jun 16, 2011, at 2:51 PM, Mateusz Paprocki wrote: > >> Hi, >> >> On 16 June 2011 13:28, Aaron Meurer <[email protected]> wrote: >> This is related to >> http://code.google.com/p/sympy/issues/detail?id=2032. The polys >> pretend that they can work in K[x, 1/x], but they actually do not >> implement things properly, which can lead to wrong results: >> >> In [3]: Poly(exp(-x)) >> Out[3]: Poly(exp(-x), exp(-x), domain='ZZ') >> >> In [4]: Poly(exp(-x))*exp(x) >> Out[4]: Poly(exp(x)*exp(-x), exp(-x), domain='ZZ[exp(x)]') >> >> In [5]: Poly(exp(-x))*exp(x) - 1 >> Out[5]: Poly(exp(x)*exp(-x) - 1, exp(-x), domain='ZZ[exp(x)]') >> >> In [6]: (Poly(exp(-x))*exp(x) - 1).is_zero >> Out[6]: False >> >> And I could go further, as an incorrect is_zero result can easily be >> exacerbated to larger wrong results, but you get the idea. >> >> Result [6] is actually a valid result, because generators of polynomial >> expressions are assumed to be algebraically independent, so in [6] you can >> write (Poly(u)*v - 1).is_zero (which of course is False) as well. If you >> want [6] to give True, then you would have to imply an algebraic relation >> between exp(-x) and exp(x) (really u and v). A reasonable solution I see to >> this problem is to implement LaurentPoly which would allow negative >> exponents, which in turn would mean that exp(-x) and exp(x) would result in >> the same generator, with exponent -1 and 1 respectively. >> Another solution would be to add postprocessor to Poly, which would apply >> known algebraic relations of generators, after performing the main >> computation. > > This is a great idea. You see, there is actually a complete algorithm to > find algebraic relations on elementary functions, called the Risch structure > theorems. > > The transcendental Risch algorithm requires that each elementary extension be > transcendental over the others, or else certain results (namely nonelementary > results) will not be correct. The algorithm not only determines if a function > is transcendental over a differential extension, but provides the algebraic > relation when it is not. > > So it was necessary to implement this algorithm. So far, only the algorithm > for exponentials and logarithms has been implemented, but there is also an > algorithm for tangents, arctangents, and nonelementary extensions (functions > defined in terms of a nonelementary integral of a transcendental function, > like erf()). You can handle cosine and sine by first converting them to > tangents. > > If you want to play around with this, the easiest way is to import > DifferentialExtension from sympy.integrals.risch in my integration3 branch > and create a differential extension out of the expression. For example, you > can see it handles the above example with ease > > In [1]: from sympy.integrals.risch import DifferentialExtension > > In [2]: DifferentialExtension(exp(x) + exp(-x), x) > Out[2]: (Poly(_t0**2 + 1, _t0, domain='ZZ'), Poly(_t0, _t0, domain='ZZ'), > [Poly(1, x, domain='ZZ'), Poly(_t0, _t0, domain='ZZ')], [x, _t0], [Lambda(_i, > exp(_i))], [], [1], [x], [], []) > > The output consists of a lot of stuff, but the main thing is that the first > two Polys are the numerator and denominator of the resulting expression, and > the list of Lambdas near the end gives the function definition of t0, t1, t2, > etc. So in this case, we see that it has rewritten exp(x) + exp(-x) as > (exp(x)**2 + 1)/(exp(x)). Note that I have implemented heuristics to take > the exponentials in such a way to avoid algebraic relations, like > > In [3]: DifferentialExtension(exp(x) + exp(x/2) + exp(2*x), x) > Out[3]: (Poly(_t0**4 + _t0**2 + _t0, _t0, domain='ZZ'), Poly(1, _t0, > domain='ZZ'), [Poly(1, x, domain='ZZ'), Poly(1/2*_t0, _t0, domain='QQ')], [x, > _t0], [Lambda(_i, exp(_i/2))], [], [1], [x/2], [], []) > > Notice that t0 is exp(x/2). If we had made t0 exp(x), then exp(x/2) would > have been sqrt(t0), which is no good for the Risch algorithm (and also in the > polys), because we can't handle algebraic relations yet. Note that the > algorithms here are very powerful, and can handle much more nontrivial cases > than any heuristic would > > In [4]: DifferentialExtension(log(exp(x) + 1) + log(exp(-x) + 1), x) > Out[4]: > (Poly(2*_t1 - x, _t1, domain='ZZ[x]'), Poly(1, _t1, domain='ZZ'), [Poly(1, x, > domain='ZZ'), Poly(_t0, _t0, domain='ZZ'), Poly(_t0/(1 + _t0), _t1, > domain='ZZ(_t0)')], [x, _t0, _t1], [Lambda(_i, exp(_i) > ), Lambda(_i, log(1 + _t0))], [], [1], [x], [2], [1 + _t0]) > > (rewrote log(exp(x) + 1) + log(exp(-x) + 1) as 2*log(exp(x) + 1) - x. The > nature of the algorithms is to always rewrite the expression in terms of some > part of the rest. Sometimes, it is necessary to rewrite things in a more > canonical form before integrating, for example > > In [5]: DifferentialExtension(2**x, x) > Out[5]: > (Poly(_t0, _t0, domain='ZZ'), Poly(1, _t0, domain='ZZ'), [Poly(1, x, > domain='ZZ'), Poly(log(2)*_t0, _t0, domain='EX')], [x, _t0], [Lambda(_i, > exp(_i*log(2)))], [(exp(x*log(2)), 2**x)], [1], [x*log(2)] > , [], []) > > 2**x is rewritten as exp(x*log(2)). But the object contains an attribute > backsubs which puts it back to the way it was after integration > > In [6]: DifferentialExtension(2**x, x).backsubs > Out[6]: > ⎡⎛ x⋅log(2) x⎞⎤ > ⎣⎝ℯ , 2 ⎠⎦ > > Finally, I want to note that the algorithm for sin, cos, and tan is very > nontrivial. It's described in Bronstein's paper, "Simplification of Real > Elementary Functions," and involves rewriting all sin's and cos's as tangents > (using the half angle formulas), then the dependence on the tangents is given > by Flory's theorem (rewriting tangents in terms of symmetric polynomials). > > Right now, this code is all written with integration in mind, but it would > not be too difficult to adapt it to other uses, like finding algebraic > dependancies in Poly). > > By the way, this all applies only to the univariate case. I don't know if > any of it can be extended to the multivariate case. > > Aaron Meurer > > >> >> >> Aaron Meurer >> >> On Thu, Jun 16, 2011 at 12:15 PM, Mateusz Paprocki <[email protected]> wrote: >> > Hi, >> > >> > On 16 June 2011 20:09, smichr <[email protected]> wrote: >> >> >> >> The exp(x)*exp(-x) term in the Poly should cancel, shouldn't it? >> >> >>> Poly(exp(x) + exp(-x) - y)*exp(x) >> >> Poly(-y*exp(x) + exp(-x)*exp(x) + exp(x)**2, y, exp(-x), exp(x), >> >> domain='ZZ') >> > >> > Polynomials can't contain negative exponents, so exp(x) and exp(-x) = >> > 1/exp(x) are treated as two different generators. >> > >> >> >> >> -- >> >> You received this message because you are subscribed to the Google Groups >> >> "sympy" group. >> >> To post to this group, send email to [email protected]. >> >> To unsubscribe from this group, send email to >> >> [email protected]. >> >> For more options, visit this group at >> >> http://groups.google.com/group/sympy?hl=en. >> >> >> > >> > Mateusz >> > >> > -- >> > You received this message because you are subscribed to the Google Groups >> > "sympy" group. >> > To post to this group, send email to [email protected]. >> > To unsubscribe from this group, send email to >> > [email protected]. >> > For more options, visit this group at >> > http://groups.google.com/group/sympy?hl=en. >> > >> >> -- >> You received this message because you are subscribed to the Google Groups >> "sympy" group. >> To post to this group, send email to [email protected]. >> To unsubscribe from this group, send email to >> [email protected]. >> For more options, visit this group at >> http://groups.google.com/group/sympy?hl=en. >> >> >> Mateusz >> >> -- >> You received this message because you are subscribed to the Google Groups >> "sympy" group. >> To post to this group, send email to [email protected]. >> To unsubscribe from this group, send email to >> [email protected]. >> For more options, visit this group at >> http://groups.google.com/group/sympy?hl=en. > -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
