diff --git a/src/dteval.c b/src/dteval.c index 398f6c5422..55e64e95d8 100644 --- a/src/dteval.c +++ b/src/dteval.c @@ -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); @@ -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 ); @@ -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 ); @@ -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); diff --git a/src/integer.c b/src/integer.c index 7989bc2109..41dc7d2a73 100644 --- a/src/integer.c +++ b/src/integer.c @@ -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: 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: 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); } /**************************************************************************** diff --git a/src/integer.h b/src/integer.h index 754bea812b..4a8b155525 100644 --- a/src/integer.h +++ b/src/integer.h @@ -335,6 +335,14 @@ extern Obj GcdInt( Obj opL, Obj opR ); extern Obj AInvInt( Obj op ); +/**************************************************************************** +** +*F InverseModInt( ) . . . . 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 diff --git a/src/objcftl.c b/src/objcftl.c index 42051bc930..31511f20f1 100644 --- a/src/objcftl.c +++ b/src/objcftl.c @@ -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); } } } @@ -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); } } } diff --git a/src/permutat.c b/src/permutat.c index d6e4dd1370..bb8c097449 100644 --- a/src/permutat.c +++ b/src/permutat.c @@ -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); @@ -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); diff --git a/src/rational.c b/src/rational.c index 33b5c3ed1d..7fbe0f4e46 100644 --- a/src/rational.c +++ b/src/rational.c @@ -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); @@ -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 */ @@ -611,10 +611,10 @@ Obj QuoRat ( /**************************************************************************** ** -*F ModRat( , ) . . . . . . . . remainder of fraction mod integer +*F ModRat( , ) . . . . . . . . remainder of fraction mod integer ** ** 'ModRat' returns the remainder of the fraction modulo the integer -** . The remainder is always an integer. +** . The remainder is always an integer. ** ** ' / mod ' yields the remainder of the fraction '

/ ' ** modulo the integer '', where '

/ ' is the reduced form of @@ -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

/ represent / in reduced form */ - /* invert the denominator modulo with Euclids algorithm */ - a = opR; aL = INTOBJ_INT( 0L ); /* a = */ - b = DEN_RAT(opL); bL = INTOBJ_INT( 1L ); /* b = */ - 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 really was invertible mod */ - if ( a != INTOBJ_INT( 1L ) ) { - opR = ErrorReturnObj( - "ModRat: for / mod , /gcd(,) and must be coprime", - 0L, 0L, - "you can replace the integer via 'return ;'" ); - a = opR; aL = INTOBJ_INT( 0L ); /* a = */ - b = DEN_RAT(opL); bL = INTOBJ_INT( 1L ); /* b = */ - } + // invert the denominator + Obj d = InverseModInt( DEN_RAT(opL), n ); + + // check whether the denominator of really was invertible mod */ + if ( d == Fail ) { + ErrorMayQuit( + "ModRat: for / mod , /gcd(,) and 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 ); } @@ -712,28 +689,28 @@ Obj PowRat ( /* if 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 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 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 */ } diff --git a/tst/testinstall/intarith.tst b/tst/testinstall/intarith.tst index 2341bfa219..8f18bd6cbb 100644 --- a/tst/testinstall/intarith.tst +++ b/tst/testinstall/intarith.tst @@ -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; @@ -1251,7 +1257,9 @@ gap> INVMODINT(1,0); Error, InverseModInt: 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 diff --git a/tst/testinstall/rationals.tst b/tst/testinstall/rationals.tst index 00ad69a80a..22e6ede093 100644 --- a/tst/testinstall/rationals.tst +++ b/tst/testinstall/rationals.tst @@ -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; @@ -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 / mod , /gcd(,) and must be coprime + +#