Hi all,
I managed to improve the situation quite a lot. The new code is this
time attached as a patch, or alternatively see my (ness01) branch trigsimp.
Executive summary: run mytrigsimp, test the function mytrigsimp. Ignore
the debugging output, pass order=grlex or grevlex for sensible results,
pass n=2, n=3 etc to improve the results.
The function mytrigsimp
-----------------------
As outlined in the previous mail, mytrigsimp can be used to simplify a
fraction involving cos and sin terms. As currently designed it
specifically simplifes cos(n*x) and sin(n*x) stuff. (You can pass pretty
much anything, but nothing will happen to other terms except for slowing
things down :D). There is currently quite a lot of debug output which
you can ignore. I also added some short docstrings. Specifically,
mytrigsimp is invoked as follows:
mytrigsimp(expr, **opts)
Opts is mostly passed through to polys. You probably shouldn't use the
gens=.. option. It is probably a good idea to pass order=grlex or
order=grevlex -- for reasons that are not clear to me the default order
in polys is lex. In any case, unless you pass a graded order (i.e. grlex
or grevlex), the result will *not* be of minimal total degree, in the
naive sense.
There is one non-polys option, which is currently called n. It can be
used to extend the search space. Namely, suppose you want so simplify
2*sin(x)*cos(x). We want the answer to be sin(2*x). But currently
mytrigsimp is rather conservative, and since it sees no sin(2*x) it does
not try to find an expression involving that. Passing n=2 will double
the search space, in this case it will try sin(2*x) and cos(2*x).
You will notice that performance degrades quickly as you increase either
the number of generators or n. Also, if n is "big" funny things can
happen (e.g. try n=4 for sin(x)*cos(x) :D). The reason is that
ratsimpmodprime returns a fraction of minimal degree, but does no other
selection.
Changes to ratsimpmodprime
--------------------------
In the last mail I reported a bug in ratsimpmodprime. That was really
more a bug in my code and a failure to properly signal a problem in
ratsimpmodprime. (Specifically, there is a certain situation in which
ratsimpmodprime divides by zero, but this can only happen if the ideal
is not prime, or if the generating set is not a groebner basis.) There
also *is* a bug in ratsimpmodprime, but it is fairly minor (results
correct up to a rational number as prefactor -- fixed in my branch).
I also implemented a (tiny) change which resulted in huge performance
improvements. I'm not 100% sure it is correct, so I'd be glad if someone
could verify the following:
Ratsimpmodprime simplifies f/g mod I, essentially by writing out a
general fraction p/q of degree n in the generators (say (a + b*x +
c*y)/(e + f*x + g*y) for n=1, generators x,y), and then checking if
there exists a choice of coefficients such that f*q - p*g is in I. What
the original code did was to change the ground field from QQ to QQ(a, b,
c, e, f, g) and then reduce f*q-p*g modulo the groebner basis (one may
easily see that changing the ground field in this way does indeed
preserve groebner bases). This will result in an expression linear in
the new coefficients, and finding a solution is a matter of linear
algebra. The problem is that computations over this field of fractions
are slow.
My new method adds a, b, c, e, f, g to the generators of our polynomial
ring, that is we now work in k[x, y, a, .., g] instead of k[x, y]. Again
groebner bases are preserved under this operaton, and I believe a (the)
"completely reduced normal form" has exactly the properties we want, as
outlined above.
Discussion and Future
---------------------
It seems to me that this method is quite powerful. I believe the most
important thing now that this proof of concept code exists is to
investigate what it can and cannot do. So please test it thoroughly and
report all surprises :-). I believe we can get a fairly good trigsimp
along the following lines:
1. Replace all trig expressions by sin and cos equivalents.
2. Run mytrigsimp, probably some clever improved version.
3. Run some clever heuristics to get tan and similar back.
4. Run the old trigsimp.
I see the following problems/limitations with the code, at least in its
current form:
1. It only handles sin and cos.
2. It quickly gets slow.
3. Answers can be "funny".
Let me look at these in the opposite order *g*.
Regarding (3), I think the solutions ratsimpmodprime can be more
complicated than necessary (i.e. of minimal degree, but it could be
possible using fewer terms).
For example, consider
expr = (4*sin(x)*cos(x) + 12*sin(x) + 5*cos(x)**3 + 21*cos(x)**2 +
23*cos(x) + 15)/(-sin(x)*cos(x)**2 + 2*sin(x)*cos(x) + 15*sin(x) +
7*cos(x)**3 + 31*cos(x)**2 + 37*cos(x) + 21)
Using n=1 you get something quite short. For n=2 you get something
longer (but of the same total degree). The reason here is that the
normal form the algorithm *starts* with depends on the ideal, and may be
nasty. But it is of degree (1,1), so the algorithm only tries (1, 0) and
(0, 1) - neither of which works. Even if we change it to also try (1,1),
the result returned is still not that nice. The reason is that in the
algorithm, when presented with choice, has currently no heuristics (I'll
try to change that).
[This expr is actually quite fun. Try simplifying various powers with
various choices of n ...]
Regarding (2), I don't see much to do about that. We could try to
precompute groebner bases, or at least cache them intelligently. This
would probably speed up common cases. As the exrpession gets larger (in
particular higher degree or more variables), the reduction computations
also quickly become slow.
Regarding (1), there are actually some possibilities. ratsimpmodprime is
quite general, so various simplification strategies are possible. The
problem is that we have to produce a prime ideal of polynomial
relations, which is not quite trivial...
Anyhow, I have spent much longer than I intended already on writing this
mail. Let me know what you think!
Best,
Tom
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sympy@googlegroups.com.
To unsubscribe from this group, send email to
sympy+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
diff --git a/mytrigsimp.py b/mytrigsimp.py
new file mode 100644
index 0000000..02dbae4
--- /dev/null
+++ b/mytrigsimp.py
@@ -0,0 +1,64 @@
+from sympy import *
+
+def build_ideal(x, n):
+ """
+ Build an ideal to "accomodate" sin(n*x). Return new generators.
+
+ Basically, this decides what expressions to consider in simplifying.
+ Currently, we just return sin(x), sin(2*x), .., sin(n*x), cos(x), ...,
cos(n*x).
+ The main tradeoff here is performance: the more expressions we introduce,
+ the slower the simplification process (but the better the possible
results).
+ """
+ I = [cos(x)**2 + sin(x)**2 - 1]
+ g = [cos(x), sin(x)]
+ y = Dummy('y')
+ for k in range(2, n + 1):
+ cn = cos(k*y).expand(trig=True).subs(y, x)
+ sn = sin(k*y).expand(trig=True).subs(y, x)
+ I.append(cos(k*x) - cn)
+ I.append(sin(k*x) - sn)
+ g.extend([cos(k*x), sin(k*x)])
+ return I, g
+
+def analyse_gens(gens, m):
+ """
+ Build an ideal to simplify generators ``gens``.
+
+ Return generators for the ideal, and new generators.
+ ``m`` is a hint: if sin(n*x) is encountered, try simplifying up to
+ sin(n*m*x).
+ """
+ trigterms = [g.args[0].as_coeff_mul() for g in gens if g.func in [sin,
cos]]
+ newgens = [g for g in gens if g.func not in [sin, cos]]
+ trigdict = {}
+ for val, key in trigterms:
+ trigdict.setdefault(key, []).append(val)
+ res = []
+ for key, val in trigdict.iteritems():
+ gcd = reduce(igcd, val)
+ n = max(x / gcd for x in val) * m
+ r, n = build_ideal(gcd*Mul(*key), n)
+ res.extend(r)
+ newgens.extend(n)
+ return res, newgens
+
+def mytrigsimp(expr, **opts):
+ """
+ Simplify the sine and cosine terms in expr.
+
+ Useful options:
+
+ - order=grlex or grevlex - probably faster
+ - n=2, 3, ... - increase search space
+ """
+ n = opts.pop('n', 1)
+ num, denom = cancel(expr).as_numer_denom()
+ polys, opt = parallel_poly_from_expr([num, denom], **opts)
+ print 'initial gens:', opt.gens
+ ideal, gens = analyse_gens(opt.gens, n)
+ print 'ideal:', ideal
+ print 'new gens:', gens, " -- len", len(gens)
+ opts['gens'] = gens
+ G = groebner(ideal, **opts)
+ print 'groebner basis:', list(G), " -- len", len(G)
+ return ratsimpmodprime(expr, list(G), **opts)
diff --git a/sympy/simplify/simplify.py b/sympy/simplify/simplify.py
index 336a098..477e75a 100644
--- a/sympy/simplify/simplify.py
+++ b/sympy/simplify/simplify.py
@@ -847,14 +847,19 @@ def _ratsimpmodprime(a, b, N=0, D=0):
Cs = symbols("c:%d" % len(M1))
Ds = symbols("d:%d" % len(M2))
+ ng = Cs + Ds
- c_hat = Poly(sum([Cs[i] * M1[i] for i in xrange(len(M1))]),
opt.gens)
- d_hat = Poly(sum([Ds[i] * M2[i] for i in xrange(len(M2))]),
opt.gens)
+ c_hat = Poly(sum([Cs[i] * M1[i] for i in xrange(len(M1))]),
opt.gens+ng)
+ d_hat = Poly(sum([Ds[i] * M2[i] for i in xrange(len(M2))]),
opt.gens+ng)
- r = reduced(a * d_hat - b * c_hat, G, opt.gens, order=opt.order,
polys=True)[1]
+ print "start reduce"
+ r = reduced(a * d_hat - b * c_hat, G, opt.gens+ng,
order=opt.order, polys=True)[1]
+ print "end"
- S = r.coeffs()
+ S = Poly(r, gens=opt.gens).coeffs()
+ print "solving:", len(S), len(Cs + Ds)
sol = solve(S, Cs + Ds)
+ print "solved"
# If nontrivial solutions exist, solve will give them
# parametrized, i.e. the values of some keys will be
@@ -863,7 +868,7 @@ def _ratsimpmodprime(a, b, N=0, D=0):
for key in sol.keys():
sol[key] = sol[key].subs(dict(zip(Cs + Ds, [1] * (len(Cs) +
len(Ds)))))
- if sol and not all([s == 0 for s in sol.itervalues()]):
+ if sol and not all(s == 0 for var, s in sol.iteritems()):
c = c_hat.subs(sol)
d = d_hat.subs(sol)
@@ -875,6 +880,8 @@ def _ratsimpmodprime(a, b, N=0, D=0):
c = Poly(c, opt.gens)
d = Poly(d, opt.gens)
+ if d == 0:
+ raise ValueError('Ideal not prime?')
break
@@ -896,10 +903,10 @@ def _ratsimpmodprime(a, b, N=0, D=0):
c, d = _ratsimpmodprime(Poly(num, opt.gens), Poly(denom, opt.gens))
if not domain.has_Field:
- c = c.clear_denoms(convert=True)[1]
- d = d.clear_denoms(convert=True)[1]
+ cn, c = c.clear_denoms(convert=True)
+ dn, d = d.clear_denoms(convert=True)
- return c/d
+ return c/d*(dn/cn)
def trigsimp(expr, deep=False, recursive=False):
"""
