Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add compatibility with other mathematical software #44

Open
LauraBMo opened this issue Sep 16, 2020 · 2 comments
Open

Add compatibility with other mathematical software #44

LauraBMo opened this issue Sep 16, 2020 · 2 comments

Comments

@LauraBMo
Copy link

LauraBMo commented Sep 16, 2020

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 LauraBMo changed the title Add compativility with other mathematical software Add compatibility with other mathematical software Sep 17, 2020
@chakravala
Copy link
Owner

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.

@LauraBMo
Copy link
Author

LauraBMo commented Sep 18, 2020

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants