diff --git a/Project.toml b/Project.toml index eb4f4e0..1c77ea0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "LinearAlgebraX" uuid = "9b3f67b0-2d00-526e-9884-9e4938f8fb88" -version = "0.2.9" +version = "0.2.10" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index 0c5c09b..ae4a3ea 100644 --- a/README.md +++ b/README.md @@ -155,7 +155,7 @@ julia> A*B Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1) - ``` +``` ## Characteristic polynomial @@ -225,7 +225,7 @@ julia> rrefx(A) 0//1 1//1 0//1 1//1 0//1 -883//325 158//325 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 - ``` +``` ## Homogeneous Vectors diff --git a/docs/build/.documenter-siteinfo.json b/docs/build/.documenter-siteinfo.json index c675f40..5a93ce1 100644 --- a/docs/build/.documenter-siteinfo.json +++ b/docs/build/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-08-30T07:38:13","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-08-30T09:14:15","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/docs/build/index.html b/docs/build/index.html index c2c877d..b8b5fe0 100644 --- a/docs/build/index.html +++ b/docs/build/index.html @@ -1,5 +1,5 @@ -LinearAlgebraX · LinearAlgebraX

LinearAlgebraX

This module implements basic linear algebra methods for matrices with exact entries (e.g., Rational{Int} values). The function names typically match the standard ones in Julia but with an x (for "exact") appended.

The functions in this module work for all types of Integer, Rational, Complex{Integer}, Complex{Rational}, and Mod entries in matrices. Other exact numbers may work as well, but are not tested.

As the goal of this module is always to give exact answers and, at the same time, be type stable, the results of many of these functions are big. That is, the detx of an integer matrix returns a BigInt.

Functions

These functions in this module end with the letter x and have the same definitions as their counterparts that do not have an x. For exact types (such as Ints) these functions give exact results.

  • detx – exact determinant
  • cofactor_det– slower exact determinant (via cofactor expansion)
  • nullspacex – exact nullspace
  • rankx – exact rank
  • invx – exact inverse
  • rrefx – row reduced echelon form
  • eye – lovingly restored
  • char_poly – characteristic polynomial
  • permanent – permanent of a square matrix

Examples follow.

Determinant

julia> A = ones(Int,10,10)+eye(Int,10);
+LinearAlgebraX · LinearAlgebraX

LinearAlgebraX

This module implements basic linear algebra methods for matrices with exact entries (e.g., Rational{Int} values). The function names typically match the standard ones in Julia but with an x (for "exact") appended.

The functions in this module work for all types of Integer, Rational, Complex{Integer}, Complex{Rational}, and Mod entries in matrices. Other exact numbers may work as well, but are not tested.

As the goal of this module is always to give exact answers and, at the same time, be type stable, the results of many of these functions are big. That is, the detx of an integer matrix returns a BigInt.

Functions

These functions in this module end with the letter x and have the same definitions as their counterparts that do not have an x. For exact types (such as Ints) these functions give exact results.

  • detx – exact determinant
  • cofactor_det– slower exact determinant (via cofactor expansion)
  • nullspacex – exact nullspace
  • rankx – exact rank
  • invx – exact inverse
  • rrefx – row reduced echelon form
  • eye – lovingly restored
  • char_poly – characteristic polynomial
  • permanent – permanent of a square matrix

Examples follow.

Determinant

julia> A = ones(Int,10,10)+eye(Int,10);
 
 julia> det(A)
 10.999999999999998
@@ -84,15 +84,39 @@
  Mod{11}(0)  Mod{11}(1)  Mod{11}(0)  Mod{11}(0)  Mod{11}(0)
  Mod{11}(0)  Mod{11}(0)  Mod{11}(1)  Mod{11}(0)  Mod{11}(0)
  Mod{11}(0)  Mod{11}(0)  Mod{11}(0)  Mod{11}(1)  Mod{11}(0)
