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

Add compatibility with other mathematical software

Open LauraBMo opened this issue 5 years ago • 2 comments

Thanks for all the work in Reduce.jl, it is a great tool. I need to complement it with Polymake.jl (to study some properties of the Newton polytope of a polynomial, it is for CRNT.jl. So, I need the exponents of a polynomial as Julia vectors of integers, and ideally another Julia vector, with the same order, containing the coefficients. For this, I started with the exponent of a term

function extractexponentsofterm(T::RExpr)
    return v->Algebra.deg(T,v)
end

julia> z=R"2*x*y^2+4*x^2-1"

      2        2
4*x  + 2*x*y  - 1
    
julia> extractexponentsofterm(Algebra.lterm(z,:y)).([:x,:y,:z])
3-element Array{RExpr,1}:
 1
 2
 0

But I am not sure how to iterate over all the term of a polynomial, and how to extract the coefficient neither. One can iterate over z:=z-lterm(z) but requiring arithmetic operations seems not optimal.

LauraBMo avatar Sep 16 '20 12:09 LauraBMo

Perhaps it may be beneficial to convert the polynomial to a different representation in Julia and to then use this new representation of a polynomial as a data structure after converting from RExpr representation.

chakravala avatar Sep 17 '20 21:09 chakravala

Thanks for the comment. I do not know whether String is the structure you had in mind... For exponents I came up with:

function extractexponentsofterm(T::RExpr)
    return v->Algebra.deg(T,v)
end

function listofterms(P::RExpr)
    return RExpr.(split(String(P),r"\+|\-"))
end

function exponentsofpoly(P::RExpr, vars)
    return [extractexponentsofterm(t).(vars) for t in listofterms(P)]
end

which seems to work quite well. But for coefficients it seems more complicated. I created, iteratively appling Algebra.coeff, a multidimensional array representing the polynomial. With this actually, I get all the information, the exponents are the coordinates of the non-zero entries and their coefficients are such non-zero entries. But I would say it is overkilling and surely my implementation is very naive (I would really appreciate any comment).

function VectorOfCoeff(P::RExpr, var)
    return collect(RExpr.(parse(Algebra.coeff(P,var))))
end

function degOfPoly(P::RExpr)
    return v-> parse(Algebra.deg(P,v))
end

function PolyToTensor(P::RExpr,vars)
    dims = Tuple(1 .+ degOfPoly(P).(vars))
    N = size(dims,1)-1
    out = zeros(RExpr,dims)
    out[:,fill(1,N)...] = VectorOfCoeff(P,vars[1])
    for n in 1:N
        for I in CartesianIndices(dims[1:n])
            l = VectorOfCoeff(out[(Tuple(I)..., fill(1, N-n+1)...)...],vars[n+1])
            out[(Tuple(I)..., 1:size(l,1), fill(1, N-n)...)...] = l
        end
    end
    return out
end

Some examples:

julia> Algebra.on(:exp)

julia> p = R"2*x*y^2+4*x^2-1"

   2        2
4*x  + 2*x*y  - 1

julia> exponentsofpoly(p,[:x,:y,:w])
3-element Array{Array{RExpr,1},1}:
 [1, 2, 0]
 [2, 0, 0]
 [0, 0, 0]

julia> PolyToTensor(p,[:x,:y,:w])
3×3×1 Array{Int64,3}:
[:, :, 1] =
 -1  0  0
  0  0  2
  4  0  0

julia> q = R"4*x^2+3*x-y^2-w*x*y^2-1"

                2             2
(4*x + 3)*x - (y  + 1) - w*x*y

julia> exponentsofpoly(q,[:x,:y,:w])
5-element Array{Array{RExpr,1},1}:
 [2, 0, 0]
 [1, 0, 0]
 [0, 2, 0]
 [1, 2, 1]
 [0, 0, 0]

julia> PolyToTensor(q,[:x,:y,:w])
3×3×2 Array{Int64,3}:
[:, :, 1] =
 -1  0  -1
  3  0   0
  4  0   0

[:, :, 2] =
 0  0   0
 0  0  -1
 0  0   0

LauraBMo avatar Sep 18 '20 01:09 LauraBMo