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()).
>
> 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
>
> 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).
>
> 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.
