From 85d63046125880787489b73c730819729868ffea Mon Sep 17 00:00:00 2001 From: Max Horn Date: Sat, 21 Jan 2017 15:47:04 +0100 Subject: [PATCH] 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 use PVALUATION_INT, and the Valuation method for integers. 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/integer.gi | 17 +++---------- lib/numtheor.gi | 16 ++++-------- src/gmpints.c | 49 ++++++++++++++++++++++++++++++++++++ tst/testinstall/intarith.tst | 13 ++++++++++ tst/testinstall/numtheor.tst | 4 +++ 5 files changed, 74 insertions(+), 25 deletions(-) diff --git a/lib/integer.gi b/lib/integer.gi index f6e2583f84..7f5e829961 100644 --- a/lib/integer.gi +++ b/lib/integer.gi @@ -1809,24 +1809,13 @@ InstallMethod( StandardAssociateUnit, InstallOtherMethod( Valuation, "for two integers", IsIdenticalObj, - [ IsInt, - IsInt ], + [ IsInt, IsInt ], 0, - function( n, m ) - local val; - if n = 0 then - val := infinity; - else - val := 0; - while n mod m = 0 do - val := val + 1; - n := n / m; - od; + return infinity; fi; - return val; - + return PVALUATION_INT( n, m ); end ); diff --git a/lib/numtheor.gi b/lib/numtheor.gi index a8ac64ccbf..4e53db23d2 100644 --- a/lib/numtheor.gi +++ b/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 IsInt(p) or not IsRat(n) or p = 0 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 ########################################################################## diff --git a/src/gmpints.c b/src/gmpints.c index d8f7647ad4..b3d911a86e 100644 --- a/src/gmpints.c +++ b/src/gmpints.c @@ -1983,6 +1983,7 @@ Obj JacobiInt ( Obj opL, Obj opR ) return INTOBJ_INT( result ); } + /**************************************************************************** ** */ @@ -1993,6 +1994,51 @@ 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:

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 ); +} + + /**************************************************************************** ** */ @@ -2237,6 +2283,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" }, diff --git a/tst/testinstall/intarith.tst b/tst/testinstall/intarith.tst index fc90d2c9a1..3e23552a8d 100644 --- a/tst/testinstall/intarith.tst +++ b/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); diff --git a/tst/testinstall/numtheor.tst b/tst/testinstall/numtheor.tst index 48e591585c..39407b3c76 100644 --- a/tst/testinstall/numtheor.tst +++ b/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);