Add iroot and check parameter to root for Integer. (#939)

* Add iroot and check parameter to root for Integer.

* Switch to default check=true and fix bug.

* Add ispower and ispower_with_root for BigInt (untested).

* Implement ispower and ispower_with_root for integer types.

Author: wbhart

docs/src/integer.md

@@ -147,6 +147,16 @@ AbstractAlgebra.sqrt(a::BigInt)
 issquare(a::BigInt)
 ```
 
+```@docs
+root(a::BigInt)
+iroot(a::BigInt)
+```
+
+```@docs
+ispower(a::BigInt)
+ispower_with_root(a::BigInt)
+```
+
 ```@docs
 AbstractAlgebra.exp(a::BigInt)
 ```

src/AbstractAlgebra.jl

@@ -439,7 +439,7 @@ import .Generic: abs_series, abs_series_type,
                  isone, isreverse, isrimhook, issubmodule,
                  isunit,
                  laurent_ring, lcm, leading_coefficient, leading_monomial,
-		 leading_term, length,
+                 leading_term, length,
                  leglength, main_variable,
                  main_variable_extract, main_variable_insert,
                  map1, map2, map_from_func,
@@ -451,25 +451,25 @@ import .Generic: abs_series, abs_series_type,
                  monomial_iszero, monomial_set!,
                  MPolyBuildCtx, mullow_karatsuba,
                  ngens, norm, normalise,
-		 num_coeff, one,
+                 num_coeff, one,
                  order, ordering, parity, partitionseq, Perm, perm,
                  permtype, @perm_str, polcoeff, poly, poly_ring,
                  precision, preimage, preimage_map,
-		 prime, push_term!,
+                 prime, push_term!,
                  rand_ordering, reduce!,
                  rels, rel_series, rel_series_type,
-		 rescale!, retraction_map, reverse,
+                 rescale!, retraction_map, reverse,
                  right_kernel, rowlength, section_map, setcoeff!,
                  set_exponent_vector!, set_limit!,
                  setpermstyle, size,
                  sort_terms!, summands,
                  supermodule, term, terms, total_degree,
                  to_univariate, trailing_coefficient,
-		 truncate, upscale, var_index,
+                 truncate, upscale, var_index,
                  zero,
        # Moved from Hecke into Misc
-		 Loc, Localization, LocElem,
-		 roots, sturm_sequence
+                 Loc, Localization, LocElem,
+                 roots, sturm_sequence
 
 # Do not export inv, div, divrem, exp, log, sqrt, numerator and denominator as we define our own
 export abs_series, abs_series_type,
@@ -499,7 +499,7 @@ export abs_series, abs_series_type,
                  issubmodule, issymmetric,
                  isterm_recursive, isunit, iszero,
                  lcm, leading_coefficient, leading_monomial, leading_term,
-		 length,
+                 length,
                  main_variable, main_variable_extract, main_variable_insert,
                  mat, matrix_repr, max_fields, mod,
                  monomial, monomial!, monomials,
@@ -508,18 +508,18 @@ export abs_series, abs_series_type,
                  mul_ks, mul_red!, mullow_karatsuba, mulmod,
                  needs_parentheses, newton_to_monomial!, ngens,
                  normalise, nullspace, num_coeff,
-		 one, order, ordering,
+                 one, order, ordering,
                  parent_type, parity, partitionseq, Perm, perm, permtype,
                  @perm_str, polcoeff, polynomial, poly,
-		 poly_ring, pow_multinomial,
+                 poly_ring, pow_multinomial,
                  ppio, precision, preimage, prime,
-		 push_term!, rank,
+                 push_term!, rank,
                  rand_ordering, reduce!,
                  renormalize!, rel_series, rel_series_type, rels,
-		 resultant, resultant_ducos,
+                 resultant, resultant_ducos,
                  resultant_euclidean, resultant_subresultant,
                  resultant_sylvester, resx, reverse, rowlength,
-		 setcoeff!, set_exponent_vector!,
+                 setcoeff!, set_exponent_vector!,
                  setpermstyle,
                  size, sort_terms!, subst, summands, supermodule,
                  sylvester_matrix, term, terms, to_univariate,
@@ -528,8 +528,8 @@ export abs_series, abs_series_type,
                  MatrixElem, PolynomialElem,
        # Moved from Hecke into Misc
                  divexact_low, divhigh,
-		 ismonic, Loc, Localization, LocElem, mulhigh_n,
-		 PolyCoeffs, roots, sturm_sequence
+                 ismonic, Loc, Localization, LocElem, mulhigh_n,
+                 PolyCoeffs, roots, sturm_sequence
 
 # TODO remove these two once removed from dependent packages (Hecke)
 export displayed_with_minus_in_front, show_minus_one

src/julia/Integer.jl

@@ -4,6 +4,8 @@
 #
 ###############################################################################
 
+export iroot, ispower, ispower_with_root, root
+
 ###############################################################################
 #
 #   Data type and parent object methods
@@ -177,6 +179,159 @@ function issquare(a::T) where T <: Integer
    return a == s*s
 end
 
+###############################################################################
+#
+#   Root
+#
+###############################################################################
+
+function root(a::BigInt, n::Int; check::Bool=true)
+    a < 0 && iseven(n) && throw(DomainError((a, n),
+                      "Argument `a` must be positive if exponent `n` is even"))
+    n <= 0 && throw(DomainError(n, "Exponent must be positive"))
+    z = BigInt()
+    exact = Bool(ccall((:__gmpz_root, :libgmp), Cint,
+                  (Ref{BigInt}, Ref{BigInt}, Cint), z, a, n))
+    check && !exact && error("Not a perfect n-th power (n = $n)")
+    return z
+end
+
+@doc Markdown.doc"""
+    root(a::T, n::Int; check::Bool=true) where T <: Integer
+
+Return the $n$-th root of $a$. If `check=true` the function will test if the
+input was a perfect $n$-th power, otherwise an exception will be raised. We
+require $n > 0$ and also $a \geq 0$ if $n$ is even.
+"""
+function root(a::T, n::Int; check::Bool=true) where T <: Integer
+   if n == 2
+      a < 0 && throw(DomainError((a, n),
+                      "Argument `a` must be positive if exponent `n` is even"))
+      s = isqrt(a)
+      exact = true
+      if check
+         r = a - s*s
+         exact = r == 0
+         !exact && error("Not a perfect n-th power (n = $n)")
+      end
+      return s
+   else
+      return T(root(BigInt(a), n; check=check))
+   end
+end
+
+moduli3 = [7, 8, 13]
+residues3 = [[0, 1, 6], [0, 1, 3, 5, 7], [0, 1, 5, 8, 12]]
+
+moduli5 = [8, 11, 31]
+residues5 = [[0, 1, 3, 5, 7], [0, 1, 10], [0, 1, 5, 6, 25, 26, 30]]
+
+moduli7 = [8, 29, 43]
+residues7 = [[0, 1, 3, 5, 7], [0, 1, 12, 17, 28], [0, 1, 6, 7, 36, 37, 42]]
+
+function ispower_moduli(a::Integer, n::Int)
+   if mod(n, 3) == 0
+      for i = 1:length(moduli3)
+         if !(mod(a, moduli3[i]) in residues3[i])
+            return false
+         end
+      end
+   elseif (n % 5) == 0
+      for i = 1:length(moduli5)
+         if !(mod(a, moduli5[i]) in residues5[i])
+            return false
+         end
+      end
+   elseif (n % 3) == 0
+      for i = 1:length(moduli7)
+         if !(mod(a, moduli7[i]) in residues7[i])
+            return false
+         end
+      end
+   elseif isodd(n)
+      if !(mod(a, moduli5[1]) in residues5[1])
+            return false
+      end
+   end
+   return true
+end
+
+@doc Markdown.doc"""
+    ispower_with_root(a::T, n::Int) where T <: Integer
+
+Return `true, q` if $a$ is a perfect $n$-th power with $a = q^n$. Otherwise
+return `false, 0`. We require $n > 0$.
+"""
+function ispower_with_root(a::T, n::Int) where T <: Integer
+   n <= 0 && throw(DomainError(n, "exponent n must be positive"))
+   if n == 1 || a == 0 || a == 1
+      return (true, a)
+   elseif a == -1
+      return isodd(n) ? (true, a) : (false, zero(T))
+   elseif mod(n, 2) == 0 && a < 0
+      return false, zero(BigInt)
+   elseif !ispower_moduli(a, n)
+      return (false, zero(BigInt))
+   end
+      
+   q = BigInt()
+   r = BigInt()
+   ccall((:__gmpz_rootrem, :libgmp), Nothing,
+                     (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}, Int), q, r, a, n)
+   return iszero(r) ? (true, T(q)) : (false, zero(T))
+end
+
+@doc Markdown.doc"""
+    ispower(a::T, n::Int) where T <: Integer
+
+Return `true` if $a$ is a perfect $n$-th power, i.e. if there is an integer $b$
+such that $a = b^n$. We require $n > 0$.
+"""
+function ispower(a::T, n::Int) where T <: Integer
+   n <= 0 && throw(DomainError(n, "n is not positive"))
+   if n == 1 || a == 0 || a == 1
+      return true
+   elseif a == -1
+      return isodd(n) ? true : false
+   elseif mod(n, 2) == 0 && a < 0
+      return false
+   elseif !ispower_moduli(a, n)
+      return false
+   end
+   
+   q = BigInt()
+   r = BigInt()
+   ccall((:__gmpz_rootrem, :libgmp), Nothing,
+                     (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}, Int), q, r, BigInt(a), n)
+
+   return r == 0
+end
+
+function iroot(a::BigInt, n::Int)
+    a < 0 && iseven(n) && throw(DomainError((a, n),
+                      "Argument `a` must be positive if exponent `n` is even"))
+    n <= 0 && throw(DomainError(n, "Exponent must be positive"))
+    z = BigInt()
+    ccall((:__gmpz_root, :libgmp), Cint,
+                  (Ref{BigInt}, Ref{BigInt}, Cint), z, a, n)
+    return z
+end
+
+@doc Markdown.doc"""
+    iroot(a::T, n::Int) where T <: Integer
+
+Return the truncated integer part of the $n$-th root of $a$ (round towards
+zero). We require $n > 0$ and also $a \geq 0$ if $n$ is even.
+"""
+function iroot(a::T, n::Int) where T <: Integer
+   if n == 2
+       a < 0 && throw(DomainError((a, n),
+                      "Argument `a` must be positive if exponent `n` is even"))
+       return isqrt(a)
+   end
+   return T(iroot(BigInt(a), n))
+end
+
 ###############################################################################
 #
 #   Exponential

test/julia/Integers-test.jl

@@ -133,6 +133,73 @@ end
    end
 end
 
+@testset "Julia.Integers.root" begin
+   @test root(BigInt(1000), 3) == 10
+   @test root(-BigInt(27), 3) == -3
+   @test root(BigInt(27), 3; check=true) == 3
+   @test root(BigInt(16), 2; check=true) == 4
+
+   @test_throws DomainError root(-BigInt(1000), 4)
+   @test_throws DomainError root(BigInt(1000), -3)
+   @test_throws DomainError root(BigInt(-16), 2)
+
+   @test_throws ErrorException root(BigInt(1100), 3)
+   @test_throws ErrorException root(-BigInt(40), 3)
+
+   @test iroot(BigInt(1000), 3) == 10
+   @test iroot(BigInt(1100), 3) == 10
+   @test iroot(-BigInt(40), 3) == -3
+   @test iroot(BigInt(17), 2) == 4
+
+   @test_throws DomainError iroot(-BigInt(1000), 4)
+   @test_throws DomainError iroot(BigInt(1000), -3)
+   @test_throws DomainError iroot(-BigInt(16), 2)
+
+   @test root(1000, 3) == 10
+   @test root(-27, 3) == -3
+   @test root(27, 3) == 3
+   @test root(16, 2) == 4
+
+   @test_throws DomainError root(-1000, 4)
+   @test_throws DomainError root(1000, -3)
+   @test_throws DomainError root(-16, 2)
+
+   @test_throws ErrorException root(1100, 3)
+   @test_throws ErrorException root(-40, 3)
+
+   @test iroot(1000, 3) == 10
+   @test iroot(1100, 3) == 10
+   @test iroot(-40, 3) == -3
+   @test iroot(17, 2) == 4
+
+   @test_throws DomainError iroot(-1000, 4)
+   @test_throws DomainError iroot(1000, -3)
+   @test_throws DomainError iroot(-16, 2)
+
+   for T in [BigInt, Int]
+      for i = 1:1000
+         n = rand(1:20)
+         a = BigInt(rand(-1000:1000))
+         if iseven(n)
+            a = abs(a)
+         end
+         p = a^n
+         if T == BigInt || ndigits(p; base=2) < ndigits(typemax(T); base=2)
+            p = T(p)
+
+            @test ispower(p, n)
+
+            flag, q = ispower_with_root(p, n)
+
+            @test flag && q == a
+         end
+      end
+
+      @test_throws DomainError ispower(T(5), -1)
+      @test_throws DomainError ispower_with_root(T(5), 0)
+   end
+end
+
 @testset "Julia.Integers.exp" begin
    @test AbstractAlgebra.exp(0) == 1
    @test_throws DomainError AbstractAlgebra.exp(1)