#define RequireNonzero(funcname, op, argname) \ do { \ if (op == INTOBJ_INT(0)) { \
RequireArgumentEx(funcname, op, "<" argname ">", \ "must be a nonzero integer"); \
} \
} while (0)
GAP_STATIC_ASSERT( sizeof(mp_limb_t) == sizeof(UInt), "gmp limb size incompatible with GAP word size");
/* This ensures that all memory underlying a bag is actually committed **tophysicalmemoryandcanbewrittento. **ThisisaworkaroundtoabugspecifictoCygwin64-bitandbad **interactionwithGMP,sothisisonlyneededspecificallyfornew **bagscreatedinthismoduletoholdtheoutputsofGMProutines. ** **Thus,anytimeNewBagiscalled,itisalsonecessarytocall **ENSURE_BAG(bag)onthenewlycreatedbagifsomeGMPfunctionwillbe **thefirstplacethatbag'sdataiswrittento. ** **Togiveacounter-example,ENSURE_BAGis*not*neededinObjInt_Int, **becauseitjustcreatesabagtoholdasinglemp_limb_t,and **immediatelyassignsitavalue. ** **Thebugthisworksaroundisexplainedmorein **https://github.com/gap-system/gap/issues/3434
*/ staticinlinevoid ENSURE_BAG(Bag bag)
{ // Note: This workaround is only required with the original GMP and not with // MPIR #ifdefined(SYS_IS_CYGWIN32) && defined(SYS_IS_64_BIT) && \
!defined(__MPIR_VERSION)
memset(PTR_BAG(bag), 0, SIZE_BAG(bag)); #endif
}
// for fallbacks to library static Obj String; static Obj IsIntFilt;
// UPDATE_FAKEMPZ is a helper function for the MPZ_FAKEMPZ macro staticinlinevoid UPDATE_FAKEMPZ( fake_mpz_t fake )
{
fake->v->_mp_d = fake->obj ? (mp_ptr)ADDR_INT(fake->obj) : &fake->tmp;
}
/**************************************************************************** ** **ThisisahelperfunctionfortheCHECK_INTmacro,whichchecksthat **thegivenintegerobject<op>isnormalizedandreduced. **
*/ #if DEBUG_GMP staticBOOL IS_NORMALIZED_AND_REDUCED(Obj op, constchar * func, int line)
{
mp_size_t size; if ( IS_INTOBJ( op ) ) { returnTRUE;
} if ( !IS_LARGEINT( op ) ) { // ignore non-integers returnFALSE;
} for ( size = SIZE_INT(op); size != (mp_size_t)1; size-- ) { if ( CONST_ADDR_INT(op)[(size - 1)] != 0 ) { break;
}
} if ( size < SIZE_INT(op) ) {
Pr("WARNING: non-normalized gmp value (%s:%d)\n",(Int)func,line);
} if ( SIZE_INT(op) == 1) { if ( ( VAL_LIMB0(op) <= INT_INTOBJ_MAX ) ||
( IS_INTNEG(op) && VAL_LIMB0(op) == -INT_INTOBJ_MIN ) ) { if ( IS_INTNEG(op) ) {
Pr("WARNING: non-reduced negative gmp value (%s:%d)\n",(Int)func,line); returnFALSE;
} else {
Pr("WARNING: non-reduced positive gmp value (%s:%d)\n",(Int)func,line); returnFALSE;
}
}
} returnTRUE;
} #endif
/**************************************************************************** ** *FObjInt_Int(<cint>)..........convertcinttointegerobject ** **'ObjInt_Int'takestheCinteger<cint>andreturnstheequivalent **GMPobjorintobj,accordingtothevalueof<cint>. **
*/
Obj ObjInt_Int( Int i )
{
Obj gmp;
if (INT_INTOBJ_MIN <= i && i <= INT_INTOBJ_MAX) { return INTOBJ_INT(i);
} elseif (i < 0 ) {
gmp = NewBag( T_INTNEG, sizeof(mp_limb_t) );
i = -i;
} else {
gmp = NewBag( T_INTPOS, sizeof(mp_limb_t) );
}
SET_VAL_LIMB0( gmp, i ); return gmp;
}
Obj ObjInt_UInt( UInt i )
{
Obj gmp; if (i <= INT_INTOBJ_MAX) { return INTOBJ_INT(i);
} else {
gmp = NewBag( T_INTPOS, sizeof(mp_limb_t) );
SET_VAL_LIMB0( gmp, i ); return gmp;
}
}
// initialize to -i
Obj ObjInt_UIntInv( UInt i )
{
Obj gmp; // we need to test INT_INTOBJ_MIN <= -i; to express this with unsigned // values, we must avoid all negative terms, which leads to this equivalent // check: if (i <= -INT_INTOBJ_MIN) { return INTOBJ_INT(-i);
} else {
gmp = NewBag( T_INTNEG, sizeof(mp_limb_t) );
SET_VAL_LIMB0( gmp, i ); return gmp;
}
}
Obj ObjInt_Int8( Int8 i )
{ #ifdef SYS_IS_64_BIT return ObjInt_Int(i); #else if (i == (Int4)i) { return ObjInt_Int((Int4)i);
}
// we need two limbs to store this integer
assert( sizeof(mp_limb_t) == 4 );
Obj gmp; if (i >= 0) {
gmp = NewBag( T_INTPOS, 2 * sizeof(mp_limb_t) );
} else {
gmp = NewBag( T_INTNEG, 2 * sizeof(mp_limb_t) );
i = -i;
}
Obj ObjInt_UInt8( UInt8 i )
{ #ifdef SYS_IS_64_BIT return ObjInt_UInt(i); #else if (i == (UInt4)i) { return ObjInt_UInt((UInt4)i);
}
// we need two limbs to store this integer
assert( sizeof(mp_limb_t) == 4 );
Obj gmp = NewBag( T_INTPOS, 2 * sizeof(mp_limb_t) );
UInt *ptr = ADDR_INT(gmp);
ptr[0] = (UInt4)i;
ptr[1] = ((UInt8)i) >> 32; return gmp; #endif
}
/**************************************************************************** ** **ConvertGAPIntegerstovariousCtypes--seeheaderfile
*/ Int Int_ObjInt(Obj i)
{
UInt sign = 0; if (IS_INTOBJ(i)) return INT_INTOBJ(i); // must be a single limb if (TNUM_OBJ(i) == T_INTPOS)
sign = 0; elseif (TNUM_OBJ(i) == T_INTNEG)
sign = 1; else
RequireArgument("Conversion error", i, "must be an integer"); if (SIZE_BAG(i) != sizeof(mp_limb_t))
ErrorMayQuit("Conversion error: integer too large", 0, 0);
// now check if val is small enough to fit in the signed Int type // that has a range from -2^N to 2^N-1 so we need to check both ends // Since -2^N is the same bit pattern as the UInt 2^N (N is 31 or 63) // we can do it as below which avoids some compiler warnings
UInt val = VAL_LIMB0(i); #ifdef SYS_IS_64_BIT if ((!sign && (val > INT64_MAX)) || (sign && (val > (UInt)INT64_MIN))) #else if ((!sign && (val > INT32_MAX)) || (sign && (val > (UInt)INT32_MIN))) #endif
ErrorMayQuit("Conversion error: integer too large", 0, 0); return sign ? -(Int)val : (Int)val;
}
UInt UInt_ObjInt(Obj i)
{ if (IS_NEG_INT(i))
ErrorMayQuit("Conversion error: cannot convert negative integer to unsigned type", 0, 0); if (IS_INTOBJ(i)) return (UInt)INT_INTOBJ(i); if (TNUM_OBJ(i) != T_INTPOS)
RequireArgument("Conversion error", i, "must be a non-negative integer");
// must be a single limb if (SIZE_INT(i) != 1)
ErrorMayQuit("Conversion error: integer too large", 0, 0); return VAL_LIMB0(i);
}
Int8 Int8_ObjInt(Obj i)
{ #ifdef SYS_IS_64_BIT // in this case Int8 is Int return Int_ObjInt(i); #else if (IS_INTOBJ(i)) return (Int8)INT_INTOBJ(i);
UInt sign = 0; if (TNUM_OBJ(i) == T_INTPOS)
sign = 0; elseif (TNUM_OBJ(i) == T_INTNEG)
sign = 1; else
RequireArgument("Conversion error", i, "must be an integer");
// must be at most two limbs if (SIZE_INT(i) > 2)
ErrorMayQuit("Conversion error: integer too large", 0, 0);
UInt vall = VAL_LIMB0(i);
UInt valh = (SIZE_INT(i) == 1) ? 0 : CONST_ADDR_INT(i)[1];
UInt8 val = (UInt8)vall + ((UInt8)valh << 32); // now check if val is small enough to fit in the signed Int8 type // that has a range from -2^63 to 2^63-1 so we need to check both ends // Since -2^63 is the same bit pattern as the UInt8 2^63 we can do it // this way which avoids some compiler warnings if ((!sign && (val > INT64_MAX)) || (sign && (val > (UInt8)INT64_MIN)))
ErrorMayQuit("Conversion error: integer too large", 0, 0); return sign ? -(Int8)val : (Int8)val; #endif
}
UInt8 UInt8_ObjInt(Obj i)
{ #ifdef SYS_IS_64_BIT // in this case UInt8 is UInt return UInt_ObjInt(i); #else if (IS_NEG_INT(i))
ErrorMayQuit("Conversion error: cannot convert negative integer to unsigned type", 0, 0); if (IS_INTOBJ(i)) return (UInt8)INT_INTOBJ(i); if (TNUM_OBJ(i) != T_INTPOS)
RequireArgument("Conversion error", i, "must be a non-negative integer"); if (SIZE_INT(i) > 2)
ErrorMayQuit("Conversion error: integer too large", 0, 0);
UInt vall = VAL_LIMB0(i);
UInt valh = (SIZE_INT(i) == 1) ? 0 : CONST_ADDR_INT(i)[1]; return (UInt8)vall + ((UInt8)valh << 32); #endif
}
// This function returns an immediate integer, or // an integer object with exactly one limb, and returns // its absolute value as an unsigned integer. staticinline UInt AbsOfSmallInt(Obj x)
{ if (!IS_INTOBJ(x)) {
GAP_ASSERT(SIZE_INT(x) == 1); return VAL_LIMB0(x);
} Int val = INT_INTOBJ(x); return val > 0 ? val : -val;
}
/**************************************************************************** ** *FPrintInt(<op>)................printanintegerobject ** **'PrintInt'printstheinteger<op>intheusualdecimalnotation.
*/ void PrintInt ( Obj op )
{ // print a small integer if ( IS_INTOBJ(op) ) {
Pr("%>%d%<", INT_INTOBJ(op), 0);
}
// print a large integer elseif ( SIZE_INT(op) < 1000 ) {
CHECK_INT(op);
/* Use GMP to print the integer to a buffer. We are looking at an integerwithlessthan1000limbs,henceindecimalnotation,it willtakeupatmostLogInt(2^(1000*GMP_LIMB_BITS),10) digits.Since1000*Log(2)/Log(10)=301.03,wegetthefollowing estimateforthebuffersize(theoverestimateisbigenoughto
include space for a sign and a null terminator). */ Char buf[302 * GMP_LIMB_BITS];
mpz_t v;
v->_mp_alloc = SIZE_INT(op);
v->_mp_size = IS_INTPOS(op) ? v->_mp_alloc : -v->_mp_alloc;
v->_mp_d = (mp_ptr)ADDR_INT(op);
mpz_get_str(buf, 10, v);
// print the buffer, %> means insert '\' before a linebreak
Pr("%>%s%<",(Int)buf, 0);
} else {
Obj str = CALL_1ARGS( String, op );
Pr("%>", 0, 0);
PrintString1(str);
Pr("%<", 0, 0); /* for a long time Print of large ints did not follow the general idea *thatPrintshouldproducesomethingthatcanbereadbackintoGAP:
Pr("<<an integer too large to be printed>>", 0, 0); */
}
}
/**************************************************************************** ** *FStringIntBase(<op>,<base>) ** **Converttheinteger<op>toastringrelativetothegivenbase<base>. **Here,basemayrangefrom2to36.
*/ static Obj StringIntBase(Obj op, int base)
{ int len;
Obj res;
fake_mpz_t v;
GAP_ASSERT(IS_INT(op));
CHECK_INT(op);
GAP_ASSERT(2 <= base && base <= 36);
// 0 is special if (op == INTOBJ_INT(0)) {
res = NEW_STRING(1);
CHARS_STRING(res)[0] = '0'; return res;
}
// convert integer to fake_mpz_t
FAKEMPZ_GMPorINTOBJ(v, op);
// allocate the result string
len = mpz_sizeinbase( MPZ_FAKEMPZ(v), base ) + 2;
res = NEW_STRING( len );
// ask GMP to perform the actual conversion
mpz_get_str( CSTR_STRING( res ), -base, MPZ_FAKEMPZ(v) );
// we may have to shrink the string int real_len = strlen( CONST_CSTR_STRING(res) ); if ( real_len != GET_LEN_STRING(res) ) {
SET_LEN_STRING(res, real_len);
}
Obj IntHexString(Obj str)
{
Obj res; Int i, len, sign, nd;
mp_limb_t n; const UInt1 *p;
UInt *limbs;
GAP_ASSERT(IS_STRING_REP(str));
len = GET_LEN_STRING(str); if (len == 0) {
res = INTOBJ_INT(0); return res;
}
p = CONST_CHARS_STRING(str); if (*p == '-') {
sign = -1;
i = 1;
} else {
sign = 1;
i = 0;
}
while (p[i] == '0' && i < len)
i++;
len -= i;
if (len*4 <= NR_SMALL_INT_BITS) {
n = hexstr2int( p + i, len );
res = INTOBJ_INT(sign * n); return res;
}
else { /* Each hex digit corresponds to 4 bits, and each GMP limb has sizeof(UInt) bytes,thus2*sizeof(UInt)hexdigitsfitintoonelimb.Weusethis
to compute the number of limbs minus 1: */
nd = (len - 1) / (2*sizeof(UInt));
res = NewBag( (sign == 1) ? T_INTPOS : T_INTNEG, (nd + 1) * sizeof(mp_limb_t) );
// update pointer, in case a garbage collection happened
p = CONST_CHARS_STRING(str) + i;
limbs = ADDR_INT(res);
// if len is not divisible by 2*sizeof(UInt), then take care of the extra bytes
UInt diff = len - nd * (2*sizeof(UInt)); if ( diff ) {
n = hexstr2int( p, diff );
p += diff;
len -= diff;
limbs[nd--] = n;
}
while ( len ) {
n = hexstr2int( p, 2*sizeof(UInt) );
p += 2*sizeof(UInt);
len -= 2*sizeof(UInt);
limbs[nd--] = n;
}
Int res = 0;
UInt b;
b = a >> 32; if (b) { res+=32; a=b; }
b = a >> 16; if (b) { res+=16; a=b; }
b = a >> 8; if (b) { res+= 8; a=b; } return res + LogTable256[a]; #endif
}
Int CLog2Int(Int a)
{ if (a == 0) return -1; if (a < 0) a = -a; return CLog2UInt(a);
}
if (IS_INTOBJ(n)) { return INTOBJ_INT(CLog2Int(INT_INTOBJ(n)));
}
UInt len = SIZE_INT(n) - 1;
UInt a = CLog2UInt(CONST_ADDR_INT(n)[len]);
CHECK_INT(n);
#ifdef SYS_IS_64_BIT return INTOBJ_INT(len * GMP_LIMB_BITS + a); #else /* The final result is len * GMP_LIMB_BITS - d, which may not
fit into an immediate integer (at least on a 32bit system) */ return SumInt(ProdInt(INTOBJ_INT(len), INTOBJ_INT(GMP_LIMB_BITS)),
INTOBJ_INT(a)); #endif
}
Obj IntStringInternal(Obj string, constChar *str)
{
Obj val; // value = <upp> * <pow> + <low>
Obj upp; // upper part Int pow; // power Int low; // lower part Int sign; // is the integer negative
UInt i; // loop variable
// if <string> is given, then we ignore <str> if (string)
str = CONST_CSTR_STRING(string);
// get the sign, if any
sign = 1;
i = 0; if (str[i] == '-') {
sign = -sign;
i++;
}
// collect the digits in groups of 8, for improved performance // note that 2^26 < 10^8 < 2^27, so the intermediate // values always fit into an immediate integer
low = 0;
pow = 1;
upp = INTOBJ_INT(0); while (str[i] != '\0') { if (str[i] < '0' || str[i] > '9') { return Fail;
}
low = 10 * low + str[i] - '0';
pow = 10 * pow; if (pow == 100000000L) {
upp = ProdInt(upp, INTOBJ_INT(pow));
upp = SumInt(upp, INTOBJ_INT(sign*low)); // refresh 'str', in case the arithmetic operations triggered // a garbage collection if (string)
str = CONST_CSTR_STRING(string);
pow = 1;
low = 0;
}
i++;
}
// check if 0 char does not mark the end of the string if (string && i < GET_LEN_STRING(string)) return Fail;
// compose the integer value if (upp == INTOBJ_INT(0)) {
val = INTOBJ_INT(sign * low);
} elseif (pow == 1) {
val = upp;
} else {
upp = ProdInt(upp, INTOBJ_INT(pow));
val = SumInt(upp, INTOBJ_INT(sign * low));
}
// if at least one input is a small int, a naive equality test suffices if (IS_INTOBJ(opL) || IS_INTOBJ(opR)) return opL == opR;
// compare the sign and size // note: at this point we know that opL and opR are proper bags, so // we can use TNUM_BAG instead of TNUM_OBJ if (TNUM_BAG(opL) != TNUM_BAG(opR) || SIZE_INT(opL) != SIZE_INT(opR)) return0;
/**************************************************************************** ** *FLtInt(<opL>,<opR>)......testifanintegerislessthananother ** **'LtInt'returns1iftheinteger<opL>isstrictlylessthanthe **integer<opR>and0otherwise.
*/ Int LtInt(Obj opL, Obj opR)
{ Int res;
CHECK_INT(opL);
CHECK_INT(opR);
// compare two small integers if (ARE_INTOBJS(opL, opR)) return (Int)opL < (Int)opR;
// a small int is always less than a positive large int, // and always more than a negative large int if (IS_INTOBJ(opL)) return IS_INTPOS(opR); if (IS_INTOBJ(opR)) return IS_INTNEG(opL);
// at this point, both inputs are large integers, and we compare their // signs first if (TNUM_OBJ(opL) != TNUM_OBJ(opR)) return IS_INTNEG(opL);
// signs are equal; compare sizes and absolute values if (SIZE_INT(opL) < SIZE_INT(opR))
res = -1; elseif (SIZE_INT(opL) > SIZE_INT(opR))
res = +1; else
res = mpn_cmp((mp_srcptr)CONST_ADDR_INT(opL),
(mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opL));
// if both arguments are negative, flip the result if (IS_INTNEG(opL))
res = -res;
// add or subtract if (sign == 1)
mpz_add( MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) ); else
mpz_sub( MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
// convert result to GAP object and return it
CHECK_FAKEMPZ(mpzResult);
CHECK_FAKEMPZ(mpzL);
CHECK_FAKEMPZ(mpzR);
result = GMPorINTOBJ_FAKEMPZ( mpzResult );
CHECK_INT(result); return result;
}
// convert result to GAP object and return it
CHECK_FAKEMPZ(mpzResult);
CHECK_FAKEMPZ(mpzL);
CHECK_FAKEMPZ(mpzR);
prd = GMPorINTOBJ_FAKEMPZ( mpzResult );
CHECK_INT(prd); return prd;
}
/**************************************************************************** ** *FProdIntObj(<n>,<op>)........productofanintegerandanobject
*/ static Obj ProdIntObj ( Obj n, Obj op )
{
Obj res = 0; // result
UInt i, k; // loop variables
mp_limb_t l; // loop variable
CHECK_INT(n);
// if the integer is zero, return the neutral element of the operand if ( n == INTOBJ_INT(0) ) {
res = ZERO_SAMEMUT( op );
}
/* if the integer is one, return the object if immutable - ifmutable,addtheobjecttoitsZeroSameMutabilityto
ensure correct mutability propagation */ elseif ( n == INTOBJ_INT(1) ) { if (IS_MUTABLE_OBJ(op))
res = SUM(ZERO_SAMEMUT(op),op); else
res = op;
}
// if the integer is minus one, return the inverse of the operand elseif ( n == INTOBJ_INT(-1) ) {
res = AINV_SAMEMUT( op );
}
// if the integer is negative, invert the operand and the integer elseif ( IS_NEG_INT(n) ) {
res = AINV_SAMEMUT( op ); if ( res == Fail ) {
ErrorMayQuit("Operations: <obj> must have an additive inverse", 0, 0);
}
res = PROD(AINV_SAMEMUT(n), res);
}
// if the integer is small, compute the product by repeated doubling // the loop invariant is <result> = <k>*<res> + <l>*<op>, <l> < <k> // <res> = 0 means that <res> is the neutral element elseif ( IS_INTOBJ(n) && INT_INTOBJ(n) > 1 ) {
res = 0;
k = (Int)1 << NR_SMALL_INT_BITS;
l = INT_INTOBJ(n); while ( 0 < k ) {
res = (res == 0 ? res : SUM( res, res )); if ( k <= l ) {
res = (res == 0 ? op : SUM( res, op ));
l = l - k;
}
k = k / 2;
}
}
// if the integer is large, compute the product by repeated doubling elseif ( IS_INTPOS(n) ) {
res = 0; for ( i = SIZE_INT(n); 0 < i; i-- ) {
k = 8*sizeof(mp_limb_t);
l = CONST_ADDR_INT(n)[i-1]; while ( 0 < k ) {
res = (res == 0 ? res : SUM( res, res ));
k--; if ( (l >> k) & 1 ) {
res = (res == 0 ? op : SUM( res, op ));
}
}
}
}
// power with a large exponent elseif ( ! IS_INTOBJ(opR) ) {
ErrorMayQuit("Integer operands: <exponent> is too large", 0, 0);
}
// power with a negative exponent elseif ( INT_INTOBJ(opR) < 0 ) {
pow = QUO( INTOBJ_INT(1),
PowInt( opL, INTOBJ_INT( -INT_INTOBJ(opR)) ) );
}
// findme - can we use the gmp function mpz_n_pow_ui?
// power with a small positive exponent, do it by a repeated squaring else {
pow = INTOBJ_INT(1);
i = INT_INTOBJ(opR); while ( i != 0 ) { if ( i % 2 == 1 ) pow = ProdInt( pow, opL ); if ( i > 1 ) opL = ProdInt( opL, opL );
TakeInterrupt();
i = i / 2;
}
}
// return the power
CHECK_INT(pow); return pow;
}
/**************************************************************************** ** *FPowObjInt(<op>,<n>)..........powerofanobjectandaninteger
*/ static Obj PowObjInt(Obj op, Obj n)
{
Obj res = 0; // result
UInt i, k; // loop variables
mp_limb_t l; // loop variable
CHECK_INT(n);
// if the integer is zero, return the neutral element of the operand if ( n == INTOBJ_INT(0) ) { return ONE_SAMEMUT( op );
}
// if the integer is one, return a copy of the operand elseif ( n == INTOBJ_INT(1) ) {
res = CopyObj( op, 1 );
}
// if the integer is minus one, return the inverse of the operand elseif ( n == INTOBJ_INT(-1) ) {
res = INV_SAMEMUT( op );
}
// if the integer is negative, invert the operand and the integer elseif ( IS_NEG_INT(n) ) {
res = INV_SAMEMUT( op ); if ( res == Fail ) {
ErrorMayQuit("Operations: <obj> must have an inverse", 0, 0);
}
res = POW(res, AINV_SAMEMUT(n));
}
// if the integer is small, compute the power by repeated squaring // the loop invariant is <result> = <res>^<k> * <op>^<l>, <l> < <k> // <res> = 0 means that <res> is the neutral element elseif ( IS_INTOBJ(n) && INT_INTOBJ(n) > 0 ) {
res = 0;
k = (Int)1 << NR_SMALL_INT_BITS;
l = INT_INTOBJ(n); while ( 0 < k ) {
res = (res == 0 ? res : PROD( res, res )); if ( k <= l ) {
res = (res == 0 ? op : PROD( res, op ));
l = l - k;
}
k = k / 2;
}
}
// if the integer is large, compute the power by repeated squaring elseif ( IS_INTPOS(n) ) {
res = 0; for ( i = SIZE_INT(n); 0 < i; i-- ) {
k = 8*sizeof(mp_limb_t);
l = CONST_ADDR_INT(n)[i-1]; while ( 0 < k ) {
res = (res == 0 ? res : PROD( res, res ));
k--; if ( (l >> k) & 1 ) {
res = (res == 0 ? op : PROD( res, op ));
}
}
}
}
/**************************************************************************** ** *FModInt(<opL>,<opR>)..representativeofresidueclassofaninteger ** **'ModInt'returnsthesmallestpositiverepresentativeoftheresidue **classoftheinteger<opL>modulotheinteger<opR>.
*/
Obj ModInt(Obj opL, Obj opR)
{ Int i; // loop count, value for small int Int k; // loop count, value for small int
UInt c; // product of two digits
Obj mod; // handle of the remainder bag
Obj quo; // handle of the quotient bag
CHECK_INT(opL);
CHECK_INT(opR);
// pathological case first
RequireNonzero("Integer operations", opR, "divisor");
// compute the remainder of two small integers if ( ARE_INTOBJS( opL, opR ) ) {
// get the integer values
i = INT_INTOBJ(opL);
k = INT_INTOBJ(opR);
// compute the remainder, make sure we divide only positive numbers
i %= k; if (i < 0)
i += k > 0 ? k : -k;
mod = INTOBJ_INT(i);
}
// compute the remainder of a small integer by a large integer elseif ( IS_INTOBJ(opL) ) {
// the small int -(1<<28) mod the large int (1<<28) is 0 if ( opL == INTOBJ_MIN
&& ( IS_INTPOS(opR) )
&& ( SIZE_INT(opR) == 1 )
&& ( VAL_LIMB0(opR) == -INT_INTOBJ_MIN ) )
mod = INTOBJ_INT(0);
// in all other cases the remainder is equal the left operand elseif ( 0 <= INT_INTOBJ(opL) )
mod = opL; elseif ( IS_INTPOS(opR) )
mod = SumOrDiffInt( opL, opR, 1 ); else
mod = SumOrDiffInt( opL, opR, -1 );
}
// compute the remainder of a large integer by a small integer elseif ( IS_INTOBJ(opR) ) {
// get the integer value, make positive
i = INT_INTOBJ(opR); if ( i < 0 ) i = -i;
// check whether right operand is a small power of 2 if ( !(i & (i-1)) ) {
c = VAL_LIMB0(opL) & (i-1);
}
// otherwise use the gmp function to divide else {
c = mpn_mod_1( (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL), i );
}
// now c is the absolute value of the actual result. Thus, if the left // operand is negative, and c is non-zero, we have to adjust it. if (IS_INTPOS(opL) || c == 0)
mod = INTOBJ_INT( c ); else // even if opR is INT_INTOBJ_MIN, and hence i is INT_INTOBJ_MAX+1, we // have 0 <= i-c <= INT_INTOBJ_MAX, so i-c fits into a small integer
mod = INTOBJ_INT( i - (Int)c );
}
// compute the remainder of a large integer modulo a large integer else {
// trivial case first if ( SIZE_INT(opL) < SIZE_INT(opR) ) { if ( IS_INTPOS(opL) ) return opL; elseif ( IS_INTPOS(opR) )
mod = SumOrDiffInt( opL, opR, 1 ); else
mod = SumOrDiffInt( opL, opR, -1 ); #if DEBUG_GMP
assert( !IS_NEG_INT(mod) ); #endif
CHECK_INT(mod); return mod;
}
mod = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+1)*sizeof(mp_limb_t) );
ENSURE_BAG(mod);
quo = NewBag( T_INTPOS,
(SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) );
ENSURE_BAG(quo);
// and let gmp do the work
mpn_tdiv_qr( (mp_ptr)ADDR_INT(quo), (mp_ptr)ADDR_INT(mod), 0,
(mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL),
(mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opR) );
// reduce to small integer if possible, otherwise shrink bag
mod = GMP_NORMALIZE( mod );
// make the representative positive if ( IS_NEG_INT(mod) ) { if ( IS_INTPOS(opR) )
mod = SumOrDiffInt( mod, opR, 1 ); else
mod = SumOrDiffInt( mod, opR, -1 );
}
/**************************************************************************** ** *FQuoInt(<opL>,<opR>).............quotientoftwointegers ** **'QuoInt'returnstheintegerpartofthetwointegers<opL>and<opR>. ** **Notethatthisroutineisnotcalledfrom'EvalQuo',thedivisionoftwo **integersyieldsarationalandisthereforeperformedin'QuoRat'.This **operationishoweveravailablethroughtheinternalfunction'Quo'.
*/
Obj QuoInt(Obj opL, Obj opR)
{ Int i; // loop count, value for small int Int k; // loop count, value for small int
Obj quo; // handle of the result bag
Obj rem; // handle of the remainder bag
CHECK_INT(opL);
CHECK_INT(opR);
// pathological case first
RequireNonzero("Integer operations", opR, "divisor");
// divide two small integers if ( ARE_INTOBJS( opL, opR ) ) {
// the small int -(1<<28) divided by -1 is the large int (1<<28) if ( opL == INTOBJ_MIN && opR == INTOBJ_INT(-1) ) {
quo = NewBag( T_INTPOS, sizeof(mp_limb_t) );
SET_VAL_LIMB0( quo, -INT_INTOBJ_MIN ); return quo;
}
// get the integer values
i = INT_INTOBJ(opL);
k = INT_INTOBJ(opR);
// divide, truncated towards zero; this is also what section 6.5.5 of the // C99 standard guarantees for integer division
quo = INTOBJ_INT(i / k);
}
// divide a small integer by a large one elseif ( IS_INTOBJ(opL) ) {
// the small int -(1<<28) divided by the large int (1<<28) is -1 if ( opL == INTOBJ_MIN
&& IS_INTPOS(opR) && SIZE_INT(opR) == 1
&& VAL_LIMB0(opR) == -INT_INTOBJ_MIN )
quo = INTOBJ_INT(-1);
// in all other cases the quotient is of course zero else
quo = INTOBJ_INT(0);
}
// divide a large integer by a small integer elseif ( IS_INTOBJ(opR) ) {
k = INT_INTOBJ(opR);
// allocate a bag for the result and set up the pointers if ( IS_INTNEG(opL) == ( k < 0 ) )
quo = NewBag( T_INTPOS, SIZE_OBJ(opL) ); else
quo = NewBag( T_INTNEG, SIZE_OBJ(opL) );
ENSURE_BAG(quo);
if ( k < 0 ) k = -k;
// use gmp function for dividing by a 1-limb number
mpn_divrem_1( (mp_ptr)ADDR_INT(quo), 0,
(mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL),
k );
}
// divide a large integer by a large integer else {
// trivial case first if ( SIZE_INT(opL) < SIZE_INT(opR) ) return INTOBJ_INT(0);
// create a new bag for the remainder
rem = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+1)*sizeof(mp_limb_t) );
ENSURE_BAG(rem);
// allocate a bag for the quotient if ( TNUM_OBJ(opL) == TNUM_OBJ(opR) )
quo = NewBag( T_INTPOS,
(SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) ); else
quo = NewBag( T_INTNEG,
(SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) );
ENSURE_BAG(quo);
/**************************************************************************** ** *FRemInt(<opL>,<opR>)............remainderoftwointegers ** **'RemInt'returnstheremainderofthequotientoftheintegers<opL> **and<opR>. ** **Notethattheremainderisdifferentfromthevaluereturnedbythe'mod' **operatorwhichisalwayspositive,whiletheresultof'RemInt'has **thesamesignas<opL>.
*/
Obj RemInt(Obj opL, Obj opR)
{ Int i; // loop count, value for small int Int k; // loop count, value for small int
UInt c;
Obj rem; // handle of the remainder bag
Obj quo; // handle of the quotient bag
CHECK_INT(opL);
CHECK_INT(opR);
// pathological case first
RequireNonzero("Integer operations", opR, "divisor");
// compute the remainder of two small integers if ( ARE_INTOBJS( opL, opR ) ) {
// get the integer values
i = INT_INTOBJ(opL);
k = INT_INTOBJ(opR);
// compute the remainder with sign matching that of the dividend i; // this matches the sign of i % k as specified in section 6.5.5 of // the C99 standard, which indicates division truncates towards zero, // and the invariant i%k == i - (i/k)*k holds
rem = INTOBJ_INT(i % k);
}
// compute the remainder of a small integer by a large integer elseif ( IS_INTOBJ(opL) ) {
// the small int -(1<<28) rem the large int (1<<28) is 0 if ( opL == INTOBJ_MIN
&& IS_INTPOS(opR) && SIZE_INT(opR) == 1
&& VAL_LIMB0(opR) == -INT_INTOBJ_MIN )
rem = INTOBJ_INT(0);
// in all other cases the remainder is equal the left operand else
rem = opL;
}
// compute the remainder of a large integer by a small integer elseif ( IS_INTOBJ(opR) ) {
// get the integer value, make positive
i = INT_INTOBJ(opR); if ( i < 0 ) i = -i;
// check whether right operand is a small power of 2 if ( !(i & (i-1)) ) {
c = VAL_LIMB0(opL) & (i-1);
}
// otherwise use the gmp function to divide else {
c = mpn_mod_1( (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL), i );
}
// adjust c for the sign of the left operand if ( IS_INTPOS(opL) )
rem = INTOBJ_INT( c ); else
rem = INTOBJ_INT( -(Int)c );
}
// compute the remainder of a large integer modulo a large integer else {
// trivial case first if ( SIZE_INT(opL) < SIZE_INT(opR) ) return opL;
rem = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+1)*sizeof(mp_limb_t) );
ENSURE_BAG(rem);
quo = NewBag( T_INTPOS,
(SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) );
ENSURE_BAG(quo);
// and let gmp do the work
mpn_tdiv_qr( (mp_ptr)ADDR_INT(quo), (mp_ptr)ADDR_INT(rem), 0,
(mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL),
(mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opR) );
// reduce to small integer if possible, otherwise shrink bag
rem = GMP_NORMALIZE( rem );
}
// compute the gcd
mpz_gcd( MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
// convert result to GAP object and return it
CHECK_FAKEMPZ(mpzResult);
CHECK_FAKEMPZ(mpzL);
CHECK_FAKEMPZ(mpzR);
result = GMPorINTOBJ_FAKEMPZ( mpzResult );
CHECK_INT(result); return result;
}
// compute the gcd
mpz_lcm(MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR));
// convert result to GAP object and return it
CHECK_FAKEMPZ(mpzResult);
CHECK_FAKEMPZ(mpzL);
CHECK_FAKEMPZ(mpzR);
result = GMPorINTOBJ_FAKEMPZ(mpzResult);
CHECK_INT(result); return result;
}
static Obj FuncLCM_INT(Obj self, Obj a, Obj b)
{
RequireInt(SELF_NAME, a);
RequireInt(SELF_NAME, b); return LcmInt(a, b);
}
// convert mpzResult into a GAP integer object.
Obj result = GMPorINTOBJ_MPZ(mpzResult);
// free mpzResult
mpz_clear(mpzResult);
return result;
}
/**************************************************************************** **
*/
Obj BinomialInt(Obj n, Obj k)
{ Int negate_result = 0;
// deal with k <= 1 if (k == INTOBJ_INT(0)) return INTOBJ_INT(1); if (k == INTOBJ_INT(1)) return n; if (IS_NEG_INT(k)) return INTOBJ_INT(0);
// deal with n < 0 if (IS_NEG_INT(n)) { // use the identity Binomial(n,k) = (-1)^k * Binomial(-n+k-1, k)
negate_result = IS_ODD_INT(k);
n = DiffInt(DiffInt(k,n), INTOBJ_INT(1));
}
// deal with n <= k if (n == k) return negate_result ? INTOBJ_INT(-1) : INTOBJ_INT(1); if (LtInt(n, k)) return INTOBJ_INT(0);
// deal with n-k < k <=> n < 2k
Obj k2 = DiffInt(n, k); if (LtInt(k2, k))
k = k2;
// From here on, we only support single limb integers for k. Anything else // would lead to output too big for storage anyway, at least on a 64 bit // system. To be specific, in that case n >= 2k and k >= 2^60. Thus, in // the *best* case, we are trying to compute the central binomial // coefficient Binomial(2k k) for k = 2^60. This value is approximately // (4^k / sqrt(pi*k)); taking the logarithm and dividing by 8 yields that // we need about k/4 = 2^58 bytes to store this value. No computer in the // foreseeable future will be able to store the result (and much less // compute it in a reasonable time). // // On 32 bit systems, the limit is k = 2^28, and then the central binomial // coefficient Binomial(2k,k) takes up about 64 MB, so that would still be // feasible. However, GMP does not support computing binomials when k is // larger than a single limb, so we'd have to implement this on our own. // // Since GAP previously was effectively unable to compute such binomial // coefficients (unless you were willing to wait for a few days or so), we // simply do not implement this on 32bit systems, and instead return Fail, // jut as we do on 64 bit. If somebody complains about this, we can still // look into implementing this (and in the meantime, tell the user to use // the old GAP version of this function).
if (SIZE_INT_OR_INTOBJ(n) == 1 && SIZE_INT_OR_INTOBJ(p) == 1) {
UInt N = AbsOfSmallInt(n);
UInt P = AbsOfSmallInt(p); if (N == 0 || P == 1) return INTOBJ_INT(0);
k = 0; while (N % P == 0) {
N /= P;
k++;
} return INTOBJ_INT(k);
}
/* For certain values of p, mpz_remove replaces its "dest" argument andtriestodeallocatetheoriginalmpz_tinit.Thismeans wecannotuseafake_mpz_tforit.However,wearenotreally
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 );
if (!IS_POS_INT(k))
ErrorMayQuit("Root: <k> must be a positive integer", 0, 0); if (IS_NEG_INT(n) && IS_EVEN_INT(k))
ErrorMayQuit("Root: <n> is negative but <k> is even", 0, 0);
if (k == INTOBJ_INT(1) || n == INTOBJ_INT(0)) return n;
if (!IS_INTOBJ(k)) { #ifdef SYS_IS_64_BIT // if k is not immediate, i.e., k >= 2^60, then the root is 1 unless // n >= 2^k >= 2^(2^60), which means storage for n must be at least // 2^60 bits, i.e. 2^57 bytes. That's more RAM than anybody has for // the near future. return IS_NEG_INT(n) ? INTOBJ_INT(-1) : INTOBJ_INT(1); #else // if k is not immediate, i.e., k >= 2^28, then the root is 1 unless // n >= 2^k >= 2^(2^28), which means storage for n must be at least // 2^28 bits, i.e., 2^25 bytes, or 2^23 words (each 32 bits). if (SIZE_INT_OR_INTOBJ(n) < (1 << 23)) return IS_NEG_INT(n) ? INTOBJ_INT(-1) : INTOBJ_INT(1); else return Fail; // return fail so that high level code can handle // this #endif
}
if (mod == INTOBJ_INT(1) || mod == INTOBJ_INT(-1)) return INTOBJ_INT(0); if (base == INTOBJ_INT(0) || mod == INTOBJ_INT(0)) return Fail;
// handle small inputs separately if (IS_INTOBJ(mod)) {
Int a = INT_INTOBJ(mod); if (a < 0)
a = -a;
Int b = INT_INTOBJ(ModInt(base, mod));
Int aL = 0; // cofactor of a Int bL = 1; // cofactor of b
// 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);
}
RequireNonzero(SELF_NAME, mod, "mod"); if ( mod == INTOBJ_INT(1) || mod == INTOBJ_INT(-1) ) return INTOBJ_INT(0);
if ( IS_NEG_INT(exp) ) {
base = InverseModInt( base, mod ); if (base == Fail)
ErrorMayQuit("PowerModInt: negative <exp> but <base> is not invertible modulo <mod>", 0, 0);
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 );
/**************************************************************************** ** *FRandomIntegerMT(<mtstr>,<nrbits>) ** **Returnsanintegerwithatmost<nrbits>bitsinuniformdistribution. **<nrbits>mustbeasmallinteger.<mtstr>isastringasreturnedby **InitRandomMT. ** **Implementationdetailsareabittrickytoobtainthesamerandom **integerson32bitand64bitmachines(whichhavedifferentrangesof **integers). **
*/ static Obj FuncRandomIntegerMT(Obj self, Obj mtstr, Obj nrbits)
{
Obj res; Int i, n, q, r, qoff, len;
UInt4 *mt;
UInt4 *pt;
RequireStringRep(SELF_NAME, mtstr); if (GET_LEN_STRING(mtstr) < 2500) {
ErrorMayQuit( "RandomIntegerMT: <mtstr> must be a string with at least 2500 characters", 0, 0);
}
RequireNonnegativeSmallInt(SELF_NAME, nrbits);
n = INT_INTOBJ(nrbits);
// small int case if (n <= NR_SMALL_INT_BITS) {
mt = (UInt4 *)(ADDR_OBJ(mtstr) + 1); #ifdef SYS_IS_64_BIT if (n <= 32) {
res = INTOBJ_INT((Int)(nextrandMT_int32(mt) & ((UInt4)-1 >> (32-n))));
} else {
UInt8 rd;
rd = nextrandMT_int32(mt);
rd += (UInt8) ((UInt4) nextrandMT_int32(mt) &
((UInt4)-1 >> (64-n))) << 32;
res = INTOBJ_INT((Int)rd);
} #else
res = INTOBJ_INT((Int)(nextrandMT_int32(mt) & ((UInt4)-1 >> (32-n)))); #endif
} else { // large int case
q = n / 32;
r = n - q * 32; // qoff = number of 32 bit words we need
qoff = q + (r==0 ? 0:1); // len = number of limbs we need (limbs currently are either 32 or 64 bit wide)
len = (qoff*4 + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
res = NewBag( T_INTPOS, len*sizeof(mp_limb_t) );
pt = (UInt4*) ADDR_INT(res);
mt = (UInt4 *)(ADDR_OBJ(mtstr) + 1); for (i = 0; i < qoff; i++, pt++) {
*pt = nextrandMT_int32(mt);
} if (r != 0) { // we generated too many random bits -- chop of the extra bits
pt = (UInt4*) ADDR_INT(res);
pt[qoff-1] = pt[qoff-1] & ((UInt4)(-1) >> (32-r));
} #ifdefined(SYS_IS_64_BIT) && defined(WORDS_BIGENDIAN) // swap the halves of the 64bit words to match the // little endian resp. 32 bit versions of this code
pt = (UInt4 *)ADDR_INT(res); for (i = 0; i < qoff; i += 2, pt += 2) {
SWAP(UInt4, pt[0], pt[1]);
} #endif
res = GMP_NORMALIZE(res);
}
/**************************************************************************** ** *FInitInfoInt()..................tableofinitfunctions
*/ static StructInitInfo module = { // init struct using C99 designated initializers; for a full list of // fields, please refer to the definition of StructInitInfo
.type = MODULE_BUILTIN,
.name = "integer",
.initKernel = InitKernel,
.initLibrary = InitLibrary,
};
¤ Diese beiden folgenden Angebotsgruppen bietet das Unternehmen0.82Angebot
(Wie Sie bei der Firma Beratungs- und Dienstleistungen beauftragen können 2026-09-28)
¤
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung und die Messung sind noch experimentell.