Ask Your Question
1

Extracting terms of polynomial of certain powers

asked 2020-09-05 19:35:39 +0100

whatupmatt gravatar image

updated 2020-09-05 19:36:02 +0100

Let's say I have a polynomial f. I want to extract the monomials where every power of variable is divisible by a number, say 5. I have two ideas on how to do this.

Method 1:

R.<x,y,z> = QQ[]
f=x^(5)*y^(2)+z^(2)*x^(2)+y^(5)*x^(10)

The end result should give me y^(5)x^(10) as the first term has a power 2 which is not divisible by 5 and obviously the 2nd term is not. I am stuck here as I want to somewhat say the following.

f.exponents()

gives

[(10, 5, 0), (5, 2, 0), (2, 0, 2)]

Then I want to say

f.exponents()[i] % (5,5,5) == (0,0,0) is False, remove the ith term

There are 2 problems with this. The first is when I do f.exponents(), you see the tuples they give me are rearranged so the "ith" term is no longer what I think it is. For example, (10,0,5) is the 3rd term of f, but they put it as the first tuple. The second thing is f.exponents()[i] % (5,5,5) is not valid (i.e. gives error) but hopefully you get what I want. I basically want to divide component-wise and get all remainder 0. What is the correct code for this and the coding language for "remove the ith term"?

Method 2: I can consider f(x^(1/5),y^(1/5),z^(1/5)) and then check whether f.exponents()[i] belongs to Z^(3)-cartesian product of integers. The issue is again I don't know how to code this and f(x^(1/5),y^(1/5),z^(1/5)) gives error as I lie in the polynomial ring. The polynomial ring doesn't allow fractional exponents. I looked online and maybe I am suppose to work with Puiseux polynomials?

https://stackoverflow.com/questions/5...

edit retag flag offensive close merge delete

3 Answers

Sort by » oldest newest most voted
2

answered 2020-09-05 20:14:26 +0100

slelievre gravatar image

updated 2020-09-13 01:54:09 +0100

Pure n-th power monomials of a polynomial

Two methods to extract the "pure n-th power monomials" in a multivariate polynomial in Sage:

  • "dictionary of monomials and coefficients" corresponds to "method 1" in the question

  • "taking n-th roots of monomials" corresponds to "method 2" in the question

Dictionary of monomials and coefficients

A multivariate polynomial has a method to produce a dictionary whose keys are exponent tuples and whose values are coefficients for the corresponding monomials.

A polynomial ring can in turn convert such a dictionary into a polynomial.

With these conversions in mind, the process becomes easy!

Define a polynomial ring and a polynomial:

sage: R.<x, y, z> = QQ[]
sage: f = x^10*y^5 + x^5*y^2 + x^2*z^2
sage: f
x^10*y^5 + x^5*y^2 + x^2*z^2

Turn it into a dictionary:

sage: d = f.dict()
sage: d
{(10, 5, 0): 1, (5, 2, 0): 1, (2, 0, 2): 1}

Extract a dictionary with the congruence condition:

sage: dd = {ee: c for ee, c in d.items() if all(e % 5 == 0 for e in ee)}
sage: dd
{(10, 5, 0): 1}

Turn the extracted dictionary back into a polynomial:

sage: R(dd)
x^10*y^5

These steps can be combined into a function.

Readable version:

def only_perfect_power_monomials(f, n):
    R = f.parent()
    d = f.dict()
    dd = {ee: c for ee, c in d.items() if all(e % n == 0 for e in ee)}
    return R(dd)

Compact version:

def only_perfect_power_monomials(f, n):
    return f.parent()({ee: c for ee, c in f.dict().items()
                             if all(e % n == 0 for e in ee)})

Use the function:

sage: R.<x, y, z> = QQ[]
sage: f = x^10*y^5 + x^5*y^2 + x^2*z^2
sage: only_perfect_power_monomials(f, 5)
x^10*y^5

Taking n-th roots of monomials

We cannot take the fifth root of a polynomial variable, but if a monomial happens to be a perfect fifth power, then we can take its fifth root.

Let us illustrate this and then use it.

Define a polynomial ring and a polynomial:

sage: R.<x, y, z> = QQ[]
sage: f = x^10*y^5 + x^5*y^2 + x^2*z^2
sage: f
x^10*y^5 + x^5*y^2 + x^2*z^2

Extract its monomials and give them names:

sage: f.monomials()
[x^10*y^5, x^5*y^2, x^2*z^2]
sage: a, b, c = f.monomials()

Taking fifth roots of the monomials sometimes works:

sage: a^(1/5)
x^2*y

and sometimes fails:

sage: b^(1/5)
Traceback (most recent call last)
...
TypeError: no conversion of this rational to integer
...
During handling of the above exception, another exception occurred:
...
Traceback (most recent call last)
...
ValueError: not a 5th power

sage: c^(1/5)
Traceback (most recent call last)
...
TypeError: no conversion of this rational to integer
...
During handling of the above exception, another exception occurred:
...
Traceback (most recent call last)
...
ValueError: not a 5th power

