Wow, that's a big expression!

Unfortunately, as you've discovered, expand() is not exactly the
fastest part of SymPy.  The real solution here would be to implement a
better trigsimp algorithm.  The way you're doing it will result in
smaller expressions, but it's inefficient, as each replacement will
replace single terms with the sum of terms, requiring expansion to
simplify.

Disabling caching might help with the memory, but it could make things
go so slow that it's not worth it (especially expand, which uses it
heavily).

Looking at your expressions, I think the main simplifications that are
happening are these:

a*sin(x)**2 + a*cos(x)**2 => a
a*sin(x)*cos(x) => a/2*sin(2*x)
a*cos(x)**2 - a*sin(x)**2 => a*cos(2*x)

Would you agree?

Note that I've written these in a form that makes them useful as
transformations.

I think you may want to try writing a custom simplification routine
that uses these.  It would go something like this:

First, expand the original expression and get a list of the monomails.
 You can use poly() to do this.  For example:

In [100]: a = expressions.smallExpr

If you don't know before hand which trig functions are in your
expression, you can use something like this:

In [101]: from sympy.functions.elementary.trigonometric import
TrigonometricFunction

In [102]: a.atoms(TrigonometricFunction)
Out[102]: set([cos(q4), sin(q1), sin(q4), cos(q1), cos(q2), sin(q2),
sin(q3), cos(q3)])

In [116]: p = Poly(a, sin(q1), sin(q2), sin(q3), sin(q4), cos(q1),
cos(q2), cos(q3), cos(q4), sin(2*q1), sin(2*q2), sin(2*q3), sin(2*q4),
cos(2*q1), cos(2*q2), cos(2*q3), cos(2*q4))

In [117]: p.monoms()
Out[117]:
[(4, 4, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 3, 2, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 3, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 2, 3, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 2, 2, 1, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 2, 2, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 2, 1, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 2, 1, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0),
 (4, 2, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),...

Note that I only send the generators to poly() that are needed.  dq1,
etc. are treated as constants.  I include the double argument terms
because these will be needed after the transformations.

This list means that you have a terms with the given monomials.  For
example, (4, 4, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0) means that
there is a term sin(q1)**4*sin(q2)**4*sin(q3)*cos(q4).  You can obtain
the coefficient of this term by using p.coeffs() (they are in the same
order as the monomials, you can put them together with zip).

The idea would be to work on the monomials, which is much more space
efficient than working on the expressions directly.  The above rules
would then be something like (sorry, I've switched from numbering
starting at 1 to numbering starting at 0):

== a*sin(x)**2 + a*cos(x)**2 => a ==

- If (2 + s0, s1, s2, s3, c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0,
dc1, dc2, dc3) and (s0, s1, s2, s3, 2 + c0, c1, c2, c3, ds0, ds1, ds2,
ds3, dc0, dc1, dc2, dc3) are both monomails with the same coefficient
a, then they can be replaced with a single monomial (s0, s1, s2, s3,
c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2, dc3) with
coefficient a.

- If (s0, 2 + s1, s2, s3, c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0,
dc1, dc2, dc3) and (s0, s1, s2, s3, c0, 2 + c1, c2, c3, ds0, ds1, ds2,
ds3, dc0, dc1, dc2, dc3) are both monomails with the same coefficient
a, then they can be replaced with a single monomial (s0, s1, s2, s3,
c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2, dc3) with
coefficient a.

...

These rules can probably be generalized in to a single rule using some
more advanced logic.  Note how I allow for things like
sin(q1)**4*cos(q2)**2 + cos(q1)**2*sin(q1)**2*cos(q2)**2 =>
sin(q1)**2*cos(q2)**2.  2 + s0 would be implemented as "if s0 >= 2:".

== a*sin(x)*cos(x) => a/2*sin(2*x) ==

- If (n, s1, s2, s3, n, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2,
dc3) is a monomial with coefficient a, it can be replaced with (0, s1,
s2, s3, 0, c1, c2, c3, ds0 + n, ds1, ds2, ds3, dc0, dc1, dc2, dc3)
with coefficient a/2**n.

- If (s0, n, s2, s3, c0, n, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2,
dc3) is a monomial with coefficient a, it can be replaced with (0, s1,
s2, s3, 0, c1, c2, c3, ds0, ds1 + n, ds2, ds3, dc0, dc1, dc2, dc3)
with coefficient a/2**n.

...

Note that you have to use a/2**n (not a/2!).

== a*cos(x)**2 - sin(x)**2 => a*cos(2*x) ==

- If (s0, s1, s2, s3, 2 + c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0,
dc1, dc2, dc3) is a monomial with coefficient a (2 + s0, s1, s2, s3,
c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2, dc3) is a monomial
with coefficient -a, then they can be replaced with a single monomial
(s0, s1, s2, s3, c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0 + 1, dc1,
dc2, dc3) with coefficient a.

...

The algorithm would then be to iterate these rules until no changes
can be made.  The resulting expression is the fully simplified one.
You can probably be smart about the order that the rules are applied
in order to minimize the number of passes.

Be aware that I may have made typos or small logical errors above, but
you get the idea.  This algorithm could probably be generalized and
put into trigsimp(), but even a slightly more custom version that fits
your needs should not be too hard to write, and would be much faster
(the only thing that will bog it down memory wise will be the
expansion of the original expression).

It also could be generalized to handle 2*n angle arguments by
repeating the double angle rules over the functions that already have
double angles, if you find that is necessary.  And as you can imagine,
it can also be generalized in a bunch of other directions as well.

This logic could probably also be applied by using some fancy replace
calls, but the monomial version will be more precise, less bug prone
(because replace logic doesn't always work right for things like x**6
=> x**2*x**4), will be faster, and will use much less memory.

Sorry if this is more of a description of a solution rather than a
solution itself.  I don't really feel like implementing this right
now, though perhaps someone else will (potential GSoC students: this
would be a great way to fulfill your patch requirement).

Aaron Meurer

On Thu, Feb 16, 2012 at 1:54 PM, Vinzenz Bargsten <[email protected]> wrote:
> Hi,
>
> I'd like to symbolically simplify trigonometric expressions, which may
> become quite large.
> In another thread Aaron Meurer pointed out that
> expand(expand(a.rewrite(exp)).rewrite(cos)) may be a solution.
>
> The result is good for the small example expression (takes reasonable time).
> But for larger expressions the expand after rewrite(exp) won't finish or
> take too much memory. I tried to disable caching, to split the expression
> with e.g. Add.make_args and to disable some of expand's hints (also
> force=True or using expand_mul). That was rather not successful yet.
>
> Here is example code: http://www.ovgu.de/bargsten/download/simptest.tar.bz2
> (8KB)
>
> What would you recommend?
>
> Note that I have problems to run this code with my python3 sympy
> installation (immediately recursion limit exceeded). This is not a problem
> with python 2.7 .
>
> Regards,
> Vinzenz
>
> --
> 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.

Reply via email to