Added char_poly

Project.toml

@@ -1,12 +1,16 @@
 name = "LinearAlgebraX"
 uuid = "9b3f67b0-2d00-526e-9884-9e4938f8fb88"
-version = "0.0.1"
+version = "0.0.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SimplePolynomials = "cc47b68c-3164-5771-a705-2bc0097375a0"
+SimpleRationalFunctions = "1a520dc8-4f4e-523b-a9bd-3b3d46c5454b"
 
 [compat]
 julia = "1"
+SimplePolynomials = "0"
+SimpleRationalFunctions = "0"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -29,6 +29,7 @@ For exact types (such as `Int`s) these functions give exact results.
 * `invx` -- exact inverse
 * `rrefx` -- row reduced echelon form
 * `eye` -- lovingly restored
+* `char_poly` -- characteristic polynomial
 
 ## Examples
 
@@ -133,8 +134,31 @@ julia> A*B
  Mod{11}(0)  Mod{11}(0)  Mod{11}(0)  Mod{11}(0)  Mod{11}(1)
  ```
 
+ #### Characteristic polynomial
+
+```
+julia> using SimplePolynomials, LinearAlgebra
+
+julia> x = getx()
+x
+
+julia> A = triu(ones(Int,5,5))
+5×5 Array{Int64,2}:
+ 1  1  1  1  1
+ 0  1  1  1  1
+ 0  0  1  1  1
+ 0  0  0  1  1
+ 0  0  0  0  1
+
+julia> char_poly(A)
+-1 + 5*x - 10*x^2 + 10*x^3 - 5*x^4 + x^5
+
+julia> ans == (x-1)^5
+true
+```
+
  #### Row reduced echelon form
- 
+
  ```
  julia> A = rand(Int,4,6) .% 10
 4×6 Array{Int64,2}:
@@ -164,3 +188,7 @@ julia> rrefx(A)
  0//1  0//1  1//1   1//1  0//1   551//650  512//325
  0//1  0//1  0//1   0//1  1//1  -379//325  204//325
  ```
+
+ ## To do
+
+ Still having some issues with integer overflow.

src/LinearAlgebraX.jl

@@ -1,5 +1,5 @@
 module LinearAlgebraX
-using LinearAlgebra
+using LinearAlgebra, SimplePolynomials, SimpleRationalFunctions
 
 
 # IntegerX is any sort of real or Gaussian integer
@@ -8,6 +8,8 @@ IntegerX = Union{S,Complex{S}} where S<:Integer
 # RationalX is a Rational or Complex Rational based on integers
 RationalX = Union{Rational{S},Complex{Rational{S}}} where S<:Integer
 
+TypeX = Union{IntegerX, RationalX}
+
 
 function _recip(x::T) where T <: IntegerX
     return 1//x
@@ -22,7 +24,7 @@ include("rrefx.jl")
 include("invx.jl")
 include("rankx.jl")
 include("nullspacex.jl")
-
+include("char_poly.jl")
 
 
 end # module

src/char_poly.jl

@@ -0,0 +1,15 @@
+export char_poly
+
+function char_poly(A::Matrix{T}) where T<:TypeX
+    r,c = size(A)
+    @assert r==c "Matrix must be square"
+
+    xI = zeros(SimplePolynomial,r,r)
+    x = getx()
+    for i=1:r
+        xI[i,i] = x
+    end
+
+    f = detx(xI-A)
+    p = integerize(numerator(f))
+end

test/runtests.jl

@@ -15,3 +15,15 @@ H = hilbert(12)
 A = ones(Int,3,5)
 N = nullspacex(A)
 @test all(0 .== A*N)
+
+
+using SimplePolynomials
+
+A = triu(ones(Int,5,5))
+p = char_poly(A)
+x = getx()
+@test p == (x-1)^5
+
+A = ones(Int,5,5)
+p = char_poly(A)
+@test p == x^4 * (x-5)