More tests. Added tarce und norm for AlgebraicNumbers and NumberFields (#23)

Author: KlausC

src/CommutativeRings.jl

@@ -38,7 +38,7 @@ export groebnerbase, SPOL, lextend
 export varnames, varname, factors
 
 export characteristic_polynomial, minimal_polynomial, field_polynomial, adjugate, companion
-
+export field_matrix, field_polynomial, tr, norm
 export compose_inv, Li, Ein, lin1p, lin1pe, ver
 
 export rational_normal_form, matrix, transformation, polynomials
@@ -52,7 +52,7 @@ import Base: copy, show, promote_rule, convert, abs, isless, length, iterate, el
 import Base: Rational, numerator, denominator, precision
 
 import Primes: factor, isprime
-import LinearAlgebra: checksquare, det
+import LinearAlgebra: checksquare, det, norm, tr
 
 # Re-exports (of non-Base functions)
 export det, isprime, factor

src/algebraic.jl

@@ -44,7 +44,7 @@ function AlgebraicNumber(a::Complex)
         id = 2
         ca = conj(ca)
     end
-    AlgebraicNumber(x^2 - 2*ra * x + ia^2 + ra^2, [ca, conj(ca)], id, a, NOCHECK)
+    AlgebraicNumber(x^2 - 2 * ra * x + ia^2 + ra^2, [ca, conj(ca)], id, a, NOCHECK)
 end
 
 function AlgebraicNumber(p::UnivariatePolynomial, a = -Inf)
@@ -215,7 +215,7 @@ function lincomb(a::T, b::T, aa, bb) where T<:UnivariatePolynomial
         lincomb(a, aa)
     elseif iszero(aa)
         lincomb(b, bb)
-    elseif a == b
+    elseif a == b # there should be better performance possible in this case - also check if a and b are collinear
         ca = companion(a)
         ea = I(deg(a))
         characteristic_polynomial(kron(ca, ea * aa) + kron(ea * bb, ca))
@@ -231,6 +231,8 @@ end
 function lincomb(p::T, aa) where T<:UnivariatePolynomial
     if isone(aa)
         p
+    elseif iszero(aa)
+        monom(T)
     else
         q = copy(p)
         qq = q.coeff
@@ -465,8 +467,15 @@ function rationalconst(expr::Expr)
             p = rationalconst(args[2])
             isnothing(p) && return p
             return inv(p)
+        elseif fun == :sqrt
+            p = rationalconst(args[2])
+            isnothing(p) && return p
+            sn = isqrt(abs(numerator(p)))
+            sd = isqrt(denominator(p))
+            q = sn // sd
+            return q^2 == p ? q : nothing
         end
-        if fun ∉ (:(+), :(-), :(*), :(//), :(inv), (:sqrt))
+        if fun ∉ (:(+), :(-), :(*), :(//))
             return false
         end
         n = length(args)
@@ -516,21 +525,26 @@ function isalgebraic(expr::Expr)
     end
 end
 
-function AlgebraicNumber(expr::Expr)
+AlgebraicNumber(expr::Expr) = AlgebraicNumber(to_algebraic_or_rational(expr))
+
+function to_algebraic_or_rational(expr::Expr)
+    x = rationalconst(expr)
+    !isnothing(x) && return x
     isalgebraic(expr) || throw(ArgumentError("expression is not an algebraic number"))
     head = expr.head
     if head == :call
         args = expr.args
         fun = args[1]
         if fun == :(^)
-            a = AlgebraicNumber(args[2])
+            a = to_algebraic_or_rational(args[2])
             q = rationalconst(args[3])
+            denominator(q) != 1 && (a = AlgebraicNumber(a))
             return a^q
         elseif fun ∈ (:(+), :(-), :(*), :(//), :(/), :inv, :sqrt)
             n = length(args)
             eargs = Vector{Any}(undef, n - 1)
             for j = 2:n
-                eargs[j-1] = AlgebraicNumber(args[j])
+                eargs[j-1] = to_algebraic_or_rational(args[j])
             end
             return evaluate(Val(fun), eargs)
         end
@@ -539,11 +553,58 @@ function AlgebraicNumber(expr::Expr)
         throw(ArgumentError("no call, but $head"))
     end
 end
+to_algebraic_or_rational(a::Integer) = a
 
 evaluate(::Val{:(+)}, args) = +(args...)
 evaluate(::Val{:(-)}, args) = -(args...)
 evaluate(::Val{:(*)}, args) = *(args...)
-evaluate(::Val{:(/)}, args) = //(args...)
-evaluate(::Val{:(//)}, args) = //(args...)
+evaluate(::Val{:(/)}, args) = /(args...)
+evaluate(::Val{:(//)}, args) = /(args...)
 evaluate(::Val{:(inv)}, args) = inv(args...)
-evaluate(::Val{:(sqrt)}, args) = sqrt(args...)
+evaluate(::Val{:(sqrt)}, args) = sqrt(AlgebraicNumber(args[1]))
+
+"""
+    com2(p, q)
+
+calculates `mod(q^i, p)` for i = 0:n and stores coeffs in n*n matrix M and vector V.
+The solution of the linear sytem of equations delivers the coefficients of
+the field polynomial of `q(a)` when p is the minimal polynomial of `a`.
+
+Is an alternative to `det(x - companion(p, q))` which is faster sometimes.
+"""
+function com2(p::UnivariatePolynomial{T}, q) where T
+    n = deg(p)
+    S = typeof(value(p[0]))
+    M = zeros(S, n, n)
+    s = q^0
+    for i = 0:n-1
+        k = deg(s)
+        for k = 0:k
+            M[k+1, i+1] = s[k]
+        end
+        s = mod(s * q, p)
+    end
+    V = zeros(S, n)
+    k = deg(s)
+    for k = 0:k
+        V[k+1] = s[k]
+    end
+    M, V
+end
+
+"""
+    evaluate2(q::UnivariatePolynomial, a::AlgebraicNumber)
+
+Evaluates polynomial `q` for an algebraic number `a` using algorithm `com2`.
+"""
+function evaluate2(q::P, a::AlgebraicNumber) where P<:UnivariatePolynomial
+    p = minimal_polynomial(a)
+    M, V = com2(p, q)
+    m = -P(M \ V) + monom(P, deg(a))
+    AlgebraicNumber(m, q(approx(a)))
+end
+
+field_matrix(a::AlgebraicNumber) = companion(minimal_polynomial(a))
+norm(a::AlgebraicNumber) = minimal_polynomial(a)[0] * (-1)^deg(a)
+tr(a::AlgebraicNumber) = -minimal_polynomial(a)[deg(a)-1]
+discriminant(a::AlgebraicNumber) = discriminant(minimal_polynomial(a))

src/numberfield.jl

@@ -91,15 +91,33 @@ generator(::Type{N}) where N<:NumberField = monom(N)
 
 minimal_polynomial(b::NumberField) = minimal_polynomial(AlgebraicNumber(b))
 
-function field_polynomial(b::N) where N<:NumberField
-    r = b.repr
-    q = value(r)
-    p = modulus(r) # == minimal_polynomial(base(N))
-    x = monom(typeof(p))
-    m = det(x * I - companion(p, q))
-    m
+minimal_polynomial(::Type{N}) where N<:NumberField = minimal_polynomial(base(N))
+
+"""
+    field_matrix(b::NumberField)
+
+Return the field matrix for number field element `b` in the standard basis.
+"""
+function field_matrix(b::N) where N <:NumberField
+    p = minimal_polynomial(N)
+    q = value(b.repr)
+    companion(p, q)
 end
 
+"""
+    field polynomial(b::NumberField)
+
+The characteristic polynomial of any field matrix of an element `b` of a number field
+over `a`.
+
+The field matrix of `b` with the standard basis `a^0, a, ..., a^(n-1)` is the
+`companion(p, q)`, where `p` is the minimal polynomial of `a` and `q` is the
+polynomial representing `b (b = q(a))`.
+"""
+field_polynomial(b::NumberField) = characteristic_polynomial(field_matrix(b))
+tr(b::NumberField) = tr(field_matrix(b))
+norm(b::NumberField) = det(field_matrix(b))
+
 function field_polynomial(a::A, b::B) where {A<:AlgebraicNumber,B<:UnivariatePolynomial}
     field_polynomial(b(monom(NF(a))))
 end

src/univarpolynom.jl

@@ -563,7 +563,7 @@ end
 function content_primpart(p::P) where {T,X,P<:UnivariatePolynomial{QQ{T},X}}
     c = content(p)
     Z = ZZ{T}
-    pp = Z[X](Z.(numerator.((p / c).coeff)), ord(p))
+    pp = Z[X]([Z(numerator(x/c)) for x in p.coeff], ord(p))
     c, pp
 end
 
@@ -777,6 +777,8 @@ If `q` is given, return the companion matrix with respect to `q`. For `q == x` t
 identical to the original definition. If `p` is the minimal polynomial of an algebraic
 number `a`, then `det(x*I - companion(p, q))` is a multiple of the minimal polynomial
 of `q(a)`.
+
+Note, that `companion(p, q) == q(companion(p))`.
 """
 function companion(p::UnivariatePolynomial{T}) where T
     companion(T, p)

test/algebraic.jl

@@ -97,11 +97,15 @@ end
     @test approx(c) ≈ op(approx(a), approx(b))
 end
 
-@testset "Algebraic $op($a,A) arithmetic" for op in (+, -, *, /), a in (11, ZZ(11))
+@testset "Algebraic $op($a,A) arithmetic" for op in (+, -, *, /), a in (0, 11, ZZ(11))
     a = big(a)
     b = AlgebraicNumber(x^3 + x + 1)
     c = op(a, b)
-    @test deg(c) == deg(b)
+    if a != 0 || op in (+, -)
+        @test deg(c) == deg(b)
+    else
+        @test c == 0
+    end
     @test minimal_polynomial(c)(c) == 0
     @test approx(c) ≈ op(approx(a), approx(b))
 end
@@ -134,15 +138,34 @@ end
 end
 
 @testset "Algebraic - expressions" begin
+    expr2 = :(sqrt(-2))
+    @test AlgebraicNumber(expr2) == AlgebraicNumber(x^2 + 2, im)
+
+    exprq = :(sqrt((11 / 12)^2))
+    @test AlgebraicNumber(exprq) == AlgebraicNumber(x - 11//12)
+
     expr5 = :(sqrt(5) + 1)
     @test AlgebraicNumber(expr5) / 4 == cospi(QQ(1 // 5))
 
+    expr217 = :(
+        (
+            2sqrt(17 + 3sqrt(17) - sqrt(34 - 2sqrt(17)) - 2sqrt(34 + 2sqrt(17))) +
+            sqrt(34 - 2sqrt(17)) +
+            sqrt(17) - 1
+        ) / 8
+    )
+    @test AlgebraicNumber(expr217) == cospi(QQ(2 // 17)) * 2
+
     expr17 = :(
         1 - sqrt(17) +
         sqrt(34 - sqrt(68)) +
         sqrt(68 + sqrt(2448) + sqrt(2720 + sqrt(6284288)))
     )
     @test AlgebraicNumber(expr17) / 16 == cospi(QQ(1 // 17))
+
+    a = AlgebraicNumber(:((3^(1 // 3) - 6^(1 // 3) + 12^(1 // 3)) / 3))
+    b = AlgebraicNumber(:((2^(1 // 3) - 1)^(1 // 3)))
+    @test a == b
 end
 
 @testset "exactness of pure imaginary roots" begin
@@ -158,16 +181,19 @@ end
     @test string(a) == "AlgebraicNumber(x^2 + 4, 2, 0.0 + 2.0im)"
 end
 
-@testset "cardano's formula $(x^3+a*x^2+b*x+c)" for (a, b, c) in
-                                                    ((1, -8, -6), (6, 9, 8), rand(-10:10, 3))
+@testset "cardano's formula $(x^3+a*x^2+b*x+c)" for (a, b, c) in (
+    (1, -8, -6),
+    (6, 9, 8),
+    rand(-10:10, 3),
+)
     a, b, c = QQ.((a, b, c))
     pc = x^3 + a * x^2 + b * x + c
     p = b - a^2 / 3
     q = 2 * a^3 / 27 - a * b / 3 + c
 
     da = sqrt(AlgebraicNumber((q / 2)^2 + (p / 3)^3))
     ua = (da - q / 2)^(1 // 3)
-    va = (-p/3) / ua # (-da - q / 2)^(1 // 3)
+    va = (-p / 3) / ua # (-da - q / 2)^(1 // 3)
     e3 = cispi(QQ(2 // 3))
     @test ua * va == -p / 3