nemocas / abstractalgebra.jl Goto Github PK

There are no examples of computing popov form in docs/src/matrix.md, despite almost every other things being documented for that module. Also, the popov form docs have been inserted in a random location instead of following the order in src/generic/Matrix.jl. The latter would really help us to keep tabs on what is documented and what is missing, at least when it comes to the ordinary module specific functions that are implemented (excluding funky things like constructors, promote rules, etc.)

I would like to have a function that computes the total degree of an MPoly. Would you welcome such an addition?

It is clear that we don't have a consistent naming convention. This is especially true in the power series.

Probably polcoeff should be pol_coeff and renormalize should be renormalise.

Should parent(BigInt) and parent(BigFloat) be implemented?

And in a similar fashion should, e.g., PolynomialRing(BigInt, "x") work as an alternative of PolynomialRing(ZZ, "x")?

Is there a particular reason why multivariate polynomials can only be evaluated at values that are elements of their parent ring?

In my opinion it would be desirable if one could evaluate them at values belonging to any commutative algebra over the parent ring.

I think this would require a data structures for algebras. Does it exist or is it supposed to be defined?

AbstractAlgebra.jl intends to provide a generic implementation of finitely presented groups, rings, fields, modules, ideals, in Julia with no C dependencies.

Here is an approximate roadmap for future development:

Basic rings

• Julia integers
• Julia rationals
• Julia floats
• Prime fields GF(p)
• Finite fields (Conway polynomials)
• Number fields (arithmetic only)
• Padic field

Fraction fields

• Basic arithmetic

Residue rings R/(Rf)

• Basic arithmetic

Univariate polynomials

• Basic arithmetic
• Gcd
• Extended gcd
• Resultant (Ducos)
• Evaluation
• Composition
• Substitution
• Interpolation (over integral domain)
• Signature
• Half GCD
• Multimodular GCD over Euclidean domain
• Resultant over non-integral domain

Laurent polynomials

• Basic arithmetic

Relative power series

• Basic arithmetic
• Composition
• Reversion

Laurent series

• Basic arithmetic
• Division over a ring
• Composition
• Reversion

Puiseux series

• Basic arithmetic
• Composition
• Reversion

Lazy power series

• Basic arithmetic
• Composition
• Reversion

Multivariate polynomials

• Reverse order of terms and exponents to match Flint
• Basic arithmetic
• Gcd (psr)
• Resultant (Ducos)
• Evaluation at all variables
• [] Evaluation at one variable
• Composition
• Substitution

Multivariate power series

• Basic arithmetic

Multivariate Laurent series

• Basic arithmetic

Multivariate Puiseux series

• Basic arithmetic

Matrices (dense)

• Basic arithmetic
• Solving
• Inverse
• Determinant
• Rank
• Reduced row echelon form
• LU decomposition
• QR decomposition
• Gram matrix
• Hessenberg form
• Hermite normal form
• Smith normal form
• Popov form
• Frobenius form
• Jordan canonical form
• Rational canonical form
• Kronecker product
• Characteristic polynomial
• Minimal polynomial

Matrices (sparse)

• Basic arithmetic
• Solving
• Inverse
• Determinant
• Rank
• Reduced row echelon form
• LU decomposition
• QR decomposition
• Gram matrix
• Hessenberg form
• Hermite normal form
• Smith normal form
• Popov form
• Frobenius form
• Jordan canonical form
• Rational canonical form
• Kronecker product
• Characteristic polynomial
• Minimal polynomial

Groups

• Permutation groups
• Characters of S_n
• Young tableaux

Function fields

• Rational function fields
• Finite extensions

Maps and morphisms

• Composition
• Inverse
• Identity
• Functional maps
• Linear morphisms
• Partial maps

Modules

• Free modules
• Direct products

Finitely presented modules over euclidean domains

• Element arithmetic
• Submodules
• Quotient modules
• Intersection
• Homomorphisms
• kernels
• Images
• Module structural decomposition

Modules over Dedekind domains

• Element arithmetic
• Homomorphisms

Ideals

• Element arithmetic
• Quotient rings?
• Colon module [M:I]

Polynomials over non-commutative ring

• Basic arithmetic

Matrix algebras

• Basic arithmetic

What should the output of f.exps be if f is an MPoly that is zero? I am surprised by the output in the following example:

julia> using AbstractAlgebra

Welcome to AbstractAlgebra version 0.0.8

AbstractAlgebra comes with absolutely no warranty whatsoever


julia> R,vars_R = PolynomialRing(QQ,["x"]);

