Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

More GMP work #1075

Merged
merged 21 commits into from
Jan 23, 2017
Merged
Changes from 1 commit
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
1898ea7
gmpints: introduce fake_mpz_t
fingolfin Jan 21, 2017
8c018e0
gmpints: merge GMPorINTOBJ_INT and ObjInt_Int
fingolfin Jan 21, 2017
e918054
gmpints: rewrite PrintInt using mpz_get_str
fingolfin Jan 21, 2017
1c6d8ee
gmpints: rewrite HexStringInt and STRING_INT using fake_mpz_t
fingolfin Jan 21, 2017
b591686
gmpints: remove GMP_INTOBJ
fingolfin Jan 21, 2017
9a3689f
gmpints: rewrite FuncIntHexString
fingolfin Jan 21, 2017
99418d3
gmpints: replace TypLimb -> mp_limb_t, TypGMPSize -> mp_size_t
fingolfin Jan 21, 2017
b700ea4
gmpints: replace INTEGER_ALLOCATION_SIZE by INTEGER_UNIT_SIZE
fingolfin Jan 21, 2017
af567b1
gmpints: cleanup gmpints.h
fingolfin Jan 21, 2017
6a5815f
gmpints: simplify ModInt
fingolfin Jan 21, 2017
b8a188c
gmpints: remove NEW_INTPOS
fingolfin Jan 21, 2017
8454a6f
gmpints: add IS_INT, IS_POSITIVE, REQUIRE_INT_ARG
fingolfin Jan 21, 2017
1c8ae24
gmpints: change REQUIRE_INT_ARG to use ErrorMayQuit
fingolfin Jan 21, 2017
cdba5d4
gmpints: use REQUIRE_INT_ARG in two more places
fingolfin Jan 21, 2017
3e44492
gmpints: fix RemInt error handling
fingolfin Jan 21, 2017
42c4845
gmpints: add JACOBI_INT and change Jacobi() to use it
fingolfin Jan 21, 2017
0950124
gmpints: add POWERMODINT and INVMODINT
fingolfin Jan 21, 2017
9561e9d
gmpints: add PVALUATION_INT
fingolfin Jan 21, 2017
fbf656a
tests: replace legal header in intarith.tst
fingolfin Jan 21, 2017
9bc7f34
gmpints: explain the buffer in PrintInt
fingolfin Jan 23, 2017
cad52ab
gmpints: improve initial documentation comment
fingolfin Jan 23, 2017
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
gmpints: add POWERMODINT and INVMODINT
This adds a GMP based implementation of PowerModInt. The new kernel
function INVMODINT computes modular inverse, and is not currently
exposed to the library in a documented function.

As usual, we get a speedup. Old code:

gap> for b in [-100..100] do for e in [1..100] do for m in [1..100] do
>  PowerModInt(b,e,m);od;od;od; time;
2270
gap> for b in [-100..100] do for e in [1..100] do for m in [1..100] do
>  PowerModInt(b,2^e,m);od;od;od; time;
11354

New code:

gap> for b in [-100..100] do for e in [1..100] do for m in [1..100] do
>  PowerModInt(b,e,m);od;od;od; time;
625
gap> for b in [-100..100] do for e in [1..100] do for m in [1..100] do
>  PowerModInt(b,2^e,m);od;od;od; time;
1920
fingolfin committed Jan 22, 2017
commit 09501242ed7189cdab9fdb2ec8cb4832e2465a01
2 changes: 1 addition & 1 deletion lib/integer.gd
Original file line number Diff line number Diff line change
@@ -799,7 +799,7 @@ DeclareGlobalFunction( "NextPrimeInt" );
##
## <Description>
## returns <M><A>r</A>^{<A>e</A>} \pmod{<A>m</A>}</M> for integers <A>r</A>,
## <A>e</A> and <A>m</A> (<M><A>e</A> \geq 0</M>).
## <A>e</A> and <A>m</A>.
## <P/>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it obvious enough that e 0 ?

Negative values seem to work (now) if r is invertible, so that should maybe be documented explicitly?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But e can (still) be 0. And negative values were supported before, just not documented.
Documenting this more explicitly is fine by me, but I'd prefer to not do this as part of this PR (which has been in the works for very long by now)...

