############################################################################# ## ## This file is part of GAP, a system for computational discrete algebra. ## This file's authors include Thomas Breuer, Frank Celler, Stefan Kohl, Werner Nickel, Alice Niemeyer, Martin Schönert, Alex Wegner. ## ## Copyright of GAP belongs to its developers, whose names are too numerous ## to list here. Please refer to the COPYRIGHT file for details. ## ## SPDX-License-Identifier: GPL-2.0-or-later ##
InstallMethod( Coefficients, "for the canonical basis of Integers",
IsCollsElms,
[ IsBasis and IsCanonicalBasis and IsCanonicalBasisIntegersRep,
IsCyc ], 0, function( B, v ) if IsInt( v ) then
return [ v ]; else
return fail; fi; end );
if IsHPCGAP then
ShareSpecialObj(Primes2);
ShareSpecialObj(ProbablePrimes2); fi;
############################################################################# ## #F BestQuoInt( <n>, <m> ) ## ## `BestQuoInt' returns the best quotient <q> of the integers <n> and <m>. ## This is the quotient such that `<n>-<q>\*<m>' has minimal absolute value. ## If there are two quotients whose remainders have the same absolute value, ## then the quotient with the smaller absolute value is chosen. ##
InstallGlobalFunction(BestQuoInt,function ( n, m ) if0 <= m and 0 <= n then
return QuoInt( n + QuoInt( m - 1, 2 ), m ); elif0 <= m then
return QuoInt( n - QuoInt( m - 1, 2 ), m ); elif0 <= n then
return QuoInt( n - QuoInt( m + 1, 2 ), m ); else
return QuoInt( n + QuoInt( m + 1, 2 ), m ); fi; end);
############################################################################# ## #F ChineseRem( <moduli>, <residues> ) . . . . . . . . . . chinese remainder ##
InstallGlobalFunction(ChineseRem,function ( moduli, residues )
local i, c, l, g;
# combine the residues modulo the moduli
i := 1;
c := residues[1];
l := moduli[1]; while i < Length(moduli) do
i := i + 1;
g := Gcdex( l, moduli[i] ); if g.gcd <> 1 and (residues[i]-c) mod g.gcd <> 0then
Error("the residues must be equal modulo ",g.gcd); fi;
c := l * (((residues[i]-c) / g.gcd * g.coeff1) mod moduli[i]) + c;
l := moduli[i] / g.gcd * l;
od;
# reduce c into the range [0..l-1]
c := c mod l;
return c; end);
############################################################################# ## #F CoefficientsQadic( <i>, <q> ) . . . . . . <q>-adic representation of <i> ##
InstallMethod( CoefficientsQadic, "for two integers",
true, [ IsInt, IsInt ], 0, function( i, q )
local i1, res, l, qq, i2; if q <= 1then
Error("2nd argument of CoefficientsQadic should be greater than 1\n"); fi; if i < 0then # if FR package is loaded and supplies an implementation # to return a periodic list for negative i
TryNextMethod(); fi; # represent the integer <i> as <q>-adic number if i = 0then
return []; elif i < q then
return [i]; elif i < q^2then
i1 := QuoInt(i, q);
return [i - i1*q, i1]; elif Log2Int(q)*100 > Log2Int(i) then # straight forward loop for result length < 100
res := []; while i > 0do
i1 := QuoInt(i, q);
Add(res, i - i1*q);
i := i1;
od; else # divide and conquer method for large i
l := QuoInt(LogInt(i, q), 2)+1;
qq := q^l;
i2 := QuoInt(i,qq);
i1 := i - i2*qq;
res := CoefficientsQadic(i1, q); while Length(res) < l do
Add(res, 0);
od;
Append(res, CoefficientsQadic(i2, q)); fi;
return res; end);
############################################################################# ## #F CoefficientsMultiadic( ints, int ) ##
InstallGlobalFunction(CoefficientsMultiadic, function( ints, int )
local vec, i;
vec := List( ints, x -> 0 ); for i in Reversed( [1..Length(ints)] ) do
vec[i] := RemInt( int, ints[i] );
int := QuoInt( int, ints[i] );
od;
return vec; end);
############################################################################# ## #F DivisorsInt( <n> ) . . . . . . . . . . . . . . . divisors of an integer ##
BindGlobal("DivisorsIntCache",
List([[1],[1,2],[1,3],[1,2,4],[1,5],[1,2,3,6],[1,7]], Immutable)); if IsHPCGAP then
MakeImmutable(DivisorsIntCache); fi;
InstallGlobalFunction(DivisorsInt,function ( n )
local divisors, factors, divs;
# make <n> it nonnegative, handle trivial cases, and get prime factors if n < 0then n := -n; fi; if n = 0then Error("DivisorsInt: <n> must not be 0"); fi; if n <= Length(DivisorsIntCache) then
return DivisorsIntCache[n]; fi;
factors := Factors(Integers, n );
# recursive function to compute the divisors
divs := function ( i, m ) if Length(factors) < i then return [ m ]; elif m mod factors[i] = 0then return divs(i+1,m*factors[i]); else return Concatenation( divs(i+1,m), divs(i+1,m*factors[i]) ); fi; end;
############################################################################# ## #F FactorsRho( <n>, <inc>, <cluster>, <limit> ) Pollards rho factorization ## ## `FactorsInt' does trial divisions by the primes less than 1000 to detect ## all composites with a factor less than 1000 and primes less than 1000000. ## After that it calls `FactorsRho(<n>,1,16,8192)' to do the hard work. ## ## `FactorsRho' will return a list of factors and a list of composite ## number. Usually `FactorsInt' factors integers with prime factors ## $\<1000$ faster. However for integers with no factor $\<1000$ ## `FactorsRho' will be faster. ## ## `FactorsRho' uses Pollards $\rho$ method to factor the integer $n = p q$. ## For a small simple example lets assume we want to factor $667 = 23 * 29$. ## `FactorsRho' first calls `IsPrimeInt' to avoid trying to factor a prime. ## ## Then it uses the sequence defined by $x_0=1, x_{i+1}=(x_i^2+1)$ mod $n$. ## In our example this is $1, 2, 5, 26, 10, 101, 197, 124, 36, 630, .. $. ## ## Modulo $p$ it takes on at most $p-1$ different values, thus it eventually ## becomes recurrent, usually this happens after roughly $2 \sqrt{p}$ steps. ## In our example modulo 23 we get $1, 2, 5, 3, 10, 9, 13, 9, 13, 9, .. $. ## ## Thus there exist pairs $i, j$ such that $x_i = x_j$ mod $p$, i.e., such ## that $p$ divides $Gcd( n, x_j-x_i )$. With a bit of luck no other factor ## of $n$ divides $x_j - x_i$ so we find $p$ if we know such a pair. In our ## example $5, 7$ is the first pair, $x_7-x_5=23$, and $Gcd(667,23) = 23$. ## ## Now it is too expensive to check all pairs, but there also must be pairs ## of the form $2^i-1, j$ with $3*2^{i-1} <= j < 4*2^{i-1}$. In our example ## $7, 13$ is the first such pair, $x_13-x_7=506$, and $Gcd(667,506) = 23$. ## ## Thus by taking the gcds of $n$ and $x_j-x_i$ for such pairs, we will find ## the factor $p$ after approximately $2 \sqrt{p} \<= 2 \sqrt^4{n}$ steps. ## ## If $Gcd( n, x_j - x_i )$ is not a prime `FactorsRho' will call itself ## recursively with a different value for <inc>, i.e., it will try to factor ## the gcd using a different sequence $x_{i+1} = (x_i^2 + inc)$ mod $n$. ## ## Since the gcd computations are by far the most time consuming part of the ## algorithm one can save time by clustering differences and computing the ## gcd only every <cluster> iteration. This slightly increases the chance ## that a gcd is composite, but reduces the runtime by a large amount. ## ## Finally `FactorsRho' accepts an argument <limit> which is the number of ## iterations performed by `FactorsRho' before giving up. The default value ## is 8192 which corresponds to a few minutes while guaranteeing that all ## prime factors less than $10^6$ and most less than $10^9$ are found. ## ## Better descriptions of the algorithm and related topics can be found in: ## J. Pollard, A Monte Carlo Method for Factorization, BIT 15, 1975, 331-334 ## R. Brent, An Improved Monte Carlo Method for Fact., BIT 20, 1980, 176-184 ## D. Knuth, Seminumerical Algorithms (TACP II), AddiWesl, 1973, 369-371 ##
DeclareGlobalName( "FactorsRho" );
BindGlobal( "FactorsRho", function ( n, inc, cluster, limit )
local i, sign, factors, composite, x, y, k, z, g, tmp,
IsPrimeOrProbablyPrimeInt;
# make $n$ positive and handle trivial cases
sign := 1; if n < 0then sign := -sign; n := -n; fi; if n < 4then return [ [ sign * n ], [] ]; fi;
factors := [];
composite := []; while n mod 2 = 0do Add( factors, 2 ); n := n / 2; od; while n mod 3 = 0do Add( factors, 3 ); n := n / 3; od;
if ValueOption("UseProbabilisticPrimalityTest") = true then IsPrimeOrProbablyPrimeInt := IsProbablyPrimeInt; else IsPrimeOrProbablyPrimeInt := IsPrimeInt; fi;
if IsPrimeOrProbablyPrimeInt(n) then Add( factors, n ); n := 1; fi;
# initialize $x_0$
x := 1; z := 1; i := 0;
# loop until we have factored $n$ completely or run out of patience while1 < n and 2^i <= limit do
# $y = x_{2^i-1}$
y := x; i := i + 1;
# $x_{2^i}, .., x_{3*2^{i-1}-1}$ need not be compared to $x_{2^i-1}$ for k in [1..2^(i-1)] do
x := (x^2 + inc) mod n;
od;
# compare $x_{3*2^{i-1}}, .., x_{4*2^{i-1}-1}$ with $x_{2^i-1}$ for k in [1..2^(i-1)] do
x := (x^2 + inc) mod n;
z := z * (x - y) mod n;
# from time to time compute the gcd if k mod cluster = 0then
g := GcdInt( n, z );
# if it is > 1 we have found a factor which need not be prime if g > 1then
tmp := FactorsRho(g,inc+1,QuoInt(cluster+1,2),limit);
factors := Concatenation( factors, tmp[1] );
composite := Concatenation( composite, tmp[2] );
n := n / g; if IsPrimeOrProbablyPrimeInt(n) then
Add( factors, n ); n := 1; fi; fi; fi;
od;
od;
# add <n> to the list of composite numbers if1 < n then
Add( composite, n ); fi;
# sort the list of factors and composite numbers and return it
Sort(factors);
Sort(composite); if0 < Length(factors) then
factors[1] := sign * factors[1]; else
composite[1] := sign * composite[1]; fi;
return [ factors, composite ];
end );
############################################################################# ## #F FactorsInt( <n> ) . . . . . . . . . . . . . . prime factors of an integer #F FactorsInt( <n> : RhoTrials := <trials>) #F FactorsInt( <n> : quiet) ## ## In the second form, FactorsRho is called with a limit of <trials> ## on the number of trials is performs. The default is 8192. ## ## The option `quiet' makes the function return even if the `rho' ## factorization failed and return the factorization found so far. ##
InstallGlobalFunction(FactorsInt,function ( n )
local sign, factors, p, tmp, n_orig, len, rt, tmp2;
n_orig := n;
# make $n$ positive and handle trivial cases
sign := 1; if n < 0then sign := -sign; n := -n; fi; if n < 4then return [ sign * n ]; fi;
factors := [];
# do trial divisions by the primes less than 1000 # faster than anything fancier because $n$ mod <small int> is very fast for p in Primes do while n mod p = 0do Add( factors, p ); n := n / p; od; if n < (p+1)^2 and 1 < n then Add(factors,n); n := 1; fi; if n = 1then factors[1] := sign*factors[1]; return factors; fi;
od;
# do trial divisions by known primes
atomic readonly Primes2 do for p in Primes2 do while n mod p = 0do Add( factors, p ); n := n / p; od; if p^2 > n then break; fi; if n = 1then factors[1] := sign*factors[1]; return factors; fi;
od;
od;
# do trial divisions by known probable primes (and issue warning, if found)
tmp := [];
atomic readonly ProbablePrimes2 do for p in ProbablePrimes2 do while n mod p = 0do
AddSet(tmp, p);
Add( factors, p );
n := n / p;
od; if n = 1then break; fi;
od;
od; if Length(tmp) > 0then
Info(InfoPrimeInt, 1 , "FactorsInt: used the following factor(s) which are probably primes:"); for p in tmp do
Info(InfoPrimeInt, 1, " ", p);
od; fi; if n = 1then factors[1] := sign*factors[1]; return factors; fi;
# handle perfect powers
p := SmallestRootInt( n ); if p < n then while1 < n do
Append( factors, FactorsInt(p) );
n := n / p;
od;
Sort( factors );
factors[1] := sign * factors[1];
return factors; fi;
# let `FactorsRho' do the work if ValueOption("RhoTrials") <> fail then
tmp := FactorsRho( n, 1, 16, ValueOption("RhoTrials") ); else
tmp := FactorsRho( n, 1, 16, 8192 ); fi; if0 < Length(tmp[2]) then if ValueOption("quiet")<>true then
len := Length(tmp[2]); if IsPackageMarkedForLoading("FactInt", "") then ## # in general cases we should proceed with the found factors: ## while len > 0 do ## Append(tmp[1], Factors(tmp[2][len])); ## Unbind(tmp[2][len]); ## len := len-1; ## od; # but this way we miss that FactInt can detect certain numbers of # special shape for which it uses lookup tables, therefore for the # moment:
return Factors(n_orig); else
Error( "sorry, cannot factor ", tmp[2], "\ntype 'return;' to try again with a larger number of trials in\n", "FactorsRho (or use option 'RhoTrials')\n"); if ValueOption("RhoTrials") <> fail then
rt := 5 * ValueOption("RhoTrials"); else
rt := 5 * 8192; fi; while len > 0do
tmp2 := FactorsInt(tmp[2][len]: RhoTrials := rt);
Append(tmp[1], tmp2);
Unbind(tmp[2][len]);
len := len-1;
od; fi; else
factors := Concatenation( factors, tmp[2] ); fi; fi;
factors := Concatenation( factors, tmp[1] );
Sort( factors );
factors[1] := sign * factors[1];
return factors; end);
############################################################################# ## #F PrimeDivisors( <n> ) . . . . . . . . . . . . . . list of prime divisors ## ## delegating to Factors ##
InstallMethod( PrimeDivisors, "for integer", [ IsInt ], function(n) if n = 0then
Error( "<n> must be non zero" ); fi; if n < 0then
n := -n; fi; if n = 1then
return []; fi;
return Set(Factors(Integers,n)); end);
local N, sign, factors, p, k, root, rootfactors, rhotrials,
tmp, CheckAndSortFactors;
CheckAndSortFactors := function ( )
factors := SortedList(factors);
factors[1] := sign*factors[1]; if Product(factors) <> N then Error("PartialFactorization: Internal error, wrong result!"); fi; end;
N := n; if effort < 0then effort := 5; fi;
# make $n$ positive and handle trivial cases
sign := 1; if n < 0then sign := -sign; n := -n; fi; if n < 4then return [ sign * n ]; fi;
factors := [];
# least effort: do trial divisions by the primes less than 100 if effort = 0then for p in Primes{[1..25]} do while n mod p = 0do Add( factors, p ); n := n / p; od; if n < (p+1)^2 and 1 < n then Add(factors,n); n := 1; fi; if n = 1then CheckAndSortFactors(); return factors; fi;
od;
Add(factors,n); CheckAndSortFactors(); return factors; fi;
# do trial divisions by the primes less than 1000 # faster than anything fancier because $n$ mod <small int> is very fast for p in Primes do while n mod p = 0do Add( factors, p ); n := n / p; od; if n < (p+1)^2 and 1 < n then Add(factors,n); n := 1; fi; if n = 1then CheckAndSortFactors(); return factors; fi;
od;
if effort <= 1then
Add(factors,n); CheckAndSortFactors();
return factors; fi;
# do trial divisions by known primes
atomic readonly Primes2 do for p in Primes2 do while n mod p = 0do Add( factors, p ); n := n / p; od; if n = 1then CheckAndSortFactors(); return factors; fi;
od;
od;
# do trial divisions by known probable primes
tmp := [];
atomic readonly ProbablePrimes2 do for p in ProbablePrimes2 do while n mod p = 0do
AddSet(tmp, p);
Add( factors, p );
n := n / p;
od; if n = 1then break; fi;
od;
od;
if n = 1then CheckAndSortFactors(); return factors; fi;
# handle perfect powers
root := SmallestRootInt( n ); if root < n then
rootfactors := PartialFactorization(root,effort);
k := LogInt(n,root);
rootfactors := Concatenation(List([1..k],i->rootfactors));
factors := SortedList(Concatenation(factors,rootfactors));
CheckAndSortFactors();
return factors; fi;
if effort = 2 or IsProbablyPrimeInt(n) then
Add(factors,n); CheckAndSortFactors(); return factors; fi;
# if effort >= 3, use `FactorsRho' if ValueOption("RhoTrials") <> fail then
tmp := FactorsRho(n,1,16,ValueOption("RhoTrials"):
UseProbabilisticPrimalityTest); else if effort = 3then rhotrials := 256; elif effort = 4then rhotrials := 2048; elif effort >= 5then rhotrials := 8192; fi;
tmp := FactorsRho(n,1,16,rhotrials:UseProbabilisticPrimalityTest); fi;
factors := SortedList(Concatenation(factors,tmp[1],tmp[2]));
CheckAndSortFactors();
return factors; end );
############################################################################# ## #M PartialFactorization( <n> ) . . . . . partial factorization of an integer ##
InstallOtherMethod( PartialFactorization, "for integers", true, [ IsInt ], 0,
n -> PartialFactorization(n,5) );
############################################################################# ## #F Gcdex( <m>, <n> ) . . . . . . . . . . greatest common divisor of integers ##
InstallGlobalFunction(Gcdex,function ( m, n )
local f, g, h, fm, gm, hm, q; if0 <= m then f:=m; fm:=1; else f:=-m; fm:=-1; fi; if0 <= n then g:=n; gm:=0; else g:=-n; gm:=0; fi; while g <> 0do
q := QuoInt( f, g );
h := g; hm := gm;
g := f - q * g; gm := fm - q * gm;
f := h; fm := hm;
od; if n = 0then
return rec( gcd := f, coeff1 := fm, coeff2 := 0,
coeff3 := gm, coeff4 := 1 ); else
return rec( gcd := f, coeff1 := fm, coeff2 := (f - fm * m) / n,
coeff3 := gm, coeff4 := (0 - gm * m) / n ); fi; end);
############################################################################# ## #F IsEvenInt( <n> ) . . . . . . . . . . . . . . . . . . test if <n> is even ##
InstallGlobalFunction( IsEvenInt, n -> n mod 2 = 0 );
############################################################################# ## #F IsOddInt( <n> ) . . . . . . . . . . . . . . . . . . . test if <n> is odd ##
InstallGlobalFunction( IsOddInt, n -> n mod 2 = 1 );
############################################################################# ## #F IsPrimePowerInt( <n> ) . . . . . . . . . . . test for a power of a prime ##
InstallGlobalFunction( IsPrimePowerInt, function(n)
local k, r, s, p, l, q, i;
# check the argument if n > 1then k := 2; s := 1; elif n < -1then k := 3; s := -1; n := -n; else return false; fi;
# exclude small divisors, and thereby large exponents for p in Primes do if p*p > n then return true; fi; # n is prime
r := PVALUATION_INT(n, p); if r > 0then if s = -1 and IsEvenInt(r) then return false; fi;
return n = p^r; fi;
od;
l := LogInt( n, p );
# loop over the possible prime divisors of exponents # use Fermat's little theorem to cast out impossible ones: # for suppose we had r such that n = r^k. Then by Fermat, # n^((q-1)/k) = r^(q-1) is congruent 0 or 1 mod q
i := Position(Primes, k); while k <= l do
q := 2*k+1; while not IsPrimeInt(q) do q := q+2*k; od; if PowerModInt( n, (q-1)/k, q ) <= 1then
r := RootInt( n, k ); if r ^ k = n then
n := r;
l := QuoInt( l, k );
continue; fi; fi; if i <> fail and i < Length(Primes) then
i := i + 1;
k := Primes[i]; else # need more primes...
k := NextPrimeInt( k ); # since we are now beyond the primes in Primes, for which we # checked whether they divide n, we might now just as well # test if k divides n, too
r := PVALUATION_INT(n, k); if r > 0then if s = -1 and IsEvenInt(r) then return false; fi;
return n = k^r; fi; fi;
od;
return IsPrimeInt(n); end);
############################################################################# ## #F LcmInt( <m>, <n> ) . . . . . . . . . . least common multiple of integers ##
InstallGlobalFunction(LcmInt, LCM_INT);
############################################################################# ## #F LogInt( <n>, <base> ) . . . . . . . . . . . . . . logarithm of an integer ##
InstallGlobalFunction(LogInt,function ( n, base )
local log, p;
# check arguments if not IsInt(n) or n <= 0then
Error("<n> must be a positive integer"); fi; if not IsInt(base) or base <= 1then
Error("<base> must be an integer greater than 1"); fi;
# `log(b)' returns $log_b(n)$ and divides `n' by `b^log(b)' ## log := function ( b ) ## local i; ## if b > n then return 0; fi; ## i := log( b^2 ); ## if b > n then return 2 * i; ## else n := QuoInt( n, b ); return 2 * i + 1; fi; ## end; ## ## return log( base ); if n < base then
return 0; elif base = 2then
return Log2Int(n); elif base = 8then
return QuoInt(Log2Int(n), 3); elif base = 16then
return QuoInt(Log2Int(n), 4); elif IsSmallIntRep(n) then
log := 1;
p := base * base; while p <= n do
log := log + 1;
p := p * base;
od;
return log; elif base = 10then
log := QuoInt(Log2Int(n) * 10^6 , 3321929);
return log + LogInt(QuoInt(n, 10^log), 10); else
log := QuoInt(Log2Int(n), Log2Int(base)+1); if log = 0then
log := 1; fi;
return log + LogInt(QuoInt(n, base^log), base); fi; end);
############################################################################# ## #F NextPrimeInt( <n> ) . . . . . . . . . . . . . . . . . . next larger prime ##
InstallGlobalFunction(NextPrimeInt,function ( n ) if -3 = n then n := -2; elif -3 < n and n < 2then n := 2; elif n mod 2 = 0then n := n+1; else n := n+2; fi; while not IsPrimeInt(n) do if n mod 6 = 1then n := n+4; else n := n+2; fi;
od;
return n; end);
############################################################################# ## #F PowerModInt(<r>,<e>,<m>) . . . . . . power of one integer modulo another ##
InstallGlobalFunction(PowerModInt, POWERMODINT);
############################################################################# ## #F PrevPrimeInt( <n> ) . . . . . . . . . . . . . . . previous smaller prime ## ## `PrevPrimeInt' returns the largest prime which is strictly smaller than ## the integer <n>. ##
InstallGlobalFunction(PrevPrimeInt,function ( n ) if3 = n then n := 2; elif -2 < n and n < 3then n := -2; elif n mod 2 = 0then n := n-1; else n := n-2; fi; while not IsPrimeInt(n) do if n mod 6 = 5then n := n-4; else n := n-2; fi;
od;
return n; end);
############################################################################# ## #F PrimePowerInt( <n> ) . . . . . . . . . . . . . . . . prime powers of <n> ##
InstallGlobalFunction(PrimePowersInt,function( n ) if n = 1then
return []; elif n = 0then
Error( "<n> must be non zero" ); elif n < 0then
n := -1 * n; fi;
return Flat(Collected(Factors(Integers,n)));
# get the arguments if Length(arg) = 1then n := arg[1]; k := 2; elif Length(arg) = 2then n := arg[1]; k := arg[2]; else Error("usage: `Root( <n> )' or `Root( <n>, <k> )'"); fi;
# ask the kernel to compute the root; this can only fail for # huge values of k and enormously huge values of n
r := ROOT_INT(n, k); if r <> fail then
return r; fi;
# r is the first approximation, s the second, we need: root <= s < r
r := n; s := 2^( QuoInt( LogInt(n,2), k ) + 1 ) - 1;
# do Newton iterations until the approximations stop decreasing while s < r do
r := s; t := r^(k-1); s := QuoInt( n + (k-1)*r*t, k*t );
od;
# and that's the integer part of the root
return r; end);
############################################################################# ## #F AbsInt( <n> ) . . . . . . . . . . . . . . . absolute value of an integer ## #InstallGlobalFunction( AbsInt, ABS_INT );
InstallGlobalFunction( AbsInt, ABS_RAT ); # support rationals for backwards compatibility
############################################################################# ## #F SmallestRootInt( <n> ) . . . . . . . . . . . smallest root of an integer ##
InstallGlobalFunction(SmallestRootInt,function ( n )
local k, r, s, p, l, q, i;
# check the argument if n > 0then k := 2; s := 1; elif n < 0then k := 3; s := -1; n := -n; else return 0; fi;
# exclude small divisors, and thereby large exponents for p in Primes do if p*p > n then return s * n; fi; if n mod p = 0then break; fi;
od;
l := LogInt( n, p );
# loop over the possible prime divisors of exponents # use Fermat's little theorem to cast out impossible ones: # for suppose we had r such that n = r^k. Then by Fermat, # n^((q-1)/k) = r^(q-1) is congruent 0 or 1 mod q
i := Position(Primes, k); while k <= l do
q := 2*k+1; while not IsPrimeInt(q) do q := q+2*k; od; if PowerModInt( n, (q-1)/k, q ) <= 1then
r := RootInt( n, k ); if r ^ k = n then
n := r;
l := QuoInt( l, k );
continue; fi; fi; if i <> fail and i < Length(Primes) then
i := i + 1;
k := Primes[i]; else
k := NextPrimeInt( k ); fi;
od;
return s * n; end);
############################################################################# ## #M RingByGenerators( <elms> ) . . . . . . . ring generated by some integers ##
InstallMethod( RingByGenerators, "method that catches the cases of `Integers' and subrings of `Integers'",
[ IsCyclotomicCollection ],
SUM_FLAGS, # test this before doing anything else function( elms ) if ForAll( elms, IsInt ) then # check that the number of generators is bigger than one # to avoid infinite recursion if Length( elms ) > 1then
return RingByGenerators( [ Gcd(elms) ] ); elif elms[1] = 1then
return Integers; else
TryNextMethod(); fi; else
TryNextMethod(); fi; end );
############################################################################# ## #M RingWithOneByGenerators( <elms> ) . . . . ring generated by some integers ##
InstallMethod( RingWithOneByGenerators, "method that catches the cases of `Integers'",
[ IsCyclotomicCollection ],
SUM_FLAGS, # test this before doing anything else function( elms ) if ForAll( elms, IsInt ) then
return Integers; else
TryNextMethod(); fi; end );
############################################################################# ## #M DefaultRingByGenerators( <elms> ) default ring generated by some integers ##
InstallMethod( DefaultRingByGenerators, "method that catches the cases of `(Gaussian)Integers' and cycl. fields",
[ IsCyclotomicCollection ],
SUM_FLAGS, # test this before doing anything else function( elms ) if ForAll( elms, IsInt ) then
return Integers; elif ForAll( elms, IsGaussInt ) then
return GaussianIntegers; else
return DefaultField( elms ); fi; end );
############################################################################# ## #M DefaultRingByGenerators( <mats> ) . for a list of n x n integer matrices ##
InstallMethod( DefaultRingByGenerators, "for lists of n x n integer matrices", true,
[ IsCyclotomicCollCollColl and IsFinite ],
function ( mats )
local d; if IsEmpty(mats) or not ForAll(mats,IsRectangularTable and IsMatrix) then
TryNextMethod(); fi;
d := Length( mats[1] ); if d=0then
TryNextMethod(); fi; if not ForAll( mats, m -> Length(m)=d and Length(m[1])=d ) then
TryNextMethod(); fi; if not ForAll( mats, m -> ForAll( m, r -> ForAll(r,IsInt))) then
TryNextMethod(); fi;
return FullMatrixAlgebra(Integers,d); end );
############################################################################# ## #M Enumerator( Integers ) ## ## $a_n = \frac{n}{2}$ if $n$ is even, and ## $a_n = \frac{1-n}{2}$ otherwise. ##
InstallMethod( Enumerator, "for integers",
[ IsIntegers ],
Integers -> EnumeratorByFunctions( Integers,
rec( ElementNumber := function( e, n ) if n mod 2 = 0then
return n / 2; else
return ( 1 - n ) / 2; fi; end,
NumberElement := function( e, x )
local pos; if not IsInt( x ) then
return fail; elif0 < x then
pos:= 2 * x; else
pos:= -2 * x + 1; fi;
return pos; end ) ) );
############################################################################# ## #M EuclideanDegree( Integers, <n> ) . . . . . . . . . . . . . absolute value ##
InstallMethod( EuclideanDegree, "for integers",
true,
[ IsIntegers, IsInt ], 0, function ( Integers, n ) if n < 0then
return -n; else
return n; fi; end );
############################################################################# ## #M EuclideanQuotient( Integers, <n>, <m> ) . . . . . . Euclidean quotient ##
InstallMethod( EuclideanQuotient, "for integers",
true,
[ IsIntegers, IsInt, IsInt ], 0, function ( Integers, n, m )
return QuoInt( n, m ); end );
############################################################################# ## #M EuclideanRemainder( Integers, <n>, <m> ) . . . . . . Euclidean remainder ##
InstallMethod( EuclideanRemainder, "for integers",
true,
[ IsIntegers, IsInt, IsInt ], 0, function ( Integers, n, m )
return RemInt( n, m ); end );
############################################################################# ## #M Factors( Integers, <n> ) . . . . . . . . . . factorization of an integer ##
InstallMethod( Factors, "for integers",
true,
[ IsIntegers, IsInt ], 0, function ( Integers, n )
return FactorsInt( n ); end );
############################################################################# ## #M IsIrreducibleRingElement( Integers, <n> ) ##
InstallMethod( IsIrreducibleRingElement, "for integers",
true,
[ IsIntegers, IsInt ], 0, function ( Integers, n )
return IsPrimeInt( n ); end );
############################################################################# ## #M IsPrime( Integers, <n> ) . . . . . . test whether an integer is a prime ##
InstallMethod( IsPrime, "for integers",
true,
[ IsIntegers, IsInt ], 0, function ( Integers, n )
return IsPrimeInt( n ); end );
############################################################################# ## #M Iterator( Integers ) ## ## uses the succession $0, 1, -1, 2, -2, 3, -3, \ldots$, that is, ## $a_n = \frac{n}{2}$ if $n$ is even, and $a_n = \frac{1-n}{2}$ ## otherwise. ##
InstallMethod( Iterator, "for `Integers'",
[ IsIntegers ],
Integers -> IteratorByFunctions( rec(
NextIterator := function( iter )
iter!.counter:= iter!.counter + 1; if iter!.counter mod 2 = 0then
return iter!.counter / 2; else
return ( 1 - iter!.counter ) / 2; fi; end,
IsDoneIterator := ReturnFalse,
ShallowCopy := iter -> rec( counter:= iter!.counter ),
PrintObj := function(iter)
local msg;
msg := "<iterator of Integers at "; if iter!.counter mod 2 = 0then
Append(msg, String(iter!.counter / 2)); else
Append(msg, String((1 - iter!.counter) / 2)); fi;
Append(msg,">");
Print(msg); end,
counter := 0 ) ) );
counter := 0 ) ) ); # 0, since we first increment then return
############################################################################# ## #M LcmOp( Integers, <n>, <m> ) . . . . . . least common multiple of integers ##
InstallMethod( LcmOp, "for integers",
true,
[ IsIntegers, IsInt, IsInt ], 0, function ( Integers, n, m )
return LcmInt( n, m ); end );
############################################################################# ## #M PowerMod( Integers, <r>, <e>, <m> ) . . . power of an integer mod another ##
InstallMethod( PowerMod, "for integers",
true,
[ IsIntegers, IsInt, IsInt, IsInt ], 0, function ( Integers, r, e, m )
return PowerModInt( r, e, m ); end );
############################################################################# ## #M Quotient( <Integers>, <n>, <m> ) . . . . . . . quotient of two integers ##
InstallMethod( Quotient, "for integers",
true,
[ IsIntegers, IsInt, IsInt ], 0, function ( Integers, n, m )
local q; if m = 0then
return fail; fi;
q := QuoInt( n, m ); if n <> q * m then
q := fail; fi;
return q; end );
############################################################################# ## #M QuotientMod( Integers , <r>, <s>, <m> ) . . . . . . . quotient modulo <m> ##
InstallMethod( QuotientMod, "for integers",
true,
[ IsIntegers, IsInt, IsInt, IsInt ], 0, function ( Integers, r, s, m )
local g; if m = 1then
return 0; else
r := r mod m; if r = 0then return 0; fi; # as required by QuotientMod documentation
s := s mod m; if s = 0then return fail; fi;
g := GcdInt( r, s );
r := r / g;
s := s / g;
g := GcdInt( g, m );
m := m / g; if GcdInt( s, m ) <> 1then
return fail; fi;
return r * INVMODINT(s, m) mod m; fi; end );
############################################################################# ## #M QuotientRemainder( Integers, <n>, <m> ) . . . . . . . . . . . quo and rem ##
InstallMethod( QuotientRemainder, "for integers",
true,
[ IsIntegers, IsInt, IsInt ], 0, function ( Integers, n, m )
local q;
q := QuoInt(n,m); #T kernel function should compute remainder at same time
return [ q, n - q * m ]; end );
############################################################################# ## #M Random( Integers ) . . . . . . . . . . . . . . . . . . . random integer ## ## returns pseudo random integers between $-10$ and $10$ ## distributed according to a binomial distribution. ## ## \begintt ## gap> Random( Integers ); ## 1 ## gap> Random( Integers ); ## -4 ## \endtt ## ## To generate uniformly distributed integers from a range, use the ## construct `Random( [ <low> .. <high> ] )'. ##
BindGlobal( "NrBitsInt", function ( n )
local nr, nr64;
nr64:=[0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5, 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6];
nr := 0; while0 < n do
nr := nr + nr64[ n mod 64 + 1 ];
n := QuoInt( n, 64 );
od;
return nr; end );
InstallMethodWithRandomSource(Random, "for a random source and `Integers'", true,
[IsRandomSource, IsIntegers], 0, function( rg, Integers )
return NrBitsInt( Random( rg, 0, 2^20-1 ) ) - 10; end );
############################################################################# ## #M StandardAssociate( Integers, <n> ) . . . . . . . . . . . absolute value ##
InstallMethod( StandardAssociate, "for integers",
true,
[ IsIntegers, IsInt ], 0, function ( Integers, n ) if n < 0then
return -n; else
return n; fi; end );
############################################################################# ## #M StandardAssociateUnit( Integers, <n> ) ##
InstallMethod( StandardAssociateUnit, "for integers",
true,
[ IsIntegers, IsInt ], 0, function ( Integers, n ) if n < 0then
return -1; else
return 1; fi; end );
############################################################################# ## #M Valuation( <n>, <m> ) ##
InstallOtherMethod( Valuation, "for two integers",
IsIdenticalObj,
[ IsInt, IsInt ], 0, function( n, m ) if n = 0then
return infinity; fi;
return PVALUATION_INT( n, m ); end );
############################################################################# ## #M \in( <n>, <Integers> ) . . . . . . . . . . membership test for integers ##
InstallMethod( \in, "for integers",
IsElmsColls,
[ IsCyclotomic, IsIntegers ], 0, function( n, Integers )
return IsInt( n ); end );
############################################################################# ## #M \in( <n>, <PositiveIntegers> ) ##
InstallMethod( \in, "for positive integers",
IsElmsColls,
[ IsCyclotomic, IsPositiveIntegers ], 0, function( n, PositiveIntegers )
return IsPosInt( n ); end );
############################################################################# ## #M \in( <n>, <NonnegativeIntegers> ) ##
InstallMethod( \in, "for nonnegative integers",
IsElmsColls,
[ IsCyclotomic, IsNonnegativeIntegers ], 0, function( n, NonnegativeIntegers )
return IsPosInt( n ) or IsZeroCyc( n ); end );
############################################################################# ## #F PrintFactorsInt( <n> ) . . . . . . . . print factorization of an integer ##
InstallGlobalFunction(PrintFactorsInt,function ( n )
Print( StringPP( n ) ); end);
############################################################################# ## #M Iterator( <posint> ) . . . . . . . . . . . . .give more informative error ## ## This method is mainly there to trap the "natural" error ## for i in 3 do ... od; ##
InstallOtherMethod(Iterator, "more helpful error for integers", true,
[IsPosInt], 0, function(n)
Error("You cannot loop over the integer ",n, " did you mean the range [1..",n,"]"); end);
## The behaviour of View(String) for large integers can be configured via a ## user preference.
DeclareUserPreference( rec(
name:= "MaxBitsIntView",
description:= [ "Maximal bit length of integers to <C>View</C> unabbreviated. \
Default is about <M>30</M> lines of a <M>80</M> character wide terminal. \
Set this to <C>0</C> to avoid abbreviated ints."
],
default:= 8000,
check:= val -> IsInt( val ) and 0 <= val,
) ); ## give only a short info if |n| is larger than 2^GAPInfo.MaxBitsIntView
InstallMethod(ViewString, "for integer", [IsInt], function(n)
local mb, l, start, trail;
mb := UserPreference("MaxBitsIntView"); if not IsSmallIntRep(n) and mb <> fail and
mb > 64 and Log2Int(n) > mb then if n < 0then
l := LogInt(-n, 10);
trail := String(-n mod 1000); else
l := LogInt(n, 10);
trail := String(n mod 1000); fi;
start := String(QuoInt(n, 10^(l-2))); while Length(trail) < 3do
trail := Concatenation("0", trail);
od;
return Concatenation("<integer ",start,"...",trail," (",
String(l+1)," digits)>"); else
return String(n); fi; end);
Messung V0.5 in Prozent
¤ Diese beiden folgenden Angebotsgruppen bietet das Unternehmen0.70Angebot
(Wie Sie bei der Firma Beratungs- und Dienstleistungen beauftragen können 2026-10-13)
¤
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.
