Optimize IsPrimePowerInt (speed up 10x)

Before:
gap> Number([-100000..100000], IsPrimePowerInt);; time;
1753

After:
gap> Number([-100000..100000], IsPrimePowerInt);; time;
168

hpcgap/lib/integer.gi

@@ -923,8 +923,51 @@ InstallGlobalFunction( IsOddInt, n -> n mod 2 = 1 );
 ##
 #F IsPrimePowerInt( <n> ) . . . . . . . . . . . test for a power of a prime
 ##
-InstallGlobalFunction( IsPrimePowerInt,
- n -> IsPrimeInt( SmallestRootInt( n ) ) );
+InstallGlobalFunction( IsPrimePowerInt, function(n)
+ local k, r, s, p, l, q, i;
+
+ # check the argument
+ if n > 1 then k := 2; s := 1;
+ elif n < -1 then k := 3; s := -1; n := -n;
+ else return false;
+ fi;
+
+ # exclude small divisors, and thereby large exponents
+ for p in Primes do
+ if p*p > n then return true; fi; # n is prime
+ r := PVALUATION_INT(n, p);
+ if r > 0 then
+ if s = -1 and IsEvenInt(r) then return false; fi;
+ return n = p^r;
+ fi;
+ od;
+ l := LogInt( n, p );
+
+ # loop over the possible prime divisors of exponents
+ # use Fermat's little theorem to cast out impossible ones:
+ # for suppose we had r such that n = r^k. Then by Fermat,
+ # n^((q-1)/k) = r^(q-1) is congruent 0 or 1 mod q
+ i := Position(Primes, k);
+ while k <= l do
+ q := 2*k+1; while not IsPrimeInt(q) do q := q+2*k; od;
+ if PowerModInt( n, (q-1)/k, q ) <= 1 then
+ r := RootInt( n, k );
+ if r ^ k = n then
+ n := r;
+ l := QuoInt( l, k );
+ continue;
+ fi;
+ fi;
+ if i <> fail and i < Length(Primes) then
+ i := i + 1;
+ k := Primes[i];
+ else
+ k := NextPrimeInt( k );
+ fi;
+ od;
+
+ return IsPrimeInt(n);
+end);
 
 
 #############################################################################

lib/integer.gi

@@ -910,8 +910,51 @@ InstallGlobalFunction( IsOddInt, n -> n mod 2 = 1 );
 ##
 #F IsPrimePowerInt( <n> ) . . . . . . . . . . . test for a power of a prime
 ##
-InstallGlobalFunction( IsPrimePowerInt,
- n -> IsPrimeInt( SmallestRootInt( n ) ) );
+InstallGlobalFunction( IsPrimePowerInt, function(n)
+ local k, r, s, p, l, q, i;
+
+ # check the argument
+ if n > 1 then k := 2; s := 1;
+ elif n < -1 then k := 3; s := -1; n := -n;
+ else return false;
+ fi;
+
+ # exclude small divisors, and thereby large exponents
+ for p in Primes do
+ if p*p > n then return true; fi; # n is prime
+ r := PVALUATION_INT(n, p);
+ if r > 0 then
+ if s = -1 and IsEvenInt(r) then return false; fi;
+ return n = p^r;
+ fi;
+ od;
+ l := LogInt( n, p );
+
+ # loop over the possible prime divisors of exponents
+ # use Fermat's little theorem to cast out impossible ones:
+ # for suppose we had r such that n = r^k. Then by Fermat,
+ # n^((q-1)/k) = r^(q-1) is congruent 0 or 1 mod q
+ i := Position(Primes, k);
+ while k <= l do
+ q := 2*k+1; while not IsPrimeInt(q) do q := q+2*k; od;
+ if PowerModInt( n, (q-1)/k, q ) <= 1 then
+ r := RootInt( n, k );
+ if r ^ k = n then
+ n := r;
+ l := QuoInt( l, k );
+ continue;
+ fi;
+ fi;
+ if i <> fail and i < Length(Primes) then
+ i := i + 1;
+ k := Primes[i];
+ else
+ k := NextPrimeInt( k );
+ fi;
+ od;
+
+ return IsPrimeInt(n);
+end);
 
 
 #############################################################################

tst/testinstall/intarith.tst

@@ -1414,6 +1414,25 @@ gap> List(data, SmallestRootInt);
 gap> List([-2^101,-2^100,2^100,2^101], SmallestRootInt);
 [ -2, -16, 2, 2 ]
 
+#
+# test IsPrimePowerInt
+#
+gap> P:=2^255-19;; # big prime
+gap> Filtered([-10..10], IsPrimePowerInt);
+[ -8, -7, -5, -3, -2, 2, 3, 4, 5, 7, 8, 9 ]
+gap> SetX(Primes, [1..30], {p,k}->IsPrimePowerInt(p^k));
+[ true ]
+gap> SetX(Primes, [1..10], {p,k}->IsPrimePowerInt(P*p^k));
+[ false ]
+gap> SetX(Primes, [1..10], {p,k}->IsPrimePowerInt((P*p)^k));
+[ false ]
+gap> ForAll([1..30], k->IsPrimePowerInt(P^k));
+true
+gap> ForAll([1..10], k->not IsPrimePowerInt(1009*P^k));
+true
+gap> ForAll([1..30], k->not IsPrimePowerInt((1009*P)^k));
+true
+
 #
 gap> STOP_TEST( "intarith.tst", 1);