## Note that <Ref Func="PowerModInt"/> can reduce intermediate results and
## thus will generally be faster than using
35 changes: 1 addition & 34 deletions lib/integer.gi
Original file line number Diff line number Diff line change
@@ -1161,40 +1161,7 @@ end);
##
#F PowerModInt(<r>,<e>,<m>) . . . . . . power of one integer modulo another
##
InstallGlobalFunction(PowerModInt,function ( r, e, m )
local pow, f;

# handle special cases
if m = 1 then
return 0;
elif e = 0 then
return 1;
fi;

# reduce `r' initially
r := r mod m;

# if `e' is negative then invert `r' modulo `m' with Euclids algorithm
if e < 0 then
r := 1/r mod m;
e := -e;
fi;

# now use the repeated squaring method (right-to-left)
pow := 1;
f := 2 ^ (LogInt( e, 2 ) + 1);
while 1 < f do
pow := (pow * pow) mod m;
f := QuoInt( f, 2 );
if f <= e then
pow := (pow * r) mod m;
e := e - f;
fi;
od;

# return the power
return pow;
end);
InstallGlobalFunction(PowerModInt, POWERMODINT);


#############################################################################
106 changes: 105 additions & 1 deletion src/gmpints.c
Original file line number Diff line number Diff line change
@@ -1993,6 +1993,104 @@ Obj FuncJACOBI_INT ( Obj self, Obj opL, Obj opR )
return JacobiInt( opL, opR );
}

/****************************************************************************
**
*/
Obj InverseModInt ( Obj base, Obj mod )
{
fake_mpz_t base_mpz, mod_mpz, result_mpz;
int success;

CHECK_INT(base);
CHECK_INT(mod);

if ( mod == INTOBJ_INT(0) )
ErrorMayQuit( "InverseModInt: <mod> must be nonzero", 0L, 0L );
if ( mod == INTOBJ_INT(1) || mod == INTOBJ_INT(-1) )
return INTOBJ_INT(0);
if ( base == INTOBJ_INT(0) )
return Fail;

NEW_FAKEMPZ( result_mpz, SIZE_INT_OR_INTOBJ(mod) + 1 );
FAKEMPZ_GMPorINTOBJ( base_mpz, base );
FAKEMPZ_GMPorINTOBJ( mod_mpz, mod );

success = mpz_invert( MPZ_FAKEMPZ(result_mpz),
MPZ_FAKEMPZ(base_mpz),
MPZ_FAKEMPZ(mod_mpz) );

if (!success)
return Fail;

CHECK_FAKEMPZ(result_mpz);
CHECK_FAKEMPZ(base_mpz);
CHECK_FAKEMPZ(mod_mpz);

return GMPorINTOBJ_FAKEMPZ( result_mpz );
}

/****************************************************************************
**
*/
Obj FuncINVMODINT ( Obj self, Obj base, Obj mod )
{
REQUIRE_INT_ARG( "InverseModInt", "base", base );
REQUIRE_INT_ARG( "InverseModInt", "mod", mod );
return InverseModInt( base, mod );
}


/****************************************************************************
**
*/
Obj PowerModInt ( Obj base, Obj exp, Obj mod )
{
fake_mpz_t base_mpz, exp_mpz, mod_mpz, result_mpz;

CHECK_INT(base);
CHECK_INT(exp);
CHECK_INT(mod);

if ( mod == INTOBJ_INT(0) )
ErrorMayQuit( "PowerModInt: <mod> must be nonzero", 0L, 0L );
if ( mod == INTOBJ_INT(1) || mod == INTOBJ_INT(-1) )
return INTOBJ_INT(0);

if ( IS_NEGATIVE(exp) ) {
base = InverseModInt( base, mod );
if (base == Fail)
ErrorMayQuit( "PowerModInt: negative <exp> but <base> is not invertible modulo <mod>", 0L, 0L );
exp = AInvInt(exp);
}

NEW_FAKEMPZ( result_mpz, SIZE_INT_OR_INTOBJ(mod) );
FAKEMPZ_GMPorINTOBJ( base_mpz, base );
FAKEMPZ_GMPorINTOBJ( exp_mpz, exp );
FAKEMPZ_GMPorINTOBJ( mod_mpz, mod );

mpz_powm( MPZ_FAKEMPZ(result_mpz), MPZ_FAKEMPZ(base_mpz),
MPZ_FAKEMPZ(exp_mpz), MPZ_FAKEMPZ(mod_mpz) );

CHECK_FAKEMPZ(result_mpz);
CHECK_FAKEMPZ(base_mpz);
CHECK_FAKEMPZ(exp_mpz);
CHECK_FAKEMPZ(mod_mpz);

return GMPorINTOBJ_FAKEMPZ( result_mpz );
}