julia> x = vars_R[1]
x

julia> f = 0*x
0//1

julia> f.exps
1×1 Array{UInt64,2}:
 0x0000000000000003

We probably need to define divmod and use it everywhere.

AbstractAlgebra exports a function IdealSet, but this is the name of a type in Singular.jl, so it gets overwritten in Singular.jl.

As @wbhart proposed in #99 there should be an interface (actually two) to access the coefficients and exponents of multivariate polynomials.

In fact there will be two interfaces, one for polynomials which only define an iterator (because they are implemented as linked lists), and one for polynomials that have random access to terms (because they are implemented as an array).

Such an interface would allow to implement functions such as those in PR #102 more generically.

I have exported popov_with_trafo, since it was included in the user facing docs, but there seems to be no call to it in the test code. It seems untested.

I finally updated my old code to Nemo/AbstractAlgebra split and found that for p::perm, A::MatElem the standard action p*A*inv(p) no longer works, i.e. the method in the title is not (no longer?) available

doc"""
    *(x::AbstractAlgebra.MatElem, P::Generic.perm)
> Apply the pemutation $P$ to the columns of the matrix $x$ and return the result.
"""
function *(x::AbstractAlgebra.MatElem, P::Generic.perm)
   z = similar(x)
   m = rows(x)
   n = cols(x)
   for i = 1:m
      for j = 1:n
         z[i, j] = x[i, P[j]]
      end
   end
   return z
end

should I open pull, or there were some reasons for removal?

set_length! for polynomials currently doesn't return the mutated polynomial. Therefore it cannot be used with immutable types.

The documentation for the Univariate Polynomials interface should be fixed when this is rectified.

The interface for multivariate polys describes functions exponent and exponent! for sparse distributed formats with random access terms. We should implement these functions for generic polynomials, as they are already available for Singular polynomials.

In addition, we should implement functions for getting the coefficient and exponent vector of a term specified by its monomial (rather than its index in the array), and for setting the coefficient of such. This should be documented in the main, required interface for multivariate polynomials.

We should also decide how indexing of terms in a sparse distributed, random access polynomial should work. At present, the first term is given index 0, but this seems inconsistent with Julia's convention on array indices. On the other hand, it is consistent with the indexing of terms of a univariate, except that they are numbered from the least significant, rather than most significant term (for performance reasons in both cases). Also, the i-th term there is the degree i term, whereas there is no such meaning for multivariates.

Base changed the following names:
lufact is now lu and similar for lufact!. And trace is now called tr. Should we follow?

In Julia 1.0.0 the operation A' now returns a matrix of type Adjoint instead of the same type as the original matrix. The function transpose now returns a type Transpose. Therefore we can't use either of them to actually transpose in the implementation of the parent call R([1, 2, 3, 4, 5, 6, 7, 8, 9]) where R is a 3x3 matrix space.

