Widening depend on modulus

Author: scheinerman

Project.toml

@@ -1,7 +1,7 @@
 name = "Mods"
 uuid = "7475f97c-0381-53b1-977b-4c60186c8d62"
 author = ["Edward Scheinerman <[email protected]>"]
-version = "2.2.2"
+version = "2.2.3"
 
 [compat]
 julia = "1"

README.md

@@ -34,11 +34,25 @@ likely to give incorrect results, version 1 of this module was buggy.
 
 Users who require *smaller* integer (e.g., `Int8`) types should use the latest version 1 of `Mods`.
 
-> **NEW!** For `Int8` moduli (between 2 and 127) see the [MiniMods](https://github.com/scheinerman/MiniMods.jl) module.
+> **NEW!** For small moduli (between 2 and 255) see the [MiniMods](https://github.com/scheinerman/MiniMods.jl) module.
 
 In addition, some functionality has been moved to the `extras` folder. 
 See the `README` there. 
 
+### New in 2.2.3
+
+The values of `Mod` number is held in 64-bit integers. If the moduli and values are 
+large enough, integer arithmetic might overlow and yield incorrect results. To deal
+with this, integer values are expanded to 128 bits in order to ensure correctness,
+and then reduced by the modulus.
+
+In prior versions, we always expanded values to 128 bits before arithmetic. 
+
+Starting in version 2.2.3, expansion to 128 bits only happens for moduli above 
+`typemax(Int32)` which equals `2^31 - 1 = 2,147,483,647`. This results in a 
+roughly 4 or 5 times speed compared to prior versions.
+
+
 
 ## Quick Overview
 This module supports modular values and arithmetic. The moduli are integers (at least 2)

src/arithmetic.jl

@@ -1,14 +1,31 @@
-(+)(a::Mod{N}, b::Mod{N}) where {N} = Mod{N}(widen(value(a)) + widen(value(b)))
-(+)(a::GaussMod{N}, b::GaussMod{N}) where {N} = Mod{N}(widen(value(a)) + widen(value(b)))
+# (+)(a::Mod{N}, b::Mod{N}) where {N} = Mod{N}(widen(a.val) + widen(b.val))
+# (+)(a::GaussMod{N}, b::GaussMod{N}) where {N} = Mod{N}(widen(a.val) + widen(b.val))
 
-(-)(a::Mod{N}) where {N} = Mod{N}(-value(a))
-(-)(a::GaussMod{N}) where {N} = GaussMod{N}(-value(a))
+function (+)(a::Mod{N}, b::Mod{N}) where {N}
+    N <= typemax(Int32) ? Mod{N}(a.val + b.val) : Mod{N}(widen(a.val) + widen(b.val))
+end
+function (+)(a::GaussMod{N}, b::GaussMod{N}) where {N}
+    N <= typemax(Int32) ? Mod{N}((a.val) + (b.val)) :
+    GaussMod{N}(widen(a.val) + widen(b.val))
+end
+
+
+
+(-)(a::Mod{N}) where {N} = Mod{N}(N - a.val)
+(-)(a::GaussMod{N}) where {N} = GaussMod{N}(N - a.val)
 
 (-)(a::AbstractMod, b::AbstractMod) = a + (-b)
 
-(*)(a::Mod{N}, b::Mod{N}) where {N} = Mod{N}(widen(value(a)) * widen(value(b)))
-(*)(a::GaussMod{N}, b::GaussMod{N}) where {N} =
-    GaussMod{N}(widen(value(a)) * widen(value(b)))
+# (*)(a::Mod{N}, b::Mod{N}) where {N} = Mod{N}(widemul(a.val, b.val))
+# (*)(a::GaussMod{N}, b::GaussMod{N}) where {N} = GaussMod{N}(widemul(a.val, b.val))
+
+
+function (*)(a::Mod{N}, b::Mod{N}) where {N}
+    N <= typemax(Int32) ? Mod{N}(a.val * b.val) : Mod{N}(widemul(a.val, b.val))
+end
+function (*)(a::GaussMod{N}, b::GaussMod{N}) where {N}
+    N <= typemax(Int32) ? GaussMod{N}(a.val * b.val) : GaussMod{N}(widemul(a.val, b.val))
+end
 
 """
     is_invertible(a::AbstractMod)::Bool

test/speed_test.jl

@@ -0,0 +1,71 @@
+using Mods, BenchmarkTools, LinearAlgebra, LinearAlgebraX
+
+# m = typemax(Int32) - 18   # this is a big prime but < 2^31 - 1
+# d = 500
+
+function speed_test(m, d)
+
+    A = rand(Mod{m}, d, d)
+
+    @info "modulus = $m, matrix size = $d"
+
+    @info "Determinant"
+    @btime det(A)
+
+    @info "Rank"
+    @btime rankx(A)
+
+    @info "Inverse"
+    @btime inv(A)
+
+    @info "Multiplication"
+    @btime A * A * A
+
+    AA = value.(A)
+    @show A^3 == AA^3   # this should be false if big enough
+    nothing
+end
+
+## RESULT WITH NO CHANGE IN CODE
+# [ Info: modulus = 2147483629, matrix size = 500
+# [ Info: Determinant
+#   1.260 s (4 allocations: 1.91 MiB)
+# [ Info: Rank
+#   1.685 s (256007 allocations: 1001.13 MiB)
+# [ Info: Inverse
+#   6.202 s (8 allocations: 5.73 MiB)
+# [ Info: Multiplication
+#   5.557 s (12 allocations: 3.87 MiB)
+
+## RESULTS WITH REMOVING widemul and widen
+# [ Info: modulus = 2147483629, matrix size = 500
+# [ Info: Determinant
+#   156.085 ms (4 allocations: 1.91 MiB)
+# [ Info: Rank
+#   329.689 ms (256007 allocations: 1001.13 MiB)
+# [ Info: Inverse
+#   1.444 s (8 allocations: 5.73 MiB)
+# [ Info: Multiplication
+#   1.215 s (12 allocations: 3.87 MiB)
+
+## RESULTS WITH NEW CONTINGENT CODE 
+# [ Info: modulus = 2147483629, matrix size = 500  <-- mod is < 2^31
+# [ Info: Determinant
+#   11.625 μs (3 allocations: 3.48 KiB)
+# [ Info: Rank
+#   30.292 μs (594 allocations: 103.02 KiB)
+# [ Info: Inverse
+#   128.041 μs (5 allocations: 10.19 KiB)
+# [ Info: Multiplication
+#   61.791 μs (8 allocations: 47.53 KiB)
+
+# [ Info: modulus = 8589934609, matrix size = 500  <-- mod is > 2^32
+# [ Info: Determinant
+#   71.791 μs (3 allocations: 3.48 KiB)
+# [ Info: Rank
+#   125.250 μs (594 allocations: 103.02 KiB)
+# [ Info: Inverse
+#   409.250 μs (5 allocations: 10.19 KiB)
+# [ Info: Multiplication
+#   359.625 μs (8 allocations: 47.53 KiB)
+# A ^ 3 == AA ^ 3 = false