Add compatibility with other mathematical software
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.
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.
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