From d8c67b6e756be1b3329774d3c949d0f72567e75f Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 2 Jun 2023 13:23:34 -0400 Subject: [PATCH 1/3] add kronecker function (useful for lots of number theory stuff) --- src/IntegerMathUtils.jl | 58 +++++++++++++++++++++++++++++++++++++---- tests/runtests.jl | 6 +++++ 2 files changed, 59 insertions(+), 5 deletions(-) diff --git a/src/IntegerMathUtils.jl b/src/IntegerMathUtils.jl index af54a42..3a0ae3a 100644 --- a/src/IntegerMathUtils.jl +++ b/src/IntegerMathUtils.jl @@ -1,5 +1,5 @@ module IntegerMathUtils -export iroot, ispower, rootrem, find_exponent, is_probably_prime +export iroot, ispower, rootrem, find_exponent, is_probably_prime, kronecker iroot(x::Integer, n::Integer) = iroot(x, Cint(n)) @@ -10,7 +10,7 @@ function iroot(x::BigInt, n::Cint) n <= 0 && throw(DomainError(n, "Exponent must be > 0")) x <= 0 && iseven(x) && throw(DomainError(n, "This is a math no-no")) ans = BigInt() - ccall((:__gmpz_root, :libgmp), Cint, (Ref{BigInt}, Ref{BigInt}, Cint), ans, x, n) + @ccall :libgmp.__gmpz_root(ans::Ref{BigInt}, x::Ref{BigInt}, n::Cint)::Cint ans end @@ -22,7 +22,7 @@ function rootrem(x::T, n::Integer) where {T<:Integer} x <= 0 && iseven(x) && throw(DomainError(n, "This is a math no-no")) root = BigInt() rem = BigInt() - ccall((:__gmpz_rootrem, :libgmp), Nothing,(Ref{BigInt}, Ref{BigInt}, Ref{BigInt}, Int), root, rem, x, n) + @ccall :libgmp.__gmpz_rootrem(root::Ref{BigInt}, rem::Ref{BigInt}, x::Ref{BigInt}, n::Int)::Nothing return (root, T(rem)) end @@ -30,7 +30,7 @@ end ispower(x::Integer) = ispower(big(x)) function ispower(x::BigInt) - return 0 != ccall((:__gmpz_perfect_power_p, :libgmp), Cint, (Ref{BigInt},), x) + return 0 != @ccall :libgmp.__gmpz_perfect_power_p(x::Ref{BigInt})::Cint end # TODO: Add more efficient implimentation for smaller types @@ -43,7 +43,55 @@ function find_exponent(x::Integer) end function is_probably_prime(x::Integer; reps=25) - return ccall((:__gmpz_probab_prime_p, :libgmp), Cint, (Ref{BigInt}, Cint), x, reps) != 0 + if ! x isa BigInt + x = BigInt(x) + end + return 0 != @ccall :libgmp.__gmpz_probab_prime_p(x::Ref{BigInt}, reps::Cint)::Cint +end + +function kronecker(a::BigInt, b::Clong) + return @ccall :libgmp.__gmpz_kronecker_si(a::Ref{BigInt}, b::Clong)::Cint +end +function kronecker(a::Clong, b::BigInt) + return @ccall :libgmp.__gmpz_si_kronecker(a::Clong, b::Ref{BigInt})::Cint +end +function kronecker(a, n) + @assert n != -n || n == 0 + @assert a != -a || a == 0 + t = 1 + if iszero(n) && abs(a) != 1 + return 0 + end + if n < 0 + n = abs(n) + if a < 0 + t = -t + end + end + trail = trailing_zeros(n) + if trail > 0 + n >>= trail + if iseven(a) + return 0 + elseif isodd(trail) && a&7 in (3,5) + t = -t + end + end + a = mod(a, n) + while a != 0 + while iseven(a) + a = a >> 1 + if n&7 in (3, 5) + t = -t + end + end + a, n = n, a + if a&3 == n&3 == 3 + t = -t + end + a = mod(a, n) + end + return n == 1 ? t : 0 end end diff --git a/tests/runtests.jl b/tests/runtests.jl index 57c81f5..e692d85 100644 --- a/tests/runtests.jl +++ b/tests/runtests.jl @@ -62,3 +62,9 @@ end @test is_probably_prime(big(2)^128-1) == false end +@testset "kroneker" begin + @test kronecker(4, 5) == kronecker(4, -5) == 1 + @test kronecker(1, 0) == kronecker(-1, 0) == 1 + @test kronecker(-4, -5) == -1 + @test kronecker(4, 6) == kronecker(-4, 0) == 0 +end From 4ceadff9a378a6d4982f3b390cb210fc881c654d Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 2 Jun 2023 13:25:31 -0400 Subject: [PATCH 2/3] fix tests --- src/IntegerMathUtils.jl | 6 ++-- tests/runtests.jl | 70 ----------------------------------------- 2 files changed, 3 insertions(+), 73 deletions(-) delete mode 100644 tests/runtests.jl diff --git a/src/IntegerMathUtils.jl b/src/IntegerMathUtils.jl index 3a0ae3a..b7b3e30 100644 --- a/src/IntegerMathUtils.jl +++ b/src/IntegerMathUtils.jl @@ -43,7 +43,7 @@ function find_exponent(x::Integer) end function is_probably_prime(x::Integer; reps=25) - if ! x isa BigInt + if !(x isa BigInt) x = BigInt(x) end return 0 != @ccall :libgmp.__gmpz_probab_prime_p(x::Ref{BigInt}, reps::Cint)::Cint @@ -59,8 +59,8 @@ function kronecker(a, n) @assert n != -n || n == 0 @assert a != -a || a == 0 t = 1 - if iszero(n) && abs(a) != 1 - return 0 + if iszero(n) + return Int(abs(a) == 1) end if n < 0 n = abs(n) diff --git a/tests/runtests.jl b/tests/runtests.jl deleted file mode 100644 index e692d85..0000000 --- a/tests/runtests.jl +++ /dev/null @@ -1,70 +0,0 @@ -using Test - -@testset "iroot" begin - for T in (Int32, Int64, BigInt) - @test iroot(T(100), 2) == T(10) - @test iroot(T(101), 2) == T(10) - @test iroot(T(99), 2) == T(9) - @test iroot(T(1000), 3) == T(10) - @test iroot(T(1001), 3) == T(10) - @test iroot(T(999), 3) == T(9) - @test iroot(T(10000), 4) == T(10) - @test iroot(T(10001), 4) == T(10) - @test iroot(T(9999), 4) == T(9) - end - @test iroot(big(23)^50, 50) == big(23) - @test iroot(big(23)^50 + 1, 50) == big(23) - @test iroot(big(23)^50 - 1, 50) == big(22) -end - -@testset "ispower" begin - for T in (Int32, Int64, BigInt) - @test ispower(T(100)) == true - @test ispower(T(1)) == true - @test ispower(T(0)) == true - @test ispower(T(12)) == false - @test ispower(T(2^30)) == true - @test ispower(T(5)^10) == true - @test ispower(T(2^30)+1) == false - @test ispower(T(5)^10+1) == false - end - @test ispower(big(5)^40) == true - @test ispower(big(5)^40 + 1) == false - @test ispower(6 * big(5)^40) == false -end - -@testset "find_exponent" begin - for T in (Int32, Int64, BigInt) - @test find_exponent(T(100)) === 2 - @test find_exponent(T(1)) === 1 - @test find_exponent(T(0)) === 1 - @test find_exponent(T(12)) === 1 - @test find_exponent(T(2^30)) === 30 - @test find_exponent(T(5)^10) === 10 - @test find_exponent(T(2^30)+1) === 1 - @test find_exponent(T(5)^10+1) === 1 - end - @test find_exponent(big(5)^40) === 40 - @test find_exponent(big(5)^40 + 1) === 1 - @test find_exponent(6*big(5)^40 ) === 1 -end - -@testset "is_probably_prime" begin - for T in (Int32, Int64, BigInt) - @test is_probably_prime(T(2)^7-1) == true - @test is_probably_prime(T(2)^13-1) == true - @test is_probably_prime(T(2)^19-1) == true - @test is_probably_prime(T(2)^27-1) == false - @test is_probably_prime(T(2)^23-1) == false - @test is_probably_prime(T(2)^30-1) == false - end - @test is_probably_prime(big(2)^127-1) == true - @test is_probably_prime(big(2)^128-1) == false -end - -@testset "kroneker" begin - @test kronecker(4, 5) == kronecker(4, -5) == 1 - @test kronecker(1, 0) == kronecker(-1, 0) == 1 - @test kronecker(-4, -5) == -1 - @test kronecker(4, 6) == kronecker(-4, 0) == 0 -end From e3b23f2fdc6a3fbb6af8214a4066bfe8d14a0b3a Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 2 Jun 2023 14:01:13 -0400 Subject: [PATCH 3/3] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7379318..210aeca 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "IntegerMathUtils" uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235" -version = "0.1.1" +version = "0.1.2" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"