- Mod{11}(0)  Mod{11}(0)  Mod{11}(0)  Mod{11}(0)  Mod{11}(1)
- ```
+ Mod{11}(0)  Mod{11}(0)  Mod{11}(0)  Mod{11}(0)  Mod{11}(1)

Characteristic polynomial

julia> using SimplePolynomials, LinearAlgebra
 
- ## 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 + 5x - 10x^2 + 10x^3 - 5x^4 + x^5

julia> ans == (x-1)^5 true

julia> using Mods

julia> A = rand(Mod{17},4,4) 4×4 Array{Mod{17},2}: Mod{17}(16) Mod{17}(10) Mod{17}(9) Mod{17}(12) Mod{17}(15) Mod{17}(1) Mod{17}(1) Mod{17}(6) Mod{17}(3) Mod{17}(2) Mod{17}(5) Mod{17}(11) Mod{17}(5) Mod{17}(15) Mod{17}(15) Mod{17}(7)

julia> char_poly(A) Mod{17}(1) + Mod{17}(1)x + Mod{17}(16)x^2 + Mod{17}(5)x^3 + Mod{17}(1)x^4

julia> detx(A) Mod{17}(1)


- ## Row reduced echelon form
+julia> x = getx()
+x
 
- ```
- julia> A = rand(Int,4,6) .% 10
+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
+
+julia> using Mods
+
+julia> A = rand(Mod{17},4,4)
+4×4 Array{Mod{17},2}:
+ Mod{17}(16)  Mod{17}(10)   Mod{17}(9)  Mod{17}(12)
+ Mod{17}(15)   Mod{17}(1)   Mod{17}(1)   Mod{17}(6)
+  Mod{17}(3)   Mod{17}(2)   Mod{17}(5)  Mod{17}(11)
+  Mod{17}(5)  Mod{17}(15)  Mod{17}(15)   Mod{17}(7)
+
+julia> char_poly(A)
+Mod{17}(1) + Mod{17}(1)*x + Mod{17}(16)*x^2 + Mod{17}(5)*x^3 + Mod{17}(1)*x^4
+
+julia> detx(A)
+Mod{17}(1)

Row reduced echelon form

 julia> A = rand(Int,4,6) .% 10
 4×6 Array{Int64,2}:
  6   8  0  -6  -5   4
  0  -5  2   0  -3  -4
@@ -118,24 +142,7 @@
  1//1  0//1  0//1  -1//1  0//1   -23//130  -36//65
  0//1  1//1  0//1   1//1  0//1  -883//325  158//325
  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
- ```
-
- ##  Homogeneous Vectors
-
- A point in projective space is represented by a *homogeneous vector*. This is 
- a list of numbers (like an ordinary vector) but two such vectors are equal 
- if and only if one is a nonzero multiple of the other. 
-
- We provide the `HVector` type to represent homogeneous vectors. The entries 
- in an `HVector` are scaled so that the last nonzero coordinate is `1`. 
- (Technically, we should forbid the all zero vector, but we don't implement
- that restriction.)
-
- To create an `HVector` provide either a list (one-dimensional array) of values 
- or a list of arguments:
- ```
- julia> v = HVector([1,-2,3])
+ 0//1  0//1  0//1   0//1  1//1  -379//325  204//325

Homogeneous Vectors

A point in projective space is represented by a homogeneous vector. This is a list of numbers (like an ordinary vector) but two such vectors are equal if and only if one is a nonzero multiple of the other.

