diff --git a/Project.toml b/Project.toml index 787172b..e32d702 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Mods" uuid = "7475f97c-0381-53b1-977b-4c60186c8d62" author = ["Edward Scheinerman "] -version = "2.0.3" +version = "2.1.0" [compat] julia = "1" diff --git a/README.md b/README.md index 94c338b..e0a0496 100644 --- a/README.md +++ b/README.md @@ -376,7 +376,7 @@ julia> rand(Mod{10},2,5) ### Rationals and Mods The result of `Mod{N}(a//b)` is exactly -`Mod{N}(numerator(a)) / Mod{n}(denominator(b))`. This may equal +`Mod{N}(numerator(a)) / Mod{N}(denominator(b))`. This may equal `Mod{N}(a)/Mod{N}(b)` if `a` and `b` are relatively prime to each other and to `N`. diff --git a/junk/old_GaussMods.jl b/junk/old_GaussMods.jl deleted file mode 100644 index 0f88d69..0000000 --- a/junk/old_GaussMods.jl +++ /dev/null @@ -1,55 +0,0 @@ -import Base: mod, real, imag, reim, conj, promote_rule - -export GaussMod, AbstractMod - -# mod support for Gaussian integers until officially adopted into Base -mod(z::Complex{<:Integer}, n::Integer) = Complex(mod(real(z), n), mod(imag(z), n)) - -""" -`GaussMod{N,T}` is an alias of `Mod{N,Complex{T}}`. -It is for computing Gaussian Modulus. -""" -const GaussMod{N,T} = Mod{N,Complex{T}} -GaussMod{N}(x::T) where {N,T<:Integer} = Mod{N,Complex{T}}(x) -GaussMod{N}(x::T) where {N,T<:Complex} = Mod{N,T}(x) -Mod{N}(x::Complex{Rational{T}}) where {N,T} = Mod{N}(real(x)) + Mod{N}(imag(x)) * im -Mod{N}(re::Integer, im::Integer) where {N} = Mod{N}(Complex(re, im)) - -reim(x::AbstractMod) = (real(x), imag(x)) -real(x::Mod{N}) where {N} = Mod{N}(real(x.val)) -imag(x::Mod{N}) where {N} = Mod{N}(imag(x.val)) -conj(x::Mod{N}) where {N} = Mod{N}(conj(x.val)) - -# ARITHMETIC -function (+)(x::GaussMod{N,T}, y::GaussMod{N,T}) where {N,T} - xx = widen(x.val) - yy = widen(y.val) - zz = mod(xx + yy, N) - return GaussMod{N,T}(zz) -end - -(-)(x::GaussMod{N}) where {N} = Mod{N}(-x.val) - -function (*)(x::GaussMod{N,T}, y::GaussMod{N,T}) where {N,T} - xx = widen(x.val) - yy = widen(y.val) - zz = mod(xx * yy, N) - return GaussMod{N,T}(zz) -end - -function inv(x::GaussMod{N}) where {N} - try - a, b = reim(x) - bot = inv(a * a + b * b) - aa = a * bot - bb = -b * bot - return aa + bb * im - catch - error("$x is not invertible") - end -end - -is_invertible(x::GaussMod) = is_invertible(real(x * x')) -#rand(::Type{GaussMod{N}}, args::Integer...) where {N} = rand(GaussMod{N,Int}, args...) - -rand(::Type{GaussMod{N}}, args::Integer...) where {N} = rand(Mod{N}, args...) .+ im*rand(Mod{N}, args...) diff --git a/junk/old_Mods.jl b/junk/old_Mods.jl deleted file mode 100644 index 2f936ba..0000000 --- a/junk/old_Mods.jl +++ /dev/null @@ -1,200 +0,0 @@ -module Mods - -import Base: (==), (+), (-), (*), (inv), (/), (//), (^), hash, show -import Base: rand, conj, iszero, rtoldefault, isapprox - -export Mod, modulus, value, AbstractMod -export is_invertible -export CRT - -abstract type AbstractMod <: Number end -__init__() = @warn "This version of Mods is a bit buggy. Don't upgrade to version 2.0.x yet. I'm working on it." - -""" -`Mod{N}(v)` creates a modular number in mod `N` with value `mod(v,N)`. -""" -struct Mod{N,T} <: AbstractMod - val::T -end - -# safe constructors (slower) -function Mod{N}(x::T) where {T<:Integer,N} - @assert N isa Integer && N > 1 "modulus must be an integer and at least 2" - v = Int(mod(x, N)) - Mod{Int(N),Int}(v) -end - -function Mod{N}(x::T) where {T<:Complex{<:Integer},N} - @assert N isa Integer && N > 1 "modulus must be an integer and at least 2" - S = Complex{typeof(N)} - v = mod(x, N) - Mod{Int(N),S}(v) -end - - - -# type casting -Mod{N,T}(x::Mod{N,T2}) where {T,N,T2} = Mod{N,T}(T(x.val)) -Mod{N,T}(x::Mod{N,T}) where {T,N} = x - -show(io::IO, z::Mod{N}) where {N} = print(io, "Mod{$N}($(value(z)))") -show(io::IO, ::MIME"text/plain", z::Mod{N}) where {N} = show(io, z) - -""" -`modulus(a::Mod)` returns the modulus of this `Mod` number. -``` -julia> a = Mod{13}(11); - -julia> modulus(a) -13 -``` -""" -modulus(a::Mod{N}) where {N} = N - -""" -`value(a::Mod)` returns the value of this `Mod` number. -``` -julia> a = Mod{13}(11); - -julia> value(a) -11 -``` -""" -value(a::Mod{N}) where {N} = a.val - -Base.abs(a::Mod{N,<:Real} where {N}) = abs(value(a)) - -function hash(x::Mod, h::UInt64 = UInt64(0)) - v = value(x) - m = modulus(x) - return hash(v, hash(m, h)) -end - -# Test for equality -iszero(x::Mod{N}) where {N} = iszero(x.val) -==(x::Mod, y::Mod) = modulus(x) == modulus(y) && value(x) == value(y) -# ==(x::Mod{N}, y::Mod{N}) where {N} = iszero(value(x - y)) - -# Apporximate equality -rtoldefault(::Type{Mod{N,T}}) where {N,T} = rtoldefault(T) -isapprox(x::Mod, y::Mod; kwargs...) = false -isapprox(x::Mod{N}, y::Mod{N}; kwargs...) where {N} = - isapprox(value(x), value(y); kwargs...) || isapprox(value(y), value(x); kwargs...) - -# Easy arithmetic -@inline function +(x::Mod{N}, y::Mod{N}) where {N} - t = widen(x.val) + widen(y.val) # add with added precision - return Mod{N}(mod(t, N)) -end - - -@inline function -(x::Mod{M}) where {M} - return Mod{M}(-x.val) -end - -# function -(x::Mod{M,T}) where {M,T<:Unsigned} -# return Mod{M,T}(M - value(x)) -# end - -@inline -(x::Mod, y::Mod) = x + (-y) - -@inline function *(x::Mod{N}, y::Mod{N}) where {N} - q = widemul(x.val, y.val) # multipy with added precision - return Mod{N}(q) # return with proper type -end - -# Division stuff -""" -`is_invertible(x::Mod)` determines if `x` is invertible. -""" -@inline function is_invertible(x::Mod{M})::Bool where {M} - return gcd(x.val, M) == 1 -end - - -""" -`inv(x::Mod)` gives the multiplicative inverse of `x`. -""" -@inline function inv(x::Mod{M}) where {M} - Mod{M}(_invmod(x.val, M)) -end -_invmod(x::Unsigned, m::Unsigned) = invmod(x, m) -# faster version of `Base.invmod`, only works for for signed types -@inline function _invmod(x::Signed, m::Signed) - (g, v, _) = gcdx(x, m) - if g != 1 - error("$x (mod $m) is not invertible") - end - return v -end - -function /(x::Mod{N,T}, y::Mod{N,T}) where {N,T} - return x * inv(y) -end - -(//)(x::Mod, y::Mod) = x / y -(//)(x::Number, y::Mod{N}) where {N} = x / y -(//)(x::Mod{N}, y::Number) where {N} = x / y - -Base.promote_rule(::Type{Mod{M,T1}}, ::Type{Mod{N,T2}}) where {M,N,T1,T2<:Number} = - error("can not promote types `Mod{$M,$T1}`` and `Mod{$N,$T2}`") -Base.promote_rule(::Type{Mod{M,T1}}, ::Type{Mod{M,T2}}) where {M,T1,T2<:Number} = - Mod{M,promote_type(T1, T2)} -Base.promote_rule(::Type{Mod{M,T1}}, ::Type{T2}) where {M,T1,T2<:Number} = - Mod{M,promote_type(T1, T2)} -Base.promote_rule(::Type{Mod{M,T1}}, ::Type{Rational{T2}}) where {M,T1,T2} = - Mod{M,promote_type(T1, T2)} - -# Operations with rational numbers -Mod{N}(k::Rational) where {N} = Mod{N}(numerator(k)) / Mod{N}(denominator(k)) -Mod{N,T}(k::Rational{T2}) where {N,T,T2} = Mod{N,T}(numerator(k)) / Mod{N,T}(denominator(k)) - -# Random -# rand(::Type{Mod{N}}, args::Integer...) where {N} = Mod{N}(rand(Int)) -rand(::Type{Mod{N}}, dims::Integer...) where {N} = Mod{N}.(rand(Int, dims...)) - -# Chinese remainder theorem functions -""" - `CRT([T=BigInt, ]m1, m2,...)` - -Chinese Remainder Theorem. - -``` -julia> CRT(Int, Mod{11}(4), Mod{14}(814)) -92 - -julia> 92%11 -4 - -julia> 92%14 -8 - -julia> CRT(Mod{9223372036854775783}(9223372036854775782), Mod{9223372036854775643}(9223372036854775642)) -85070591730234614113402964855534653468 -``` - -!!! note - - `CRT` uses `BigInt` by default to prevent potential integer overflow. - If you are confident that numbers do not overflow in your application, - please specify an optional type parameter as the first argument. -""" -function CRT(::Type{T}, remainders, primes) where {T} - length(remainders) == length(primes) || error("size mismatch") - isempty(remainders) && return zero(T) - primes = convert.(T, primes) - M = prod(primes) - Ms = M .÷ primes - ti = _invmod.(Ms, primes) - mod(sum(convert.(T, remainders) .* ti .* Ms), M) -end - -function CRT(::Type{T}, rs::Mod...) where {T} - CRT(T, value.(rs), modulus.(rs)) -end -CRT(rs::Mod...) = CRT(BigInt, rs...) - -include("GaussMods.jl") -include("iterate.jl") - -end # end of module Mods diff --git a/src/Mods.jl b/src/Mods.jl index 7dad260..deb982b 100644 --- a/src/Mods.jl +++ b/src/Mods.jl @@ -1,10 +1,6 @@ module Mods -__init__() = - @warn "This version of Mods is still under development. Don't upgrade to version 2.0.x yet. I'm working on it." - - import Base: (==), (+), (-), (*), (inv), (/), (//), (^) import Base: hash, show, iszero, isone, mod, abs,conj diff --git a/test/runtests.jl b/test/runtests.jl index bd95ca3..354dbd3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,12 @@ using Mods @test oneunit(Mod{17}) == Mod{17}(1) @test zero(Mod{17}) == 0 @test iszero(zero(Mod{17})) - + @test one(GaussMod{17}) == 1 + @test one(GaussMod{17}) == one(Mod{17}) + @test value(GaussMod{17}(1,2)) == 1 + 2im + @test Mod(Mod{17}(1)) == one(Mod{17}) + @test GaussMod(GaussMod{17}(2,3)) == 2 + 3im + @test Mod(GaussMod{17}(2,3)) == 2 + 3im end @@ -26,7 +31,7 @@ end @test a * b == Mod{p}(17) @test (a / b) * b == a - # @test (b // a) * (2 // 1) == b + @test (b // a) * (2 // 1) == b @test (a * 2) // 3 == (2a) * inv(Mod{p}(3)) @test is_invertible(a) @@ -35,36 +40,37 @@ end @test a^(p - 1) == 1 @test a^(-1) == inv(a) - # @test -Mod{13,Int}(typemin(Int)) == -Mod{13}(mod(typemin(Int), 13)) + @test -Mod{13}(typemin(Int)) == -Mod{13}(mod(typemin(Int), 13)) + @test Mod{13}(typemax(Int)) == typemax(Int) % 13 end -# @testset "GaussMod arithmetic" begin -# @test one(GaussMod{6}) == Mod{6}(1 + 0im) -# @test zero(GaussMod{6}) == Mod{6}(0 + 0im) -# @test GaussMod{6}(2 + 3im) == Mod{6}(2 + 3im) -# @test GaussMod{6,Int}(2 + 3im) == Mod{6}(2 + 3im) -# @test rand(GaussMod{6}) isa GaussMod -# @test eltype(rand(GaussMod{6}, 4, 4)) <: GaussMod -# p = 23 -# a = GaussMod{p}(3 - im) -# b = Mod{p}(5 + 5im) +@testset "GaussMod arithmetic" begin + @test one(GaussMod{6}) == Mod{6}(1 + 0im) + @test zero(GaussMod{6}) == Mod{6}(0 + 0im) + @test GaussMod{6}(2 + 3im) == Mod{6}(2 + 3im) + @test GaussMod{6}(2 + 3im) == Mod{6}(2 + 3im) + @test rand(GaussMod{6}) isa GaussMod + @test eltype(rand(GaussMod{6}, 4, 4)) <: GaussMod + p = 23 + a = GaussMod{p}(3 - im) + b = Mod{p}(5 + 5im) -# @test a + b == 8 + 4im -# @test a + Mod{p}(11) == Mod{p}(14, 22) -# @test -a == 20 + im -# @test a - b == Mod{p}(3 - im - 5 - 5im) + @test a + b == 8 + 4im + @test a + Mod{p}(11) == GaussMod{p}(14, 22) + @test -a == 20 + im + @test a - b == Mod{p}(3 - im - 5 - 5im) -# @test a * b == Mod{p}((3 - im) * (5 + 5im)) -# @test a / b == Mod{p}((3 - im) // (5 + 5im)) + @test a * b == Mod{p}((3 - im) * (5 + 5im)) + @test a / b == Mod{p}((3 - im) // (5 + 5im)) -# @test a^(p * p - 1) == 1 -# @test is_invertible(a) -# @test a * inv(a) == 1 + @test a^(p * p - 1) == 1 + @test is_invertible(a) + @test a * inv(a) == 1 -# @test a / (1 + im) == a / Mod{p}(1 + im) -# @test imag(a * a') == 0 -# end + @test a / (1 + im) == a / Mod{p}(1 + im) + @test imag(a * a') == 0 +end @testset "Large Modulus" begin @@ -94,33 +100,6 @@ end -# @testset "CRT" begin -# p = 23 -# q = 91 -# a = Mod{p}(17) -# b = Mod{q}(32) - -# @test CRT(Int32, [], []) === Int32(0) -# x = CRT(a, b) -# @test typeof(x) <: BigInt -# @test a == mod(x, p) -# @test b == mod(x, q) - -# c = Mod{101}(86) -# x = CRT(Int, a, b, c) -# @test typeof(x) <: Int - -# @test a == mod(x, p) -# @test b == mod(x, q) -# @test c == mod(x, 101) - -# @test CRT( -# BigInt, -# Mod{9223372036854775783}(9223372036854775782), -# Mod{9223372036854775643}(9223372036854775642), -# ) == 85070591730234614113402964855534653468 -# end - @testset "Matrices" begin M = zeros(Mod{11}, 3, 3) @@ -147,17 +126,6 @@ end end -# @testset "isapprox" begin -# @test Mod{7}(3) ≈ Mod{7}(3) -# @test Mod{7}(3) ≈ Mod{7}(10) -# @test Mod{7}(3) ≈ Mod{7}(10) atol = 1 -# @test Mod{7}(3) ≈ Mod{7}(11) atol = 1 -# @test !isapprox(Mod{7}(3), Mod{7}(11), atol = 0) -# @test !(Mod{7}(3) ≈ Mod{7}(11)) -# @test Mod{7}(3) ≈ Mod{7}(11) atol = 2 -# end - - @testset "Moduli types" begin a = Mod{17}(-1) b = Mod{17}(16)