Function that takes n-th roots of monomials if possible, zero otherwise:

def monomial_wise_nth_root_or_zero(f, n):
    R = f.parent()
    def nth_root_or_zero(m):
        try:
            return m^(1/n)
        except ValueError:
            return R.zero()
    mm = f.monomials()
    return sum([f[m] * nth_root_or_zero(m) for m in mm], R.zero())

Use the function:

sage: R.<x, y, z> = QQ[]
sage: f = x^10*y^5 + x^5*y^2 + x^2*z^2

sage: monomial_wise_nth_root_or_zero(f, 5)
x^2*y

Function that keeps only perfect n-th power monomials (choose your flavour):

def only_perfect_power_monomials(f, n):
    R = f.parent()
    def nth_root_or_zero(m):
        try:
            return m^(1/n)
        except ValueError:
            return R.zero()
    mm = f.monomials()
    return sum([f[m] * nth_root_or_zero(m)^n for m in mm], R.zero())

def only_perfect_power_monomials(f, n):
    R = f.parent()
    def keep_only_if_perfect_nth_power(m):
        try:
            return (m^(1/n))^n
        except ValueError:
            return R.zero()
    mm = f.monomials()
    return sum([f[m] * keep_only_if_perfect_nth_power(m) for m in mm], R.zero())

def only_perfect_power_monomials(f, n):
    R = f.parent()
    def keep_only_if_perfect_nth_power(m):
        try:
            m^(1/n)
            return m
        except ValueError:
            return R.zero()
    mm = f.monomials()
    return sum([f[m] * keep_only_if_perfect_nth_power(m) for m in mm], R.zero())

def only_perfect_power_monomials(f, n):
    R = f.parent()
    result = R.zero()
    for m in f.monomials():
        try:
            m^(1/n)
            result += f[m] * m
        except ValueError:
            pass
    return result

Use these functions:

sage: R.<x, y, z> = QQ[]
sage: f = x^10*y^5 + x^5*y^2 + x^2*z^2

sage: only_perfect_power_monomials(f, 5)
x^10*y^5
edit flag offensive delete link more

Comments

Great, thanks for your help.

whatupmatt gravatar imagewhatupmatt ( 2020-09-11 18:10:36 +0100 )edit

Hello, in this post, I suggested maybe taking the 5th roots. Can I take 5th roots by modifying the exponents in the dictionary? For example, in the end, my dictionary is (10,5,0): 1 . Can I extract the exponent part, divide by 5 to get (2,1,0) : 1 and convert back to polynomial to get x^(2)*y?

whatupmatt gravatar imagewhatupmatt ( 2020-09-11 23:29:38 +0100 )edit

I am asking this because in the multivariate ring, f(x^(1/5),y^(1/5),z^(1/5)) gives error for me as it says there is no 5th root.

whatupmatt gravatar imagewhatupmatt ( 2020-09-11 23:51:55 +0100 )edit

Taking fractional powers would require working with Puiseux polynomials.

So far Sage only has an implementation for univariate Puiseux series.

Progress on implementing multivariate Puiseux polyomials is tracked at

slelievre gravatar imageslelievre ( 2020-09-12 00:42:52 +0100 )edit

Yes, but I was hoping that in the case where things divide out nicely it is possible. Like even though we are composing with fractional power, we end up with integer powers. Like x^(10)y^5, maybe dividing by 4 is bad but dividing all powers by 5 felt reasonable as you end up with integer powers which remains in your polynomial ring.

whatupmatt gravatar imagewhatupmatt ( 2020-09-12 23:20:01 +0100 )edit
1

answered 2020-09-05 20:12:22 +0100

Emmanuel Charpentier gravatar image

What about :

sum([u[0]*u[1] for u in f if all([(v%5)==0 for v in u[1].exponents()[0]])])
x^10*y^5

This kind of brutal one-liner isn't very readable. If you plan to use this kind of trick often, generalizing it in a cleanly written and documented function is probably a better way...

HTH,

edit flag offensive delete link more

Comments

slelievre's answer is much better than mine...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2020-09-05 20:55:10 +0100 )edit
1

answered 2020-09-05 20:08:45 +0100

rburing gravatar image

Polynomial rings have monomial orderings (in this case it is degrevlex by default), and the terms of polynomials are automatically ordered from greatest to smallest (with respect to the monomial ordering), as you can see e.g. with print(f).

You can get the exponents you want like this:

sage: [e for e in f.exponents() if all(e_i % 5 == 0 for e_i in e)]
[(10, 5, 0)]

Or the indices of the respective terms (again, with respect to the monomial ordering):

sage: idc = [k for (k,e) in enumerate(f.exponents()) if all(e_i % 5 == 0 for e_i in e)]; idc
[0]

Or the terms themselves:

sage: sum(prod(t) for (k,t) in enumerate(list(f)) if k in idc)
x^10*y^5

A bit more efficiently:

sage: terms = list(f)
sage: sum(prod(terms[k]) for k in idc)
x^10*y^5
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2020-09-05 19:35:39 +0100

Seen: 2,276 times

Last updated: Sep 13 '20