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" diff --git a/src/IntegerMathUtils.jl b/src/IntegerMathUtils.jl index af54a42..b7b3e30 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) + return Int(abs(a) == 1) + 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 deleted file mode 100644 index 57c81f5..0000000 --- a/tests/runtests.jl +++ /dev/null @@ -1,64 +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 -