We provide the HVector type to represent homogeneous vectors. The entries in an HVector are scaled so that the last nonzero coordinate is 1. (Technically, we should forbid the all zero vector, but we don't implement that restriction.)

To create an HVector provide either a list (one-dimensional array) of values or a list of arguments:

 julia> v = HVector([1,-2,3])
 [1//3 : -2//3 : 1//1]
 
 julia> w = HVector(2,-4,6)
@@ -191,4 +198,4 @@
 3×3 Array{Rational{Int64},2}:
   0//1   4//3  -1//1
  -1//3  -4//3   2//3
- -1//1   0//1   1//1
+ -1//1 0//1 1//1
diff --git a/docs/build/objects.inv b/docs/build/objects.inv index 30402f2..3b823cf 100644 Binary files a/docs/build/objects.inv and b/docs/build/objects.inv differ diff --git a/docs/build/search_index.js b/docs/build/search_index.js index d2523b1..8d7cdcb 100644 --- a/docs/build/search_index.js +++ b/docs/build/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"#LinearAlgebraX","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"This module implements basic linear algebra methods for matrices with exact entries (e.g., Rational{Int} values). The function names typically match the standard ones in Julia but with an x (for \"exact\") appended.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"The functions in this module work for all types of Integer, Rational, Complex{Integer}, Complex{Rational}, and Mod entries in matrices. Other exact numbers may work as well, but are not tested. ","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"As the goal of this module is always to give exact answers and, at the same time, be type stable, the results of many of these functions are big. That is, the detx of an integer matrix returns a BigInt. ","category":"page"},{"location":"#Functions","page":"LinearAlgebraX","title":"Functions","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"These functions in this module end with the letter x and have the same definitions as their counterparts that do not have an x. For exact types (such as Ints) these functions give exact results.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"detx – exact determinant \ncofactor_det– slower exact determinant (via cofactor expansion)\nnullspacex – exact nullspace\nrankx – exact rank\ninvx – exact inverse\nrrefx – row reduced echelon form\neye – lovingly restored\nchar_poly – characteristic polynomial\npermanent – permanent of a square matrix","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"Examples follow.","category":"page"},{"location":"#Determinant","page":"LinearAlgebraX","title":"Determinant","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = ones(Int,10,10)+eye(Int,10);\n\njulia> det(A)\n10.999999999999998\n\njulia> detx(A)\n11\n\njulia> A = rand(Int,20,20) .% 20;\n\njulia> det(A)\n3.3905496651565455e29\n\njulia> detx(A)\n339054966515654744413389494504","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"For certain Mod matrices, there may be noninvertible nonzero elements in which case the Gaussian elimination algorithm may fail. If that happens, detx falls back to using cofactor expansion which may be very slow. Should that happen, a warning is generated.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using Mods\n\njulia> A = rand(Mod{10},5,5)\n5×5 Array{Mod{10},2}:\n Mod{10}(6) Mod{10}(1) Mod{10}(8) Mod{10}(7) Mod{10}(9)\n Mod{10}(6) Mod{10}(4) Mod{10}(6) Mod{10}(9) Mod{10}(0)\n Mod{10}(9) Mod{10}(8) Mod{10}(7) Mod{10}(8) Mod{10}(0)\n Mod{10}(9) Mod{10}(1) Mod{10}(9) Mod{10}(1) Mod{10}(3)\n Mod{10}(5) Mod{10}(4) Mod{10}(5) Mod{10}(9) Mod{10}(0)\n\njulia> detx(A)\n┌ Warning: Using cofactor expansion to calculate determinant; may be very slow.\n└ @ LinearAlgebraX ~/.julia/dev/LinearAlgebraX/src/detx.jl:41\nMod{10}(4)","category":"page"},{"location":"#Nullspace","page":"LinearAlgebraX","title":"Nullspace","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = reshape(collect(1:12),3,4)\n3×4 Array{Int64,2}:\n 1 4 7 10\n 2 5 8 11\n 3 6 9 12\n\njulia> nullspacex(A)\n4×2 Array{Rational{BigInt},2}:\n 1//1 2//1\n -2//1 -3//1\n 1//1 0//1\n 0//1 1//1\n\njulia> nullspace(A)\n4×2 Array{Float64,2}:\n -0.475185 -0.272395\n 0.430549 0.717376\n 0.564458 -0.617566\n -0.519821 0.172585","category":"page"},{"location":"#Rank","page":"LinearAlgebraX","title":"Rank","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"Consider the 12-by-12 Hibert matrix, H (see hilbert.jl in the extras folder):","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"12×12 Array{Rational{Int64},2}:\n 1//1 1//2 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12\n 1//2 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13\n 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14\n 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15\n 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16\n 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17\n 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18\n 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19\n 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20\n 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20 1//21\n 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20 1//21 1//22\n 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20 1//21 1//22 1//23","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"We compare the results of rank (from the LinearAlgebra module) and rankx (in this module):","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> rank(H)\n11\n\njulia> rankx(H)\n12","category":"page"},{"location":"#Inverse","page":"LinearAlgebraX","title":"Inverse","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using Mods\n\njulia> A = rand(Mod{11},5,5)\n5×5 Array{Mod{11},2}:\n Mod{11}(2) Mod{11}(4) Mod{11}(4) Mod{11}(0) Mod{11}(2)\n Mod{11}(9) Mod{11}(4) Mod{11}(5) Mod{11}(1) Mod{11}(10)\n Mod{11}(3) Mod{11}(4) Mod{11}(5) Mod{11}(6) Mod{11}(0)\n Mod{11}(5) Mod{11}(10) Mod{11}(4) Mod{11}(5) Mod{11}(4)\n Mod{11}(9) Mod{11}(10) Mod{11}(7) Mod{11}(8) Mod{11}(9)\n\njulia> B = invx(A)\n5×5 Array{Mod{11},2}:\n Mod{11}(4) Mod{11}(5) Mod{11}(0) Mod{11}(6) Mod{11}(8)\n Mod{11}(7) Mod{11}(4) Mod{11}(9) Mod{11}(10) Mod{11}(3)\n Mod{11}(6) Mod{11}(0) Mod{11}(2) Mod{11}(5) Mod{11}(5)\n Mod{11}(3) Mod{11}(4) Mod{11}(9) Mod{11}(10) Mod{11}(10)\n Mod{11}(9) Mod{11}(9) Mod{11}(0) Mod{11}(8) Mod{11}(9)\n\njulia> A*B\n5×5 Array{Mod{11},2}:\n Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0)\n Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0)\n Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0)\n Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0)\n Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1)\n ```\n\n ## Characteristic polynomial\n","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using SimplePolynomials, LinearAlgebra","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> x = getx() x","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"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","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> char_poly(A) -1 + 5x - 10x^2 + 10x^3 - 5x^4 + x^5","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> ans == (x-1)^5 true","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using Mods","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = rand(Mod{17},4,4) 4×4 Array{Mod{17},2}: Mod{17}(16) Mod{17}(10) Mod{17}(9) Mod{17}(12) Mod{17}(15) Mod{17}(1) Mod{17}(1) Mod{17}(6) Mod{17}(3) Mod{17}(2) Mod{17}(5) Mod{17}(11) Mod{17}(5) Mod{17}(15) Mod{17}(15) Mod{17}(7)","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> char_poly(A) Mod{17}(1) + Mod{17}(1)x + Mod{17}(16)x^2 + Mod{17}(5)x^3 + Mod{17}(1)x^4","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> detx(A) Mod{17}(1)","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"\n ## Row reduced echelon form\n\n ```\n julia> A = rand(Int,4,6) .% 10\n4×6 Array{Int64,2}:\n 6 8 0 -6 -5 4\n 0 -5 2 0 -3 -4\n 0 -4 2 -8 7 -8\n 1 -3 7 2 -6 2\n\njulia> c = A[:,1] + A[:,2] - A[:,3]\n4-element Array{Int64,1}:\n 14\n -7\n -6\n -9\n\njulia> A = [c A]\n4×7 Array{Int64,2}:\n 14 6 8 0 -6 -5 4\n -7 0 -5 2 0 -3 -4\n -6 0 -4 2 -8 7 -8\n -9 1 -3 7 2 -6 2\n\njulia> rrefx(A)\n4×7 Array{Rational{Int64},2}:\n 1//1 0//1 0//1 -1//1 0//1 -23//130 -36//65\n 0//1 1//1 0//1 1//1 0//1 -883//325 158//325\n 0//1 0//1 1//1 1//1 0//1 551//650 512//325\n 0//1 0//1 0//1 0//1 1//1 -379//325 204//325\n ```\n\n ## Homogeneous Vectors\n\n A point in projective space is represented by a *homogeneous vector*. This is \n a list of numbers (like an ordinary vector) but two such vectors are equal \n if and only if one is a nonzero multiple of the other. \n\n We provide the `HVector` type to represent homogeneous vectors. The entries \n in an `HVector` are scaled so that the last nonzero coordinate is `1`. \n (Technically, we should forbid the all zero vector, but we don't implement\n that restriction.)\n\n To create an `HVector` provide either a list (one-dimensional array) of values \n or a list of arguments:\n ```\n julia> v = HVector([1,-2,3])\n[1//3 : -2//3 : 1//1]\n\njulia> w = HVector(2,-4,6)\n[1//3 : -2//3 : 1//1]\n\njulia> v==w\ntrue","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"Entries in an HVector can be retrieved individually, or the entire vector can be converted to an array:","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> v[2]\n-2//3\n\njulia> Vector(v)\n3-element Array{Rational{Int64},1}:\n 1//3\n -2//3\n 1//1","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"However, entries cannot be assigned:","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> v[2] = 3//4\nERROR: MethodError: no method matching setindex!(::HVector{Rational{Int64}}, ::Rational{Int64}, ::Int64)","category":"page"},{"location":"#Operations-for-HVectors","page":"LinearAlgebraX","title":"Operations for HVectors","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"The product of a matrix and a homogeneous vector is a homogeneous vector:","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = rand(Int,3,3) .% 5\n3×3 Array{Int64,2}:\n -1 0 0\n 3 0 -2\n 3 -1 -2\n\njulia> A*v\n[1//1 : 3//1 : 1//1]\n\njulia> A*Vector(v)\n3-element Array{Rational{Int64},1}:\n -1//3\n -1//1\n -1//3","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"The dot product of two homogeneous vectors is a scalar. Since homogeneous vectors are unchanged by scaling, we only distinguish between zero and nonzero results. Specifically, the dot product is 0 if the two vectors are orthogonal and 1 otherwise. ","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using Mods\n\njulia> u = Mod{3}(1)\nMod{3}(1)\n\njulia> v = HVector(u,u,u)\n[Mod{3}(1) : Mod{3}(1) : Mod{3}(1)]\n\njulia> dot(v,v)\n0\n\njulia> w = HVector(-1,2,1)\n[-1//1 : 2//1 : 1//1]\n\njulia> dot(v,w)\n1","category":"page"},{"location":"#Homogeneous-Matrices","page":"LinearAlgebraX","title":"Homogeneous Matrices","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"We also provide HMatrix to represent a homogeneous matrix. These are constructed by passing an (ordinary) matrix.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = rand(Int,3,3).%5\n3×3 Array{Int64,2}:\n 0 -4 3\n 1 4 -2\n 3 0 -3\n\njulia> HMatrix(A)\nHMatrix: Rational{Int64}[0//1 4//3 -1//1; -1//3 -4//3 2//3; -1//1 0//1 1//1]\n\njulia> Matrix(ans)\n3×3 Array{Rational{Int64},2}:\n 0//1 4//3 -1//1\n -1//3 -4//3 2//3\n -1//1 0//1 1//1","category":"page"}] +[{"location":"#LinearAlgebraX","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"This module implements basic linear algebra methods for matrices with exact entries (e.g., Rational{Int} values). The function names typically match the standard ones in Julia but with an x (for \"exact\") appended.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"The functions in this module work for all types of Integer, Rational, Complex{Integer}, Complex{Rational}, and Mod entries in matrices. Other exact numbers may work as well, but are not tested. ","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"As the goal of this module is always to give exact answers and, at the same time, be type stable, the results of many of these functions are big. That is, the detx of an integer matrix returns a BigInt. ","category":"page"},{"location":"#Functions","page":"LinearAlgebraX","title":"Functions","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"These functions in this module end with the letter x and have the same definitions as their counterparts that do not have an x. For exact types (such as Ints) these functions give exact results.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"detx – exact determinant \ncofactor_det– slower exact determinant (via cofactor expansion)\nnullspacex – exact nullspace\nrankx – exact rank\ninvx – exact inverse\nrrefx – row reduced echelon form\neye – lovingly restored\nchar_poly – characteristic polynomial\npermanent – permanent of a square matrix","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"Examples follow.","category":"page"},{"location":"#Determinant","page":"LinearAlgebraX","title":"Determinant","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = ones(Int,10,10)+eye(Int,10);\n\njulia> det(A)\n10.999999999999998\n\njulia> detx(A)\n11\n\njulia> A = rand(Int,20,20) .% 20;\n\njulia> det(A)\n3.3905496651565455e29\n\njulia> detx(A)\n339054966515654744413389494504","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"For certain Mod matrices, there may be noninvertible nonzero elements in which case the Gaussian elimination algorithm may fail. If that happens, detx falls back to using cofactor expansion which may be very slow. Should that happen, a warning is generated.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using Mods\n\njulia> A = rand(Mod{10},5,5)\n5×5 Array{Mod{10},2}:\n Mod{10}(6) Mod{10}(1) Mod{10}(8) Mod{10}(7) Mod{10}(9)\n Mod{10}(6) Mod{10}(4) Mod{10}(6) Mod{10}(9) Mod{10}(0)\n Mod{10}(9) Mod{10}(8) Mod{10}(7) Mod{10}(8) Mod{10}(0)\n Mod{10}(9) Mod{10}(1) Mod{10}(9) Mod{10}(1) Mod{10}(3)\n Mod{10}(5) Mod{10}(4) Mod{10}(5) Mod{10}(9) Mod{10}(0)\n\njulia> detx(A)\n┌ Warning: Using cofactor expansion to calculate determinant; may be very slow.\n└ @ LinearAlgebraX ~/.julia/dev/LinearAlgebraX/src/detx.jl:41\nMod{10}(4)","category":"page"},{"location":"#Nullspace","page":"LinearAlgebraX","title":"Nullspace","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = reshape(collect(1:12),3,4)\n3×4 Array{Int64,2}:\n 1 4 7 10\n 2 5 8 11\n 3 6 9 12\n\njulia> nullspacex(A)\n4×2 Array{Rational{BigInt},2}:\n 1//1 2//1\n -2//1 -3//1\n 1//1 0//1\n 0//1 1//1\n\njulia> nullspace(A)\n4×2 Array{Float64,2}:\n -0.475185 -0.272395\n 0.430549 0.717376\n 0.564458 -0.617566\n -0.519821 0.172585","category":"page"},{"location":"#Rank","page":"LinearAlgebraX","title":"Rank","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"Consider the 12-by-12 Hibert matrix, H (see hilbert.jl in the extras folder):","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"12×12 Array{Rational{Int64},2}:\n 1//1 1//2 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12\n 1//2 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13\n 1//3 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14\n 1//4 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15\n 1//5 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16\n 1//6 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17\n 1//7 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18\n 1//8 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19\n 1//9 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20\n 1//10 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20 1//21\n 1//11 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20 1//21 1//22\n 1//12 1//13 1//14 1//15 1//16 1//17 1//18 1//19 1//20 1//21 1//22 1//23","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"We compare the results of rank (from the LinearAlgebra module) and rankx (in this module):","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> rank(H)\n11\n\njulia> rankx(H)\n12","category":"page"},{"location":"#Inverse","page":"LinearAlgebraX","title":"Inverse","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using Mods\n\njulia> A = rand(Mod{11},5,5)\n5×5 Array{Mod{11},2}:\n Mod{11}(2) Mod{11}(4) Mod{11}(4) Mod{11}(0) Mod{11}(2)\n Mod{11}(9) Mod{11}(4) Mod{11}(5) Mod{11}(1) Mod{11}(10)\n Mod{11}(3) Mod{11}(4) Mod{11}(5) Mod{11}(6) Mod{11}(0)\n Mod{11}(5) Mod{11}(10) Mod{11}(4) Mod{11}(5) Mod{11}(4)\n Mod{11}(9) Mod{11}(10) Mod{11}(7) Mod{11}(8) Mod{11}(9)\n\njulia> B = invx(A)\n5×5 Array{Mod{11},2}:\n Mod{11}(4) Mod{11}(5) Mod{11}(0) Mod{11}(6) Mod{11}(8)\n Mod{11}(7) Mod{11}(4) Mod{11}(9) Mod{11}(10) Mod{11}(3)\n Mod{11}(6) Mod{11}(0) Mod{11}(2) Mod{11}(5) Mod{11}(5)\n Mod{11}(3) Mod{11}(4) Mod{11}(9) Mod{11}(10) Mod{11}(10)\n Mod{11}(9) Mod{11}(9) Mod{11}(0) Mod{11}(8) Mod{11}(9)\n\njulia> A*B\n5×5 Array{Mod{11},2}:\n Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0)\n Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0)\n Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0)\n Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0)\n Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1)","category":"page"},{"location":"#Characteristic-polynomial","page":"LinearAlgebraX","title":"Characteristic polynomial","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using SimplePolynomials, LinearAlgebra\n\njulia> x = getx()\nx\n\njulia> A = triu(ones(Int,5,5))\n5×5 Array{Int64,2}:\n 1 1 1 1 1\n 0 1 1 1 1\n 0 0 1 1 1\n 0 0 0 1 1\n 0 0 0 0 1\n\njulia> char_poly(A)\n-1 + 5*x - 10*x^2 + 10*x^3 - 5*x^4 + x^5\n\njulia> ans == (x-1)^5\ntrue\n\njulia> using Mods\n\njulia> A = rand(Mod{17},4,4)\n4×4 Array{Mod{17},2}:\n Mod{17}(16) Mod{17}(10) Mod{17}(9) Mod{17}(12)\n Mod{17}(15) Mod{17}(1) Mod{17}(1) Mod{17}(6)\n Mod{17}(3) Mod{17}(2) Mod{17}(5) Mod{17}(11)\n Mod{17}(5) Mod{17}(15) Mod{17}(15) Mod{17}(7)\n\njulia> char_poly(A)\nMod{17}(1) + Mod{17}(1)*x + Mod{17}(16)*x^2 + Mod{17}(5)*x^3 + Mod{17}(1)*x^4\n\njulia> detx(A)\nMod{17}(1)","category":"page"},{"location":"#Row-reduced-echelon-form","page":"LinearAlgebraX","title":"Row reduced echelon form","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":" julia> A = rand(Int,4,6) .% 10\n4×6 Array{Int64,2}:\n 6 8 0 -6 -5 4\n 0 -5 2 0 -3 -4\n 0 -4 2 -8 7 -8\n 1 -3 7 2 -6 2\n\njulia> c = A[:,1] + A[:,2] - A[:,3]\n4-element Array{Int64,1}:\n 14\n -7\n -6\n -9\n\njulia> A = [c A]\n4×7 Array{Int64,2}:\n 14 6 8 0 -6 -5 4\n -7 0 -5 2 0 -3 -4\n -6 0 -4 2 -8 7 -8\n -9 1 -3 7 2 -6 2\n\njulia> rrefx(A)\n4×7 Array{Rational{Int64},2}:\n 1//1 0//1 0//1 -1//1 0//1 -23//130 -36//65\n 0//1 1//1 0//1 1//1 0//1 -883//325 158//325\n 0//1 0//1 1//1 1//1 0//1 551//650 512//325\n 0//1 0//1 0//1 0//1 1//1 -379//325 204//325","category":"page"},{"location":"#Homogeneous-Vectors","page":"LinearAlgebraX","title":"Homogeneous Vectors","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"A point in projective space is represented by a homogeneous vector. This is a list of numbers (like an ordinary vector) but two such vectors are equal if and only if one is a nonzero multiple of the other. ","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"We provide the HVector type to represent homogeneous vectors. The entries in an HVector are scaled so that the last nonzero coordinate is 1. (Technically, we should forbid the all zero vector, but we don't implement that restriction.)","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"To create an HVector provide either a list (one-dimensional array) of values or a list of arguments:","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":" julia> v = HVector([1,-2,3])\n[1//3 : -2//3 : 1//1]\n\njulia> w = HVector(2,-4,6)\n[1//3 : -2//3 : 1//1]\n\njulia> v==w\ntrue","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"Entries in an HVector can be retrieved individually, or the entire vector can be converted to an array:","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> v[2]\n-2//3\n\njulia> Vector(v)\n3-element Array{Rational{Int64},1}:\n 1//3\n -2//3\n 1//1","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"However, entries cannot be assigned:","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> v[2] = 3//4\nERROR: MethodError: no method matching setindex!(::HVector{Rational{Int64}}, ::Rational{Int64}, ::Int64)","category":"page"},{"location":"#Operations-for-HVectors","page":"LinearAlgebraX","title":"Operations for HVectors","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"The product of a matrix and a homogeneous vector is a homogeneous vector:","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = rand(Int,3,3) .% 5\n3×3 Array{Int64,2}:\n -1 0 0\n 3 0 -2\n 3 -1 -2\n\njulia> A*v\n[1//1 : 3//1 : 1//1]\n\njulia> A*Vector(v)\n3-element Array{Rational{Int64},1}:\n -1//3\n -1//1\n -1//3","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"The dot product of two homogeneous vectors is a scalar. Since homogeneous vectors are unchanged by scaling, we only distinguish between zero and nonzero results. Specifically, the dot product is 0 if the two vectors are orthogonal and 1 otherwise. ","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> using Mods\n\njulia> u = Mod{3}(1)\nMod{3}(1)\n\njulia> v = HVector(u,u,u)\n[Mod{3}(1) : Mod{3}(1) : Mod{3}(1)]\n\njulia> dot(v,v)\n0\n\njulia> w = HVector(-1,2,1)\n[-1//1 : 2//1 : 1//1]\n\njulia> dot(v,w)\n1","category":"page"},{"location":"#Homogeneous-Matrices","page":"LinearAlgebraX","title":"Homogeneous Matrices","text":"","category":"section"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"We also provide HMatrix to represent a homogeneous matrix. These are constructed by passing an (ordinary) matrix.","category":"page"},{"location":"","page":"LinearAlgebraX","title":"LinearAlgebraX","text":"julia> A = rand(Int,3,3).%5\n3×3 Array{Int64,2}:\n 0 -4 3\n 1 4 -2\n 3 0 -3\n\njulia> HMatrix(A)\nHMatrix: Rational{Int64}[0//1 4//3 -1//1; -1//3 -4//3 2//3; -1//1 0//1 1//1]\n\njulia> Matrix(ans)\n3×3 Array{Rational{Int64},2}:\n 0//1 4//3 -1//1\n -1//3 -4//3 2//3\n -1//1 0//1 1//1","category":"page"}] } diff --git a/docs/src/index.md b/docs/src/index.md index 0c5c09b..ae4a3ea 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -155,7 +155,7 @@ julia> A*B Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(0) Mod{11}(1) - ``` +``` ## Characteristic polynomial @@ -225,7 +225,7 @@ julia> rrefx(A) 0//1 1//1 0//1 1//1 0//1 -883//325 158//325 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 - ``` +``` ## Homogeneous Vectors