/****************************************************************************
**
*/
Obj FuncPOWERMODINT ( Obj self, Obj base, Obj exp, Obj mod )
{
REQUIRE_INT_ARG( "PowerModInt", "base", base );
REQUIRE_INT_ARG( "PowerModInt", "exp", exp );
REQUIRE_INT_ARG( "PowerModInt", "mod", mod );
return PowerModInt( base, exp, mod );
}


/****************************************************************************
**
** * * * * * * * "Mersenne twister" random numbers * * * * * * * * * * * * *
@@ -2138,7 +2236,13 @@ static StructGVarFunc GVarFuncs [] = {

{ "JACOBI_INT", 2, "gmp1, gmp2",
FuncJACOBI_INT, "src/gmpints.c:JACOBI_INT" },


{ "POWERMODINT", 3, "base, exp, mod",
FuncPOWERMODINT, "src/gmpints.c:POWERMODINT" },

{ "INVMODINT", 2, "base, mod",
FuncINVMODINT, "src/gmpints.c:INVMODINT" },

{ "HexStringInt", 1, "gmp",
FuncHexStringInt, "src/gmpints.c:HexStringInt" },

89 changes: 89 additions & 0 deletions tst/testinstall/intarith.tst
Original file line number Diff line number Diff line change
@@ -1154,6 +1154,95 @@ gap> Print(2^(64*1050), "\n");
948688599817481827330434634041214374947601458440559056883819756862448124537321\
8157000453169824579899621376

#
# test InversModInt / INVMODINT
#
gap> for m in [1..100] do
> for b in [1..100] do
> i := INVMODINT(b,m);
> g := GcdInt(b,m);
> Assert(0, (g<>1 and i = fail) or (g=1 and i in [0..m-1] and (i*b-1) mod m = 0));
> od;
> od;

#
gap> INVMODINT(1,0);
Error, InverseModInt: <mod> must be nonzero
gap> INVMODINT(2,10);
fail
gap> INVMODINT(2,10);
fail
gap> INVMODINT(2,1);
0
gap> INVMODINT(2,-1);
0

#
# test PowerModInt
#
gap> for m in [1..100] do
> for b in [-30..30] do
> for e in [0..10] do
> Assert(0, PowerModInt(b,e,m)=(b^e mod m));
> od;
> od;
> od;

#
gap> PowerModInt(1,1,0);
Error, PowerModInt: <mod> must be nonzero
gap> PowerModInt(0,-1,5);
Error, PowerModInt: negative <exp> but <base> is not invertible modulo <mod>
gap> PowerModInt(2,-1,5);
3
gap> PowerModInt(2,-1,10);
Error, PowerModInt: negative <exp> but <base> is not invertible modulo <mod>

# compare C and GAP implementation
gap> POWERMODINT_GAP := function ( r, e, m )
> local pow, f;
>
> # handle special cases
> if m = 1 then
> return 0;
> elif e = 0 then
> return 1;
> fi;
>
> # reduce `r' initially
> r := r mod m;
>
> # if `e' is negative then invert `r' modulo `m' with Euclids algorithm
> if e < 0 then
> r := 1/r mod m;
> e := -e;
> fi;
>
> # now use the repeated squaring method (right-to-left)
> pow := 1;
> f := 2 ^ (LogInt( e, 2 ) + 1);
> while 1 < f do
> pow := (pow * pow) mod m;
> f := QuoInt( f, 2 );
> if f <= e then
> pow := (pow * r) mod m;
> e := e - f;
> fi;
> od;
>
> # return the power
> return pow;
> end;;
gap> for m in [1..100] do
> for b in [-30..30] do
> for e in [0..30] do
> Assert(0, POWERMODINT_GAP(b,e,m)=PowerModInt(b,e,m));
> if GcdInt(m,b) = 1 then
> Assert(0, POWERMODINT_GAP(b,-e,m)=PowerModInt(b,-e,m));
> fi;
> od;
> od;
> od;

#
gap> STOP_TEST( "intarith.tst", 290000);