Merge pull request #24 from KlausC:krc/flint

use FLINT bigint library to replace `BigInt`.

Author: KlausC

Project.toml

@@ -6,12 +6,14 @@ repo = "https://github.com/KlausC/CommutativeRings.jl"
 version = "0.6.0"
 
 [deps]
+FLINT_jll = "e134572f-a0d5-539d-bddf-3cad8db41a82"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 Aqua = "0.8"
+FLINT_jll = "300.100"
 LinearAlgebra = "1"
 Polynomials = "1,2,3,4"
 Primes = "0.5"

src/CommutativeRings.jl

@@ -4,10 +4,17 @@ using LinearAlgebra
 using Base.Checked
 using Primes
 using Random
+using FLINT_jll
+
+const libflint = FLINT_jll.libflint
+
+function flint_abort()
+  error("Problem in the Flint-Subsystem")
+end
 
 export category_trait, isfield
 export Ring, RingInt, FractionRing, QuotientRing, Polynomial
-export ZZ, QQ, ZZmod, Frac, Quotient, UnivariatePolynomial, MultivariatePolynomial
+export ZZ, ZZZ, ZI, QQ, ZZmod, Frac, Quotient, UnivariatePolynomial, MultivariatePolynomial
 export GaloisField, GF, FSeries, AlgebraicNumber, NumberField, NF
 
 export SpecialPowerSeries, PowerSeries, O, precision, absprecision
@@ -64,6 +71,7 @@ include("typevars.jl")
 include("promoteconvert.jl")
 include("generic.jl")
 include("zz.jl")
+include("zzflint.jl")
 include("qq.jl")
 include("zzmod.jl")
 include("quotient.jl")

src/algebraic.jl

@@ -1,6 +1,6 @@
 
 import Base: *, /, inv, +, -, sqrt, ^, literal_pow, iszero, zero, one, ==, isapprox, hash
-import Base: conj, real, imag, abs, copy, isreal
+import Base: conj, real, imag, abs, copy, isreal, cispi, cospi, sinpi, tanpi
 
 # construction
 basetype(::Type{<:AlgebraicNumber}) = QQ{BigInt}
