MultivariatePolynomials.jl icon indicating copy to clipboard operation
MultivariatePolynomials.jl copied to clipboard

leading coefficient with respect to a given variable

Open lucaferranti opened this issue 4 years ago • 5 comments

Browsing the source code, I didn't find a function to compute the "leading coefficient" of a multivariate polynomial with respect to a given variable, for example if I have the polynomial x^2*y + x*y the expected behavior would be

leadingcoefficient(x^2*y + x*y, x) = y
leadingcoefficient(x^2*y + x*y, y) = x^2 + x

(sometimes that's also called initial of the polynomial)

this functionality is used e.g. when computing the pseudodivision with multivariate polynomials, which is used e.g. in the Ritt-Wu method of characteristic sets.

Would it be a nice addition to the package?

lucaferranti avatar Aug 09 '21 18:08 lucaferranti

Currently I have been using the following patch, but since I'm not very familiar with the package internals, I suspect my solution is way far from optimal 😄

function LC!(p, var)

    d = maxdegree(p, var)
    d == 0 && return zero(p)

    idx = findfirst(x -> x == var, p.x.vars)
    lc = zero(p)
    @inbounds  for term in p
        if degree(term, var) == d
            term.x.z[idx] = 0
            lc += term
        end
    end
    return lc
end

LC(p, var) = LC!(copy(p), var)

lucaferranti avatar Aug 09 '21 18:08 lucaferranti

bump 😄

lucaferranti avatar Sep 21 '21 17:09 lucaferranti

Here is a one line:

coefficient(leadingterm(isolate_variable(p, var)))

but it's not the fastest as it computes the coefficient for each degree of var, not just the maximal one. This might be faster:

function LC(p, var)
    d = maxdegree(p, var)
    d == 0 && return zero(p)
    polynomial([subs(t, var => 1) for t in terms(p) if degree(t, var) == d)
end

Your suggestion uses internals of DynamicPolynomials and will be slow because of the line: lc += term. It's best to either mutate lc using MutableArithmetics or to build a list of term and then turn it into a polynomial with polynomial.

blegat avatar Sep 21 '21 18:09 blegat

First of all, apologies for never answering you after even bumping the issue, it wasn't very polite from my side. 🤦‍♂️

Anyway, your suggested implementation seemed to give a good speedup to my naive approach! Thanks!

lucaferranti avatar Jan 16 '22 16:01 lucaferranti

No worry, I'll leave the issue open as it might be a nice additional feature of the package.

blegat avatar Jan 16 '22 19:01 blegat