It makes sense to require that all Euclidean domains provide a Euclidean function (though I'm not sure what to do when there is more than one).

It's not immediately clear if this can be done generically, for example, given a Euclidean division function one can't easily retrieve the Euclidean function. And given a Euclidean function, it is obviously possible to compute a GCD, but this might not even be the most efficient generic way of doing it, which raises the question: how does the implementation decide whether to use the Euclidean function or use the gcd function which is currently implemented in terms of the required Euclidean division function.

So the metaphysical question here is whether a Euclidean domain on a computer is a domain that provides a Euclidean division function, or whether it is a domain that provides a Euclidean function.

We should implement the following and document it for (explicit?) Euclidean rings:

ann(a) : return the annihilator ideal of a

The documentation for the general interface claims:

divides(f::MyMPoly{T}, g::MyMPoly{T}) where T <: AbstractAlgebra.RingElem

Return a tuple (flag, q) where flag is true if g
divides f, in which case q will be the exact quotient, or flag is false and q is set to zero.

However, I obtain the following:

julia> K = Nemo.GF(2)
Finite field F_2

julia> R,varsR = AbstractAlgebra.Generic.PolynomialRing( K, ["x", "y"])
(Multivariate Polynomial Ring in x, y over Finite field F_2, AbstractAlgebra.Generic.MPoly{AbstractAlgebra.gfelem{Int64}}[x, y])

julia> divides( R(1), R(0) )
ERROR: Division by zero in divides_monagan_pearce
Stacktrace:
 [1] divides_monagan_pearce(::AbstractAlgebra.Generic.MPoly{AbstractAlgebra.gfelem{Int64}}, ::AbstractAlgebra.Generic.MPoly{AbstractAlgebra.gfelem{Int64}}, ::Int64) at /home/florian/.julia/v0.6/AbstractAlgebra/src/generic/MPoly.jl:1614
 [2] divides(::AbstractAlgebra.Generic.MPoly{AbstractAlgebra.gfelem{Int64}}, ::AbstractAlgebra.Generic.MPoly{AbstractAlgebra.gfelem{Int64}}) at /home/florian/.julia/v0.6/AbstractAlgebra/src/generic/MPoly.jl:1770

Why is hash for perms so complicated? In other projects I'd use something like simple hash(a::perm, h::UInt) = hash(a.d, hash(perm, h)) delegating the hashing to standard routines (note that this does not account for type parametrizing perms, but our equality does without the type as well).

I guess this is a left-over of the flint-based implementation, but can anyone more experienced explain?

I had a feeling we had documentation for these somewhere, but I don't even see them in Nemo.jl.

Perhaps @kalmarek knows what has happened here.

I also noticed that eye(::perm) was documented, but didn't exist. The PermGroups documentation also says there's lots of ways to construct permutations, then doesn't give any. I added an obvious one to the docs, but there are surely more.

It would be theoretically possible to allow ZZ[["x"]], say, by overloading

getindex(::AbstractAlgebra.Integers{BigInt}, ::Array{String,1})

I think in the matrices should work as Julia matrices as much as possible. In particular I think rows(A) should be replaced by size(A,1) and cols(A) by size(A,2). This makes it easier to write generic code that works for both type of matrices.

The type hierarchy diagram in the help still refers to Nemo, and is out-of-date:

https://nemocas.github.io/AbstractAlgebra.jl/types.html

I am surprised to obtain true as the result of the comparison in the following example. I would either expect an error or false. Is true supposed to be the correct result? If so, why?

julia> using AbstractAlgebra
INFO: Recompiling stale cache file /home/florian/.julia/lib/v0.6/AbstractAlgebra.ji for module AbstractAlgebra.

Welcome to AbstractAlgebra version 0.0.8

AbstractAlgebra comes with absolutely no warranty whatsoever


julia> R,vars_R = PolynomialRing(QQ,["x","y"]; ordering=:deglex); x,y = vars_R;

julia> S,vars_S = PolynomialRing(QQ,["a","b"]; ordering=:deglex); a,b = vars_S;

julia> a==x
true

Executing docs/make.jl leads to the following warning:

WARNING: julia keyword to Documenter.deploydocs() not specified.
In the future julia will be a required keyword argument
instead of defaulting to "nightly".

So I think one should add the julia keyword there. I am not sure which one to use. "nightly"? "0.6"?

We currently have constructors to build a generic power series from an array of coefficients, but the corresponding constructors for arrays of integers appear to be missing.

The generic Smith normal form is not documented in docs/src/matrix.md.

We currently have a section in the docs for multivariates with random access, but not an interface for multivariates with iterator only (such as Singular.jl provides).

For elements of type MPolyElem{T} the method nvars is implemented, but vars is not, so one has to use vars(parent(p)). I think either both or neither should be implemented. (I think it would be nice to have vars available.)

I would consider it desirable to be able to convert a polynomial of type MPoly{T}, which is in fact a polynomial in only one variable (even if its parent ring has more variables), to an element of type PolyElem{T}.

We should implement the following and document it for the general Euclidean interface:

g, e,f, x, y = gcdxx(a,b)

g = gcd(a,b) = ea+fb
1 = ey-fx

Since the names of the variables of polynomials are stored and retrieved (using the vars function) as Symbols, would it be better for the constructors of polynomial rings to also use Symbols instead of Strings? I think it is more consistent.

I am surprised to find that after

using AbstractAlgebra
R,v = PolynomialRing(ZZ,["x","y"]); x=v[1]

the output of

x in Set(gens(R))

seems not to be deterministic, as the following example shows:

julia> for i=1:20 println(x in Set(gens(R))) end
false
false
false
false
false
false
true
false
false
false
true
false
false
true
false
false
false
false
false
false

I am not sure whether this behavior is indented or whether it is a bug in either Julia or AbstractAlgebra.

After

using Nemo; using AbstractAlgebra
R,varsR = PolynomialRing(Nemo.ArbField(64),["x","y"]); x,y=varsR

the output of

x + Nemo.ArbField(64)(-sqrt(2)) * y

is

x[-1.4142135623730951455 +/- 2.54e-20]*y

I think a plus should be printed between x and the coefficient of y, because I would interpret the current output to be a monomial.

I would like to suggest a function that takes a matrix as input and returns its minors. In fact, I have an implementation, but for convenience I used the function subsets from IterTools.jl.

Is such a function welcome? Is a dependency of IterTools.jl acceptable?

There seems to be a problem with showing elements in AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.MPoly{BigInt}} as the following example illustrates:

