Suppose I have, e.g., 250 linear equations with Rational{BigInt} coefficients, with 150 variables. If the system is inconsistent I want to know about it, otherwise I want the exact solution.
Is there an easier way than writing my own Gaussian elimination?
This won’t run often, so I’m not concerned with long runtimes.
Have a look at Nemo.jl; it does exact linear algebra.
Just for sake of completeness, either using AbstractAlgebra or (more efficiently) Nemo, it can be done as follows:
julia> using AbstractAlgebra # will also work with using Nemo;
julia> A = QQ[1 0; 3 0; 5 0]; # create an AbstractAlgebra matrix
julia> v = QQ[1; 2; 3]; # create a 3x1 matrix
julia> can_solve_with_solution(A, v) # solve Ax = v
(false, [1; 0])
julia> A = QQ[1 2; 3 4; 5 6];
julia> can_solve_with_solution(A, v) # solve Ax = v
(true, [0; 1//2])
julia> fl, x = can_solve_with_solution(A, v) # solve Ax = v
(true, [0; 1//2])
julia> A*x == v
true
One can also do can_solve_with_solution(A, v; side = :left) to solve xA = v.
TBH I didn’t even realize that can_solve_with_solution existed, instead i used AbstractAlgebra.rref 
hello thofma,
this interests me very much. does a QQMatrix has some elements of the type BigInt ? I want to solve a very big system A.X=Y with Rational{BigInt}. Thanks !
And, is there a fast way to convert a Matrix{Rational{BigInt}} to a QQMatrix, and back ? I could copy it term by term, eventually …
thanks a lot 
Yes, there is a fast way. Here is an example, where one solves A*X = Y where X and Y are of type Vector.
julia> AA = Rational{BigInt}[1 2; 3 4]
2×2 Matrix{Rational{BigInt}}:
1 2
3 4
julia> A = matrix(QQ, AA)
[1 2]
[3 4]
julia> YY = Rational{BigInt}[5, 6]
2-element Vector{Rational{BigInt}}:
5
6
julia> Y = QQ.(YY)
2-element Vector{QQFieldElem}:
5
6
julia> fl, X = can_solve_with_solution(A, Y; side = :right) # solves A * X = Y for X
(true, QQFieldElem[-4, 9//2])
julia> Rational{BigInt}.(X)
2-element Vector{Rational{BigInt}}:
-4
9//2
Here is an example, where X and Y are themselves of type Matrix:
julia> YY = Rational{BigInt}[5 6; 7 8]
2×2 Matrix{Rational{BigInt}}:
5 6
7 8
julia> Y = matrix(QQ, YY)
[5 6]
[7 8]
julia> fl, X = can_solve_with_solution(A, Y; side = :right) # solves A * X = Y for X
(true, [-3 -4; 4 5])
julia> XX = Matrix{Rational{BigInt}}(X)
2×2 Matrix{Rational{BigInt}}:
-3 -4
4 5
Depending on how what you mean with “very big”, it might take a while to solve your system, but give it a try!
Ok thanks, everything works well I am new to Julia, and I don’t understand “constructors” well !
“Big Matrix” means only that I might need some BigInt, in any case the coefficients explode
Thanks a lot !
Gustave
Hello
I actually I should use Matrix{Complex{Rational{Int}}} . Is there way to do so with Nemo ? Or maybe with AbstractAlgebra ? Thanks !
finally I did the code with Complex{Float64}. But I like Algebra, so that your answer interests me !!!
Bye
Sorry for the late reply. Yes, this is also possible:
julia> QQi = Nemo.QQiField()
Gaussian rational field
julia> AA = Complex{Rational{BigInt}}[1 2im; 3 4im]
2×2 Matrix{Complex{Rational{BigInt}}}:
1//1+0//1*im 0//1+2//1*im
3//1+0//1*im 0//1+4//1*im
julia> A = matrix(QQi, AA)
[1 2*im]
[3 4*im]
julia> YY = Complex{Rational{BigInt}}[5im, 6]
2-element Vector{Complex{Rational{BigInt}}}:
0//1 + 5//1*im
6//1 + 0//1*im
julia> Y = QQi.(YY)
2-element Vector{Nemo.QQiFieldElem}:
5*im
6
julia> fl, X = can_solve_with_solution(A, Y; side = :right) # solves A * X = Y for X
(true, Nemo.QQiFieldElem[6 - 10*im, 15//2 + 3*im])
julia> XX = [Complex{Rational{BigInt}}(real(y), imag(y)) for y in X]
2-element Vector{Complex{Rational{BigInt}}}:
6//1 - 10//1*im
15//2 + 3//1*im
OK thanks, it’s perfect ! I gonna try it tomorrow
Have a nice day
Gustave
It works perfectly well, thank you !!!
hello, once again
I want to solve AX=Y, for A a matrix and Y a vector, for some Rational{Polynomial{Complex{Rational{BigInt}}}}, is it possible ? I found everything in the documentation of AbstractAlgebra, but the Complex{ } part. Any help ?
I’ve also tried with A \ Y, where the awkward type is :
RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}},
Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}}
generated by “@polyvar eta_var”
but strangely it works for a small matrix but not for a big one. I think it changes the algorithm when the size increases, or I didn’t understand what happens there.
It’s frustrating because this is a simple algorithm, but I can’t find it anywhere 
Maybe I could write it by myself.
Thanks for any help !
Gustave
julia> using Nemo
Welcome to Nemo version 0.48.0
Nemo comes with absolutely no warranty whatsoever
julia> QQi = Nemo.GaussianRationals()
Gaussian rational field
julia> QQi(3, 4)
3 + 4*im
julia> QQi(3 + 4im)
3 + 4*im
Sounds like a bug. Can you provide a reproducer, so the bug could be reported?
Thank you !
For the bug, I describe what I use, what works and what does not.
using DynamicPolynomials
@polyvar eta_var
eta = (BigInt(0)//1 + (1//1 + 1im*0//1) *eta_var) / (BigInt(1)//1 +(0//1+0//1 * 1im) *eta_var)
the goal of using eta is to have a type of a rational{Polynomial{Rational{Complex}}
the type of eta is :
RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}, Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}
A = [ eta 2*eta ; 0 eta * eta]
YA = [(1+0 * eta) (1+0 * eta)] '
A \ YA works
B # Is a matrix
21×10 Matrix{RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}, Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}}:
YB # Is a vector
21-element Vector{RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}, Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}}:
B[1,2] = ((-1//8 + 0//1*im) eta_var^2) / ((1//1 + 0//1*im))
B[3,2] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[4,3] = ((-1//4 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[5,3] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[6,3] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[7,4] = ((-1//8 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[8,4] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[10,1] = ((1//1 + 0//1*im)) / ((1//1 + 0//1*im))
B[11,6] = ((-1//4 + 0//1*im)) / ((1//1 + 0//1*im))
B[12,7] = ((-1//4 + 0//1*im)) / ((1//1 + 0//1*im))
B[13,8] = ((-1//8 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[14,8] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[16,9] = ((-1//4 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[17,9] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[18,9] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[19,10] = ((-1//8 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[21,10] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
YB[12] = ((1//1 + 0//1*im)) / ((1//1 + 0//1*im))
(mind to use eta and not eta_var, while recopying, if you do. And add the * )
B \ YB gives a bug :
ERROR: MethodError: no method matching abs2(::RationalPoly{Polynomial{…}, Polynomial{…}})
The function `abs2` exists, but no method is defined for this combination of argument types.
Stacktrace:
[1] _qreltype(::Type{RationalPoly{Polynomial{DynamicPolynomials.Commutative{…}, Graded{…}, Complex{…}}, Polynomial{DynamicPolynomials.Commutative{…}, Graded{…}, Complex{…}}}})
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:341
[2] qr(A::Matrix{RationalPoly{Polynomial{…}, Polynomial{…}}}, arg::ColumnNorm; kwargs::@Kwargs{})
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:424
[3] qr(A::Matrix{RationalPoly{Polynomial{…}, Polynomial{…}}}, arg::ColumnNorm)
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:422
[4] \(A::Matrix{RationalPoly{Polynomial{…}, Polynomial{…}}}, B::Vector{RationalPoly{Polynomial{…}, Polynomial{…}}})
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1134
for the 0 elments, you should use :
0 + 0*eta
I hope it’s okay for you
Could you edit your message to use code block formatting? Instead of this:
> code
… do this:
```julia
code
```
That will preserve the * symbols. EDIT: thanks!
I realize now this isn’t a bug, it’s just \ behaving as documented. It’s doc string says that overdetermined (actually all non-square) systems are solved via:
the minimum-norm least squares solution computed by a pivoted QR factorization
… so not what you want.
I suggest trying with Nemo.jl.
To elaborate on what @nsjako wrote. Here is how you can do it using Nemo.jl:
julia> QQi = Nemo.GaussianRationals();
julia> QQieta, eta_var = QQi[:eta];
julia> iim = QQi(0, 1);
julia> B = zero_matrix(QQieta, 21, 10);
julia> YB = zeros(QQieta, 21);
julia> begin
B[1,2] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[3,2] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[4,3] = ((-1//4 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[5,3] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[6,3] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[7,4] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[8,4] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[10,1] = ((1//1 + 0//1*iim)) / ((1//1 + 0//1*iim))
B[11,6] = ((-1//4 + 0//1*iim)) / ((1//1 + 0//1*iim))
B[12,7] = ((-1//4 + 0//1*iim)) / ((1//1 + 0//1*iim))
B[13,8] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[14,8] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[16,9] = ((-1//4 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[17,9] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[18,9] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[19,10] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[21,10] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
YB[12] = QQieta(((1//1 + 0//1*iim)) / ((1//1 + 0//1*iim)))
end;
julia> solve(B, YB; side = :right)
10-element Vector{AbstractAlgebra.Generic.Poly{Nemo.QQiFieldElem}}:
0
0
0
0
0
0
-4
0
0
0
If you want to go back and forth between Complex{Rational{BigInt}} and the Nemo type, you can do
julia> a = QQi(2, 3);
julia> Complex(Rational(real(a)), Rational(imag(a)))
2//1 + 3//1*im