@@ -160,14 +160,11 @@ function root(a::AlgebraicNumber, n::Integer)
     AlgebraicNumber(p, approx(a)^(1 // n))
 end
 
-# The imaginary constant a algebraic number
-IM() = AlgebraicNumber(monom(QQ{Int}[:x], 2) + 1, im)
-
 sqrt(a::AlgebraicNumber) = ^(a, 1 // 2)
 # cbrt(a::AlgebraicNumber) = ^(a, 1 // 3) # intentionally not defined - alike Complex.
 isreal(a::AlgebraicNumber) = isreal(approx(a))
 real(a::AlgebraicNumber) = isreal(a) ? a : (a + conj(a)) / 2
-imag(a::AlgebraicNumber) = isreal(a) ? zero(a) : (conj(a) - a) * IM() / 2
+imag(a::AlgebraicNumber) = isreal(a) ? zero(a) : (conj(a) - a) * AlgebraicNumber(im) / 2
 abs(a::AlgebraicNumber) = !isreal(a) ? sqrt(a * conj(a)) : real(approx(a)) >= 0 ? a : -a
 
 function *(a::T, b::T) where T<:AlgebraicNumber
@@ -427,8 +424,8 @@ function Base.sinpi(q::QQ)
     b = inv(a)
     s = (a - b) / 2
     fs = sinpi(big(value(q))) * im
-    as = AlgebraicNumber(s, fs) / IM()
-    as
+    as = AlgebraicNumber(s, fs)
+    AlgebraicNumber(_squaremulim(minimal_polynomial(as)), approx(as) / im)
 end
 function Base.cospi(q::QQ)
     d = cispi(q)
@@ -440,6 +437,29 @@ function Base.cospi(q::QQ)
     ac = AlgebraicNumber(c, fc)
     ac
 end
+function Base.tanpi(q::QQ)
+    d = cispi(q)
+    N = NF(d)
+    a = monom(N)
+    b = inv(a)
+    c = (a - b) / (a + b)
+    fc = tanpi(big(value(q))) * im
+    ac = AlgebraicNumber(c, fc)
+    AlgebraicNumber(_squaremulim(minimal_polynomial(ac)), approx(ac) / im)
+end
+
+# only for internal usage
+# assume minimal polynomial has only powers of x^2.
+# simulate effect of multiplication with `im` in minimal polynomial
+function _squaremulim(p::UnivariatePolynomial)
+    c = copy(p)
+    n = length(c.coeff)
+    for k = n-2:-4:1
+        c.coeff[k] = -c.coeff[k]
+    end
+    c
+end
+
 
 """
     rationalconst(expr)

src/generic.jl

@@ -48,6 +48,7 @@ Base.big(a::T) where T<:Union{QQ,ZZ} = big(T)(a)
 basetype(::T) where T<:Ring = basetype(T)
 basetype(::Type{T}) where T = Union{}
 basetype(::Type{Union{}}) = Union{}
+basetype(::Type{Rational{S}}) where S = S
 
 depth(::Type{Union{}}) = -1
 depth(::Type) = 0

src/intfactorization.jl

@@ -3,12 +3,12 @@
 """
 const MINPRIME = 9999
 
-function isirreducible(p::P; p0 = MINPRIME) where P<:UnivariatePolynomial{<:ZZ}
+function isirreducible(p::P; p0 = MINPRIME) where {T<:ZI, P<:UnivariatePolynomial{T}}
     (iszero(p) || isunit(p)) && return false
     deg(p) <= 1 && return true
     iszero(p[0]) && return false
     X = varname(P)
-    Z = ZZ{BigInt}[X]
+    Z = wide_type(T)[X]
     q = convert(Z, p)
     isone(pgcd(q, derive(q))) || return false
     zassenhaus_irr(q, p0)
@@ -20,11 +20,11 @@ function isirreducible(p::P; p0 = MINPRIME) where P<:UnivariatePolynomial{<:QQ}
     isirreducible(pp; p0)
 end
 
-function factor(p::P, a::Integer = 1; p0 = MINPRIME) where P<:UnivariatePolynomial{<:ZZ}
+function factor(p::P, a::Integer=1; p0 = MINPRIME) where {T<:ZI,P<:UnivariatePolynomial{T}}
     #println("factor($p)")
     X = varname(P)
     c = content(p)
-    Z = ZZ{BigInt}[X]
+    Z = wide_type(T)[X]
     q = Z(isone(c) ? copy(p) : p / c)
     x = monom(Z)
     q, e, k = stripzeroscompress(q)
@@ -145,8 +145,7 @@ function zassenhaus2(u::UnivariatePolynomial{<:ZZ{<:Integer}}, val::Val{BO}, p0)
 end
 
 # D = []
-
-function zassenhaus2(u::P, ::Val{BO}, p0) where {BO,P<:UnivariatePolynomial{ZZ{BigInt}}}
+function zassenhaus2(u::P, ::Val{BO}, p0) where {BO,T<:Union{ZZ{BigInt},ZZZ},P<:UnivariatePolynomial{T}}
     v, p = best_prime(u, p0)
     a = allgcdx(v)
     #println(" initial v/$p = "); display([v a])
@@ -266,7 +265,7 @@ function factor_exp(u::P, a::Integer, p0) where P<:UnivariatePolynomial
 end
 
 # check if p is coprime with leading coefficient
-function compatible_with(p::Integer, u::ZZ)
+function compatible_with(p::Integer, u::ZI)
     u = abs(value(u))
     gcd(p, u) == 1
 end
@@ -335,7 +334,7 @@ function combinefactors(
     u::Z,
     vv::AbstractVector{<:UnivariatePolynomial{Zp}},
     aa::AbstractVector,
-) where {Z<:UnivariatePolynomial{<:ZZ},Zp}
+) where {Z<:UnivariatePolynomial{<:ZI},Zp}
     res = Tuple{Z,Any,Any}[]
     un = LC(u)
     unp = Zp(un)
@@ -394,10 +393,10 @@ end
 function stripmod(
     ::Type{P},
     a::UnivariatePolynomial{<:ZZmod},
-) where {Z<:ZZ,P<:UnivariatePolynomial{Z}}
+) where {Z<:ZI,P<:UnivariatePolynomial{Z}}
     P(stripmod.(Z, a.coeff), ord(a))
 end
-stripmod(::Type{Z}, a::ZZmod) where Z<:ZZ = Z(value(a))
+stripmod(::Type{Z}, a::ZZmod) where Z<:ZI = Z(value(a))
 
 
 function subset_with_a(v, d, a)
@@ -421,8 +420,8 @@ If `u` has an integer factor polynom `v` with `deg(v) == m`,
 calculate array of bounds `b` with `abs(v[i]) <= b[i+1] for i = 0:m`.
 Algorithm see TAoCP 2.Ed 4.6.2 Exercise 20.
 """
-function coeffbounds(u::UnivariatePolynomial{ZZ{T},X}, m::Integer) where {T<:Integer,X}
-    W = widen(T)
+function coeffbounds(u::UnivariatePolynomial{T,X}, m::Integer) where {T<:ZI,X}
+    W = widen(basetype(T))
     n = deg(u)
     0 <= m <= n || throw(ArgumentError("required m ∈ [0,deg(u)] but $m ∉ [0,$n]"))
     accuracy = 100 # use fixed point decimal arithmetic with accuracy 0.01 for the norm
@@ -475,7 +474,7 @@ function hensel_lift(
     v::AbstractVector{Pq},
     a::AbstractVector{Pp},
 ) where {
-    Z<:ZZ,
+    Z<:ZI,
     Zq<:ZZmod,
     Zp<:ZZmod,
     Pq<:UnivariatePolynomial{Zq},
@@ -517,10 +516,10 @@ function downmod(::Type{Zp}, f::P, q::Integer) where {Zp,X,T,P<:UnivariatePolyno
 end
 
 
-function _liftmod(::Type{Z}, a::ZZmod) where {T,Z<:ZZ{T}}
-    Z(signed(T)(value(a)))
+function _liftmod(::Type{Z}, a::ZZmod) where {Z<:ZI}
+    Z(signed(basetype(Z))(value(a)))
 end
-function _liftmod(::Type{Z}, a::ZZ) where {X,T,Z<:ZZmod{X,T}}
+function _liftmod(::Type{Z}, a::ZI) where {X,T,Z<:ZZmod{X,T}}
     Z(value(a))
 end
 function _liftmod(::Type{Z}, a::ZZmod) where {X,T,Z<:ZZmod{X,T}}

src/promoteconvert.jl

@@ -1,18 +1,18 @@
 
 # promotions
 
-_xpromote_rule(::Type{T}, ::Type{S}) where {T<:Ring,S<:RingInt} = begin
+function promote_rule(::Type{T}, ::Type{S}) where {T<:Ring,S<:RingInt}
     depth(T) < depth(S) && return Base.Bottom
     B = basetype(T)
     (S <: T || S <: B) ? T : promote_type(B, S) == B ? T : Base.Bottom
 end
 _xpromote_rule(a::Type, b::Type) = promote_type(a, b)
-
+#=
 @generated function Base.promote_rule(a::Type{<:Ring}, b::Type{<:RingInt})
     bt = _xpromote_rule(a.parameters[1], b.parameters[1])
     :($bt)
 end
-
+=#
 function promote_rule(::Type{R}, ::Type{S}) where {R<:Ring,S<:Rational}
     promote_rule(R, promote_type(basetype(R), S))
 end

src/ringtypes.jl

@@ -93,6 +93,55 @@ struct ZZ{T<:Signed} <: Ring{ZZClass{T}}
     ZZ(val::T) where T<:Signed = ZZ{T}(val)
 end
 
+mutable struct ZZZ <: Ring{ZZClass{Integer}}
+    d::Int
+
+    function ZZZ(x::ZZZ)
+        z = new()
+        ccall((:fmpz_init_set, libflint), Nothing, (Ref{ZZZ}, Ref{ZZZ}), z, x)
+        finalizer(_fmpz_clear_fn, z)
+        return z
+    end
+
+    function ZZZ(x::Int)
+        z = new()
+        ccall((:fmpz_init_set_si, libflint), Nothing, (Ref{ZZZ}, Int), z, x)
+        finalizer(_fmpz_clear_fn, z)
+        return z
+    end
+
+    function ZZZ(x::UInt)
+        z = new()
+        ccall((:fmpz_init_set_ui, libflint), Nothing, (Ref{ZZZ}, UInt), z, x)
+        finalizer(_fmpz_clear_fn, z)
+        return z
+    end
+
+    function ZZZ(x::BigInt)
+        z = new()
+        ccall((:fmpz_init, libflint), Nothing, (Ref{ZZZ},), z)
+        ccall((:fmpz_set_mpz, libflint), Nothing, (Ref{ZZZ}, Ref{BigInt}), z, x)
+        finalizer(_fmpz_clear_fn, z)
+        return z
+    end
+
+    function ZZZ(x::Float64)
+        !isinteger(x) && throw(InexactError(:convert, ZZZ, x))
+        z = new()
+        ccall((:fmpz_init, libflint), Nothing, (Ref{ZZZ},), z)
+        ccall((:fmpz_set_d, libflint), Nothing, (Ref{ZZZ}, Cdouble), z, x)
+        finalizer(_fmpz_clear_fn, z)
+        return z
+    end
+
+end
+
+const ZI = Union{ZZ,ZZZ}
+
+function _fmpz_clear_fn(a::ZZZ)
+   ccall((:fmpz_clear, libflint), Nothing, (Ref{ZZZ},), a)
+end
+
 """
     ZZmod{m,S<:Integer}
 
@@ -111,7 +160,7 @@ The ring of fractions of `R`. The elements consist of pairs `num::R,den::R`.
 During creation the values may be canceled to achieve `gcd(num, den) == one(R)`.
 The special case of `R<:Integer` is handled by `QQ{R}`.
 """
-struct Frac{P<:Union{Polynomial,ZZ}} <: FractionRing{P,FractionClass{P}}
+struct Frac{P<:Union{Polynomial,ZI}} <: FractionRing{P,FractionClass{P}}
     num::P
     den::P
     Frac{P}(num::P, den::P, ::NCT) where P = new{P}(num, den)

src/univarpolynom.jl

@@ -99,6 +99,8 @@ function _promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{S}) where {X,R,
 end
 promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{S}) where {X,R,S<:Integer} =
     UnivariatePolynomial{promote_type(R, S),X}
+promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{S}) where {X,R,S<:ZI} =
+        UnivariatePolynomial{promote_type(R, S),X}
 promote_rule(::Type{UnivariatePolynomial{R,X}}, ::Type{S}) where {X,R,S<:Rational} =
     UnivariatePolynomial{promote_type(R, S),X}
 
@@ -826,7 +828,7 @@ end
 
 import Base: show
 
-issimple(::Union{ZZ,ZZmod,QQ,Number}) = true
+issimple(::Union{ZZ,ZZmod,QQ,Number,ZZZ}) = true
 issimple(::Quotient{<:UnivariatePolynomial{S,:γ}}) where S = true
 issimple(p::Polynomial) = ismonom(p)
 issimple(::Any) = false

src/zz.jl

@@ -23,8 +23,9 @@ function ZZ{T}(a::Union{QQ{T},Frac{ZZ{T}}}) where T
     a.den != 1 && throw(InexactError(:ZZ, ZZ{T}, a))
     ZZ(a.num)
 end
-#ZZ(a::Union{QQ{T},Frac{ZZ{T}}}) where T = ZZ{T}(a)
-ZZ{S}(a::Union{QQ{T},Frac{ZZ{T}}}) where {S,T} = ZZ(promote_type(S,T)(a))
+
+ZZ(a::Union{QQ{T},Frac{ZZ{T}}}) where T = ZZ{T}(a)
+ZZ{S}(a::Union{QQ{T},Frac{ZZ{T}}}) where {S,T} = ZZ(promote_type(S, T)(a))
 
 # promotion and conversion
 promote_rule(::Type{ZZ{T}}, ::Type{ZZ{S}}) where {S,T} = ZZ{promote_type(S, T)}
@@ -78,3 +79,5 @@ Base.show(io::IO, z::ZZ) = print(io, z.val)
 function (::Type{T})(a::ZZ) where T<:Integer
     T(value(a))
 end
+
+wide_type(::Type{<:ZZ}) = ZZ{BigInt}