using AbstractAlgebra
R, gens_R = PolynomialRing(ZZ,["x"]); x = gens_R[1]
S, gens_S = PolynomialRing(R,["y"]); y = gens_S[1]
show(1+y)
yERROR: MethodError: no method matching monomial_iszero(::Array{UInt64,2}, ::Int64)
Closest candidates are:
  monomial_iszero(::Array{UInt64,2}, ::Int64, ::Int64) at /home/florian/.julia/v0.6/AbstractAlgebra/src/generic/MPoly.jl:150
Stacktrace:
 [1] isnegative(::AbstractAlgebra.Generic.MPoly{BigInt}) at /home/florian/.julia/v0.6/AbstractAlgebra/src/generic/MPoly.jl:529
 [2] show(::Base.TTY, ::AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.MPoly{BigInt}}, ::Array{String,1}) at /home/florian/.julia/v0.6/AbstractAlgebra/src/generic/MPoly.jl:449
 [3] show(::AbstractAlgebra.Generic.MPoly{AbstractAlgebra.Generic.MPoly{BigInt}}) at ./coreio.jl:3

Is there an easy way to obtain the number of variables actually involved in a multivariate polynomial of type MPoly? The function nvars reported the number of variables of the polynomial ring the polynome belongs to.

This is kinda dumb, but it can trip up a new user like me:

julia> R = PolynomialRing(K,'x')
ERROR: MethodError: no method matching PolynomialRing(::Nemo.FlintPadicField, ::Char)

If it's not difficult/bad for some reason allowing a char to be the name of a variable would be nice.

The nullspace documentation in the Generic Matrices section seems to be pulling an example from Base, or something like that. It doesn't return a tuple as the documentation describes.

Currently, divexact for generic Laurent series is implemented by inverting the second series. This is fine over a field, but over a ring, e.g. ZZ, it isn't sufficient. It can be that division is possible, even if the second series isn't invertible.

Unless these functions return the mutated series, it is not possible to implement immutable power series types.

Of course, every use of such a function must be modified if this is changed.

In the Generic matrices docs in docs/src/matrix.md, there are examples of almost every function, except hnf.

We should add a Group interface to the docs. Actually, we should add a general AbstractAlgebra object interface, since it seems like some things have to exist for every AbstractAlgebra object, regardless of what they are, e.g. parent or deepcopy, etc.

We could then factor out some stuff that would be in both the Group and Ring interface.

It's not true that a polynomial is a unit in its polynomial ring iff its constant coefficient is a unit. This has been assumed in the current implementation. For example, it isn't true over Z/nZ.

As it is actually quite expensive in general to compute this, we should look carefully where isunit has been used in generic code and remove references to it where this is not precisely what is required.

We will need abstract types AbstractAlgebra.LaurentSeriesRingElem and AbstractAlgebra.LaurentSeriesFieldElem for the Generic types of the same name to belong to. This is because when we come to implement external C implementations for various rings, we will want their element types to belong to the same abstract types as the generic implementation.

On the other hand, there seems to be no reason to have abstract types for all the different series parent types. They can just belong directly to AbstractAlgebra.Ring and AbstractAlgebra.Field. However, this should first be discussed.

In general, the abstract types exist to lower code duplication (implement the function once for an abstract type, instead of multiple times for each type belonging to it). Ring and Field parent objects generally don't have enough functions implemented for them to waste time with all the different abstract types that the parent types can belong to.

I'm currently creating a MatrixAlgebra type, but identity_matrix and matrix are now ambiguous. Should they create the element of the MatrixAlgebra or the element of the MatrixSpace.

I actually don't see an easy fix for this.

Currently, we can only put a matrix in Hessenberg form if it is square. But we allow non-square matrices in ishessenberg. I'm just wondering if it was intentional, and whether we should return false or an exception if it is not square.

I would like to propose a second form of the function AbstractAlgebra.Generic.sub that takes as input a matrix A, a list of row indices I (integers) and column indices J and returns the submatrix B of A given by B_{i,j} = A_{I[i], J[j]}.

The currently existing function sub(M::AbstractAlgebra.MatElem, r1::Int, c1::Int, r2::Int, c2::Int) would be the special case where I=[i for i=r1:r2] and J=[j for j=c1:c2].