Skip to content

Compute modular inverses of integers faster #2053

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.

Sign up for GitHub

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

Merged
fingolfin merged 4 commits into from
Jan 9, 2018
Merged

Compute modular inverses of integers faster #2053

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e415fb2
kernel: use AInvInt in some more places
fingolfin Dec 27, 2017
e683d90
kernel: clang-format InverseModInt
fingolfin Dec 28, 2017
14d9ca1
kernel: speed up InverseModInt for small inputs
fingolfin Dec 28, 2017
79ec5f3
kernel: simplify and sped up `ModRat` using `InverseModInt`
fingolfin Dec 27, 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
9 changes: 4 additions & 5 deletions src/dteval.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -332,8 +332,7 @@ Obj Power(
{
y = NEW_PLIST( T_PLIST, 0);
SET_LEN_PLIST(y, 0);
return Power( Solution(x, y, dtpols),
ProdInt(INTOBJ_INT(-1), n), dtpols );
return Power( Solution(x, y, dtpols), AInvInt(n), dtpols );
}
res = NEW_PLIST(T_PLIST, 2);
SET_LEN_PLIST(res, 0);
Expand Down Expand Up @@ -398,7 +397,7 @@ Obj Solution( Obj x,
}
else if ( CELM(x, j) < CELM(y, k) )
{
m = ProdInt( INTOBJ_INT(-1), ELM_PLIST(x, j+1) );
m = AInvInt( ELM_PLIST(x, j+1) );
SET_ELM_PLIST( res, i, ELM_PLIST(x, j) );
SET_ELM_PLIST( res, i+1, m );
CHANGED_BAG( res );
Expand All @@ -415,7 +414,7 @@ Obj Solution( Obj x,
if ( j < len1 )
while( j < len1 )
{
m = ProdInt( INTOBJ_INT(-1), ELM_PLIST( x, j+1 ) );
m = AInvInt( ELM_PLIST( x, j+1 ) );
SET_ELM_PLIST( res, i, ELM_PLIST(x, j) );
SET_ELM_PLIST( res, i+1, m );
CHANGED_BAG( res );
Expand Down Expand Up @@ -469,7 +468,7 @@ Obj Solution( Obj x,
}
else if ( !IS_INTOBJ( ELM_PLIST(xk, i) ) || CELM( xk, i ) != 0 )
{
m = ProdInt( INTOBJ_INT(-1), ELM_PLIST(xk, i) );
m = AInvInt( ELM_PLIST(xk, i) );
SET_ELM_PLIST( res, j, INTOBJ_INT(i) );
SET_ELM_PLIST( res, j+1, m );
CHANGED_BAG(res);
Expand Down
72 changes: 49 additions & 23 deletions src/integer.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2423,37 +2423,63 @@ Obj FuncPVALUATION_INT(Obj self, Obj n, Obj p)
/****************************************************************************
**
*/
Obj InverseModInt ( Obj base, Obj mod )
Obj InverseModInt(Obj base, Obj mod)
{
fake_mpz_t base_mpz, mod_mpz, result_mpz;
int success;
fake_mpz_t base_mpz, mod_mpz, result_mpz;
int success;

CHECK_INT(base);
CHECK_INT(mod);
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;
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 );
// handle small inputs separately
if (IS_INTOBJ(mod)) {

success = mpz_invert( MPZ_FAKEMPZ(result_mpz),
MPZ_FAKEMPZ(base_mpz),
MPZ_FAKEMPZ(mod_mpz) );
Int a = INT_INTOBJ(mod);
if (a < 0)
a = -a;

if (!success)
return Fail;
Int b = INT_INTOBJ(ModInt(base, mod));

CHECK_FAKEMPZ(result_mpz);
CHECK_FAKEMPZ(base_mpz);
CHECK_FAKEMPZ(mod_mpz);
Int aL = 0; // cofactor of a
Int bL = 1; // cofactor of b

return GMPorINTOBJ_FAKEMPZ( result_mpz );
// extended Euclidean algorithm
while (b != 0) {
Int hdQ = a / b;
Int c = b;
Int cL = bL;
b = a - hdQ * b;
bL = aL - hdQ * bL;
a = c;
aL = cL;
}
if (a != 1)
return Fail;
return ModInt(INTOBJ_INT(aL), mod);
}

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);
}

/****************************************************************************
Expand Down
8 changes: 8 additions & 0 deletions src/integer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -335,6 +335,14 @@ extern Obj GcdInt( Obj opL, Obj opR );
extern Obj AInvInt( Obj op );


/****************************************************************************
**
*F InverseModInt( <op> ) . . . . mult. inverse of an integer modulo another
**
*/
extern Obj InverseModInt(Obj base, Obj mod);


/****************************************************************************
**
** Compute log2 of the absolute value of an Int, i.e. the index of the highest
Expand Down
4 changes: 2 additions & 2 deletions src/objcftl.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,7 +148,7 @@ Obj CollectPolycyc (
if( (y = GET_IPOWER( h )) ) {
e = QuoInt( s, e );
if( !IS_INT_ZERO( t ) ) e = DiffInt( e, INTOBJ_INT(1) );
e = ProdInt( e, INTOBJ_INT(-1) );
e = AInvInt(e);
}
}
}
Expand Down Expand Up @@ -209,7 +209,7 @@ Obj CollectPolycyc (
ge = QuoInt( ge, e );
if( !IS_INT_ZERO( mge ) )
ge = DiffInt( ge, INTOBJ_INT(1) );
ge = ProdInt( ge, INTOBJ_INT(-1) );
ge = AInvInt(ge);
}
}
}
Expand Down
4 changes: 2 additions & 2 deletions src/permutat.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1177,7 +1177,7 @@ Obj PowPerm2Int (
ptKnown[p] = 0;

/* get pointer to the permutation and the power */
opR = ProdInt( INTOBJ_INT(-1), opR );
opR = AInvInt(opR);
ptL = CONST_ADDR_PERM2(opL);
ptP = ADDR_PERM2(pow);

Expand Down Expand Up @@ -1439,7 +1439,7 @@ Obj PowPerm4Int (
ptKnown[p] = 0;

/* get pointer to the permutation and the power */
opR = ProdInt( INTOBJ_INT(-1), opR );
opR = AInvInt(opR);
ptL = CONST_ADDR_PERM4(opL);
ptP = ADDR_PERM4(pow);

Expand Down
69 changes: 23 additions & 46 deletions src/rational.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,7 @@ Obj AInvRat (
Obj tmp;
CHECK_RAT(op);
res = NewBag( T_RAT, 2 * sizeof(Obj) );
tmp = AINV( NUM_RAT(op) );
tmp = AInvInt( NUM_RAT(op) );
SET_NUM_RAT(res, tmp);
SET_DEN_RAT(res, DEN_RAT(op));
CHANGED_BAG(res);
Expand Down Expand Up @@ -572,8 +572,8 @@ Obj QuoRat (
/* we multiply the left numerator with the right denominator */
/* so the right denominator should carry the sign of the right operand */
if ( IS_NEG_INT(numR) ) {
numR = ProdInt( INTOBJ_INT( -1L ), numR );
denR = ProdInt( INTOBJ_INT( -1L ), denR );
numR = AInvInt( numR );
denR = AInvInt( denR );
}

/* find the gcds */
Expand Down Expand Up @@ -611,10 +611,10 @@ Obj QuoRat (

/****************************************************************************
**
*F ModRat( <opL>, <opR> ) . . . . . . . . remainder of fraction mod integer
*F ModRat( <opL>, <n> ) . . . . . . . . remainder of fraction mod integer
**
** 'ModRat' returns the remainder of the fraction <opL> modulo the integer
** <opR>. The remainder is always an integer.
** <n>. The remainder is always an integer.
**
** '<r> / <s> mod <n>' yields the remainder of the fraction '<p> / <q>'
** modulo the integer '<n>', where '<p> / <q>' is the reduced form of
Expand All @@ -634,43 +634,20 @@ Obj QuoRat (
** such that $0 \<= t/s \< n$ and $r/s - t/s$ is a multiple of $n$. This is
** rarely needed while computing modular inverses is very useful.
*/
Obj ModRat (
Obj opL,
Obj opR )
Obj ModRat(Obj opL, Obj n)
{
Obj a, aL, b, bL, c, cL, hdQ;

/* make the integer positive */
if ( IS_NEG_INT(opR) ) {
opR = ProdInt( INTOBJ_INT( -1L ), opR );
}

/* let <p>/<q> represent <r>/<s> in reduced form */
/* invert the denominator <q> modulo <n> with Euclids algorithm */
a = opR; aL = INTOBJ_INT( 0L ); /* a = <n> */
b = DEN_RAT(opL); bL = INTOBJ_INT( 1L ); /* b = <q> */
while ( a != INTOBJ_INT( 1L ) ) {
while ( b != INTOBJ_INT( 0L ) ) {
hdQ = QuoInt( a, b );
c = b; cL = bL;
b = DiffInt( a, ProdInt( hdQ, b ) );
bL = DiffInt( aL, ProdInt( hdQ, bL ) );
a = c; aL = cL;
}

/* check whether the denominator <q> really was invertible mod <n> */
if ( a != INTOBJ_INT( 1L ) ) {
opR = ErrorReturnObj(
"ModRat: for <r>/<s> mod <n>, <s>/gcd(<r>,<s>) and <n> must be coprime",
0L, 0L,
"you can replace the integer <n> via 'return <n>;'" );
a = opR; aL = INTOBJ_INT( 0L ); /* a = <n> */
b = DEN_RAT(opL); bL = INTOBJ_INT( 1L ); /* b = <q> */
}
// invert the denominator
Obj d = InverseModInt( DEN_RAT(opL), n );

// check whether the denominator of <opL> really was invertible mod <n> */
if ( d == Fail ) {
ErrorMayQuit(
"ModRat: for <r>/<s> mod <n>, <s>/gcd(<r>,<s>) and <n> must be coprime",
0, 0 );
}

/* return the remainder */
return ModInt( ProdInt( NUM_RAT(opL), aL ), opR );
// return the remainder
return ModInt( ProdInt( NUM_RAT(opL), d ), n );
}


Expand Down Expand Up @@ -712,28 +689,28 @@ Obj PowRat (

/* if <opR> is negative and numerator is 1 just power the denominator */
else if ( NUM_RAT(opL) == INTOBJ_INT( 1L ) ) {
pow = PowInt( DEN_RAT(opL), ProdInt( INTOBJ_INT(-1L), opR ) );
pow = PowInt( DEN_RAT(opL), AInvInt( opR ) );
}

/* if <opR> is negative and numerator is -1 return (-1)^r * num(l) */
else if ( NUM_RAT(opL) == INTOBJ_INT( -1L ) ) {
numP = PowInt( NUM_RAT(opL), ProdInt( INTOBJ_INT( -1L ), opR ) );
denP = PowInt( DEN_RAT(opL), ProdInt( INTOBJ_INT( -1L ), opR ) );
numP = PowInt( NUM_RAT(opL), AInvInt( opR ) );
denP = PowInt( DEN_RAT(opL), AInvInt( opR ) );
pow = ProdInt(numP, denP);
}

/* if <opR> is negative do both powers, take care of the sign */
else {
numP = PowInt( DEN_RAT(opL), ProdInt( INTOBJ_INT( -1L ), opR ) );
denP = PowInt( NUM_RAT(opL), ProdInt( INTOBJ_INT( -1L ), opR ) );
numP = PowInt( DEN_RAT(opL), AInvInt( opR ) );
denP = PowInt( NUM_RAT(opL), AInvInt( opR ) );
pow = NewBag( T_RAT, 2 * sizeof(Obj) );
if ( IS_POS_INT(denP) ) {
SET_NUM_RAT(pow, numP);
SET_DEN_RAT(pow, denP);
}
else {
SET_NUM_RAT(pow, ProdInt( INTOBJ_INT( -1L ), numP ));
SET_DEN_RAT(pow, ProdInt( INTOBJ_INT( -1L ), denP ));
SET_NUM_RAT(pow, AInvInt( numP ));
SET_DEN_RAT(pow, AInvInt( denP ));
}
/* 'CHANGED_BAG' not needed, 'pow' is the youngest bag */
}
Expand Down
12 changes: 10 additions & 2 deletions tst/testinstall/intarith.tst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1242,7 +1242,13 @@ 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));
> if g = 1 then
> Assert(0, i in [0..m-1]);
> Assert(0, (i*b-1) mod m = 0); # formulated this way to support m=1
> Assert(0, i = INVMODINT(b,-m));
> else
> Assert(0, i = fail);
> fi;
> od;
> od;

Expand All @@ -1251,7 +1257,9 @@ gap> INVMODINT(1,0);
Error, InverseModInt: <mod> must be nonzero
gap> INVMODINT(2,10);
fail
gap> INVMODINT(2,10);
gap> INVMODINT(2,-10);
fail
gap> INVMODINT(2,2^80);
fail
gap> INVMODINT(2,1);
0
Expand Down
18 changes: 18 additions & 0 deletions tst/testinstall/rationals.tst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,10 @@ gap> NumeratorRat(1/2);
1
gap> DenominatorRat(1/2);
2

#
# 'Rat' -- converting things to a rational
#
gap> Rat(1/2) = 1/2;
true
gap> Rat(1) = 1;
Expand All @@ -19,3 +23,17 @@ gap> Rat(3.7e-1);
37/100

#
# modular inverses of rationals
#
gap> 1/2 mod 3;
2
gap> 1/2 mod -3;
2
gap> -1/2 mod 3;
1
gap> -1/2 mod -3;
1
gap> 1/2 mod 2;
Error, ModRat: for <r>/<s> mod <n>, <s>/gcd(<r>,<s>) and <n> must be coprime

#