On Apr 17, 2011, at 1:28 AM, Mateusz Paprocki wrote:
> Hi,
>
> On 17 April 2011 08:33, refp16 <[email protected]> wrote:
> Mateusz and Aaron,
>
> Thank you both for your responses.
>
> Precisely what I was looking for: the case where I don´t know the
> power of the polynomial.
>
> sorted(auxex.args, cmp=Basic.compare_pretty)[6]
>
> works perfectly, but,
>
> Poly(auxex).terms()[-7]
>
> gives me an error: Traceback (most recent call last):
> File "<editor selection>", line 2, in <module>
> AttributeError: 'Poly' object has no attribute 'terms'
>
> Why is this?
>
> This was added after SymPy 0.6.7 (and btw. it should be all_terms()).
Why all_terms()? That gives every term, even if the coefficient is zero. That
is not only a waste, but it makes it much harder to find the seventh term.
>
>
> Also, I managed to substitute
>
> product(1 + n*x**3**n, (n, 1, 3))
>
> for my original FOR loop, making the code shorter.
>
> This is short but for larger `n` not very efficient, e.g.:
>
> In [1]: %time f = product(1 + n*x**3**n, (n, 1, 2000))
> CPU times: user 1.72 s, sys: 0.00 s, total: 1.72 s
> Wall time: 1.78 s
>
> In [3]: %time g = Mul(*[ 1 + n*x**3**n for n in xrange(1, 2000+1) ])
> CPU times: user 0.15 s, sys: 0.00 s, total: 0.15 s
> Wall time: 0.15 s
>
> In [5]: f == g
> Out[5]: True
product should be fixed to work like Mul in this case. There's no reason why
it should be slower.
>
>
> One last question. I´m trying to expand that expression up to n=2000,
> but I get a memory error. Besides buying more physical RAM, can you
> recommend any other way of doing it?
>
> n=2000 is definitively too much for SymPy. I was able to compute this for
> n=21 in reasonable time (under a minute) using the following code:
>
> >>> from sympy.polys.monomialtools import *
> >>> from sympy.polys.groebnertools import *
>
> >>> O = monomial_lex_key
> >>> f = [((0,), ZZ(1))]
>
> >>> for n in xrange(1, 2000+1):
> ... print n, sdp_LT(f, 0, ZZ), len(f)
> ... f = sdp_mul(f, [((3**n,), ZZ(n)), ((0,), ZZ(1))], 0, O, ZZ)
>
> 1 ((0,), mpz(1)) 1
> 2 ((3,), mpz(1)) 2
> 3 ((12,), mpz(2)) 4
> 4 ((39,), mpz(6)) 8
> 5 ((120,), mpz(24)) 16
> 6 ((363,), mpz(120)) 32
> 7 ((1092,), mpz(720)) 64
> 8 ((3279,), mpz(5040)) 128
> 9 ((9840,), mpz(40320)) 256
> 10 ((29523,), mpz(362880)) 512
> 11 ((88572,), mpz(3628800)) 1024
> 12 ((265719,), mpz(39916800)) 2048
> 13 ((797160,), mpz(479001600)) 4096
> 14 ((2391483,), mpz(6227020800L)) 8192
> 15 ((7174452,), mpz(87178291200L)) 16384
> 16 ((21523359,), mpz(1307674368000L)) 32768
> 17 ((64570080,), mpz(20922789888000L)) 65536
> 18 ((193710243,), mpz(355687428096000L)) 131072
> 19 ((581130732,), mpz(6402373705728000L)) 262144
> 20 ((1743392199,), mpz(121645100408832000L)) 524288
> 21 ((5230176600L,), mpz(2432902008176640000L)) 1048576
>
> To replicate this you would need development version of SymPy + an additional
> branch.
> In the output you can see `n`, the leading term and number of non-zero terms
> of `f_{n-1}`. This implies that `f_21` has 2 million non-zero terms. This
> computation took 0.5 GiB of RAM. The case n=2**2000 won't fit in any memory.
>
> This polynomial has specific structure, so it should be possible to derive a
> formula for the terms of `f_n`. SymPy should help you with experimentation
> with `f_n` for small `n` (as you did already).
Maybe we should implement Product.coeff(), which computes this.
Aaron Meurer
>
>
> Thanks again.
>
>
>
>
> On Apr 16, 4:41 pm, "Aaron S. Meurer" <[email protected]> wrote:
> > On Apr 16, 2011, at 11:18 AM, Mateusz Paprocki wrote:
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > > Hi,
> >
> > > On 16 April 2011 18:47, refp16 <[email protected]> wrote:
> > > How do I extract the nth term of an expanded polynomial that is
> > > ordered in a particular way (increasing exponents)?
> >
> > > For example, from:
> >
> > > 1 + x**3 + 2*x**9 + 2*x**12 + 3*x**27 + 3*x**30 + 6*x**36 + 6*x**39
> >
> > > I want the 7th term.
> >
> > > I have tried args but this gives me a tuple that is not ordered in the
> > > same way as its polynomial:
> >
> > > (1, 2*x**12, 3*x**30, 6*x**36, 3*x**27, x**3, 6*x**39, 2*x**9)
> >
> > > My polynomial is type <class 'sympy.core.add.Add'> and I built it
> > > using the following code:
> >
> > > from sympy import *
> > > x,y,z = symbols('xyz')
> >
> > > aux = 1 + x**3
> > > for n in range(2, 4):
> > > baseTerm = 1 + n*x**3**n
> > > aux = aux * baseTerm
> > > auxex = aux.expand()
> > > print auxex
> >
> > > I also tried methods of the Poly class, but I believe they don´t work
> > > precisely because my polynomial is not a Poly object. Would one way be
> > > converting to Poly and then using those methods? How would I convert?
> >
> > > I will assume you use SymPy 0.6.7 (the procedure will be a little
> > > different for development version of SymPy). You can extract terms of
> > > polynomials in two ways:
> >
> > > In [1]: aux = 1 + x**3
> >
> > > In [2]: for n in range(2, 4):
> > > ...: baseTerm = 1 + n*x**3**n
> > > ...: aux = aux * baseTerm
> > > ...:
> > > ...:
> >
> > > In [3]: f = aux.expand()
> >
> > > In [4]: f
> > > Out[4]:
> > > 3 9 12 27 30 36 39
> > > 1 + x + 2⋅x + 2⋅x + 3⋅x + 3⋅x + 6⋅x + 6⋅x
> >
> > > Either using coeff() method from Basic:
> >
> > > In [5]: f.coeff(x**12)
> > > Out[5]: 2
> >
> > > In [6]: f.coeff(x**15) # this returns None, not 0
> >
> > > or convert `f` to an instance of Poly class and also use coeff():
> >
> > > In [7]: g = Poly(f, x)
> >
> > > In [8]: g.coeff(12)
> > > Out[8]: 2
> >
> > > In [9]: g.coeff(15) # now this is really zero
> > > Out[9]: 0
> >
> > This works only if you know the power of the polynomial you want.
> >
> > Strictly speaking, the way to get the nth term in the order that the
> > expression is printed is to do something like
> >
> > In [285]: sorted(auxex.args, cmp=Basic.compare_pretty)[6]
> > Out[285]:
> > 36
> > 6⋅x
> >
> > If you're using Poly, this is easier
> >
> > In [290]: Poly(auxex).terms()[-7]
> > Out[290]: ((36,), 6)
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > > The latter code takes me to another question: Sum is the sum operator,
> > > but don´t we have a product operator?
> >
> > > There are both:
> >
> > > In [1]: sum(k, (k, 1, n))
> > > Out[1]:
> > > 2
> > > n n
> > > ─ + ──
> > > 2 2
> >
> > > In [2]: product(k, (k, 1, n))
> > > Out[2]: n!
> >
> > > In [3]: sum(k, (k, 1, 10))
> > > Out[3]: 55
> >
> > > In [4]: product(k, (k, 1, 10))
> > > Out[4]: 3628800
> >
> > > However, if you have to add or multiply elements of an iterable
> > > container, then it may be simpler (and much faster) to appropriate
> > > constructors directly:
> >
> > > In [5]: Add(x, y, z)
> > > Out[5]: x + y + z
> >
> > > In [6]: Mul(x, y, z)
> > > Out[6]: x⋅y⋅z
> >
> > If you have these in a list, you should do it like
> >
> > In [291]: a = [x, y, z]
> >
> > In [292]: Add(*a)
> > Out[292]: x + y + z
> >
> > In [293]: Mul(*a)
> > Out[293]: x⋅y⋅z
> >
> > Aaron Meurer
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > > Thank you.
> >
> > > --
> > > 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
> > > athttp://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
> > > athttp://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.
