gmpints: add PVALUATION_INT

It uses mpz_remove to compute the p-valuation. In fact, it computes
arbitrary g-valuations, i.e. the second argument can be an arbitrary
nonzero integer.

We also adapt PValuation to us PVALUATION_INT. For small inputs, there
is little difference, but for larger ones we get a noticeable speedup.

To illustrate this, here are some examples, where
  ps:=[2,3,5,7,251,65537];
  lwi:=ListWithIdenticalEntries;;

Old code:

gap> x:=2^20;; ListX(lwi(10000,x),ps,PValuation);; time;
90
gap> x:=2^1000-1;; ListX(lwi(10000,x),ps,PValuation);; time;
79
gap> x:=2^1000;; ListX(lwi(10000,x),ps,PValuation);; time;
2608
gap> x:=2^1000+1;; ListX(lwi(10000,x),ps,PValuation);; time;
56
gap> x:=2^1000+2^100;; ListX(lwi(10000,x),ps,PValuation);; time;
379

New code:

gap> x:=2^20;; ListX(lwi(10000,x),ps,PValuation);; time;
40
gap> x:=2^1000-1;; ListX(lwi(10000,x),ps,PValuation);; time;
86
gap> x:=2^1000;; ListX(lwi(10000,x),ps,PValuation);; time;
64
gap> x:=2^1000+1;; ListX(lwi(10000,x),ps,PValuation);; time;
71
gap> x:=2^1000+2^100;; ListX(lwi(10000,x),ps,PValuation);; time;
67

Using PVALUATION_INT directly (returns 0 for n=0, and only
supports integer inputs):

gap> x:=2^20;; ListX(lwi(10000,x),ps,PVALUATION_INT);; time;
26
gap> x:=2^1000-1;; ListX(lwi(10000,x),ps,PVALUATION_INT);; time;
67
gap> x:=2^1000;; ListX(lwi(10000,x),ps,PVALUATION_INT);; time;
47
gap> x:=2^1000+1;; ListX(lwi(10000,x),ps,PVALUATION_INT);; time;
53
gap> x:=2^1000+2^100;; ListX(lwi(10000,x),ps,PVALUATION_INT);; time;
49

lib/numtheor.gi

@@ -1392,21 +1392,15 @@ end );
 
 InstallGlobalFunction(PValuation,function(n,p)
 local v;
-  if not IsPrimeInt(p) or not IsRat(n) then
+  if not IsPosInt(p) or not IsRat(n) then
     Error("wrong parameters");
   fi;
-  if IsZero(n) then 
+  if n = 0 then
     return infinity;
-  elif not IsInt(n) then
-    return PValuation(NumeratorRat(n),p)-PValuation(DenominatorRat(n),p);
-  elif n<0 then n:=-n;
+  elif IsInt(n) then
+    return PVALUATION_INT(n,p);
   fi;
-  v:=0;
-  while n mod p=0 do
-    v:=v+1;
-    n:=n/p;
-  od;
-  return v;
+  return PVALUATION_INT(NumeratorRat(n),p) - PVALUATION_INT(DenominatorRat(n),p);
 end);
 
 #T ##########################################################################

src/gmpints.c

@@ -1988,6 +1988,7 @@ Obj JacobiInt ( Obj opL, Obj opR )
   return INTOBJ_INT( result );
 }
 
+
 /****************************************************************************
 **
 */
@@ -1998,6 +1999,52 @@ Obj FuncJACOBI_INT ( Obj self, Obj opL, Obj opR )
   return JacobiInt( opL, opR );
 }
 
+
+/****************************************************************************
+**
+*/
+Obj PValuationInt ( Obj n, Obj p )
+{
+  fake_mpz_t mpzN, mpzP;
+  mpz_t mpzResult;
+  int k;
+
+  CHECK_INT(n);
+  CHECK_INT(p);
+
+  if ( p == INTOBJ_INT(0) ) {
+    ErrorMayQuit( "PValuationInt: <p> must be nonzero", 0L, 0L  );
+  }
+
+  /* For certain values of p, mpz_remove replaces its "dest" argument
+     and tries to deallocate the original mpz_t in it. This means
+     we cannot use a fake_mpz_t for it. However, we are not really
+     interested in it anyway. */
+  mpz_init( mpzResult );
+  FAKEMPZ_GMPorINTOBJ( mpzN, n );
+  FAKEMPZ_GMPorINTOBJ( mpzP, p );
+
+  k = mpz_remove( mpzResult, MPZ_FAKEMPZ(mpzN), MPZ_FAKEMPZ(mpzP) );
+  CHECK_FAKEMPZ(mpzN);
+  CHECK_FAKEMPZ(mpzP);
+
+  /* throw away mpzResult -- it equals m / p^k */
+  mpz_clear( mpzResult );
+
+  return INTOBJ_INT( k );
+}
+
+/****************************************************************************
+**
+*/
+Obj FuncPVALUATION_INT ( Obj self, Obj opL, Obj opR )
+{
+  REQUIRE_INT_ARG( "PValuationInt", "left", opL );
+  REQUIRE_INT_ARG( "PValuationInt", "right", opR );
+  return PValuationInt( opL, opR );
+}
+
+
 /****************************************************************************
 **
 */
@@ -2247,6 +2294,9 @@ static StructGVarFunc GVarFuncs [] = {
   { "JACOBI_INT", 2, "gmp1, gmp2",
     FuncJACOBI_INT, "src/gmpints.c:JACOBI_INT" },
 
+  { "PVALUATION_INT", 2, "n, p",
+    FuncPVALUATION_INT, "src/gmpints.c:PVALUATION_INT" },
+
   { "POWERMODINT", 3, "base, exp, mod",
     FuncPOWERMODINT, "src/gmpints.c:POWERMODINT" },
 

tst/testinstall/intarith.tst

@@ -1177,6 +1177,19 @@ gap> for m in [1..100] do
 >   od;
 > od;
 
+#
+# test PVALUATION_INT
+#
+gap> checkPValuationInt:=function(n,p)
+>   local k, m;
+>   if n = 0 then return true; fi;
+>   k:=PVALUATION_INT(n,p);
+>   m:=n/p^k;
+>   return IsInt(m) and (m mod p) <> 0;
+> end;;
+gap> ForAll([-10000 .. 10000], n-> ForAll([2,3,5,7,251], p -> checkPValuationInt(n,p)));
+true
+
 #
 gap> STOP_TEST( "intarith.tst", 290000);
 

tst/testinstall/numtheor.tst

@@ -3,6 +3,10 @@ gap> START_TEST("numtheor.tst");
 #
 gap> PValuation(0,2);
 infinity
+gap> PValuation(0,3);
+infinity
+
+#
 gap> PValuation(100,2);
 2
 gap> PValuation(100,3);