Quelle numtheor.gi Sprache: Shell

#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Martin Schönert.
##
## 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
##
## This file contains methods mainly for integer primes.
##

#############################################################################
##
#F PrimeResidues( <m> ) . . . . . . . integers relative prime to an integer
##
BindGlobal( "PrimeResiduesCache",
 List( [[],[0],[1],[1,2],[1,3],[1,2,3,4],[1,5],[1,2,3,4,5,6]], Immutable ));

if IsHPCGAP then
 MakeImmutable(PrimeResiduesCache);
fi;

InstallGlobalFunction( PrimeResidues, function ( m )
 local residues, p, i;

 # make <m> it nonnegative, handle trivial cases
 if m < 0 then m := -m; fi;
 if m < Length(PrimeResiduesCache) then
 return ShallowCopy(PrimeResiduesCache[m+1]);
 fi;

 # remove the multiples of all prime divisors
 residues := [1..m-1];
 for p in PrimeDivisors(m) do
 for i in [1..m/p-1] do
 residues[p*i] := 1;
 od;
 od;

 # return the set of residues
 return Set( residues );
end );

#############################################################################
##
#M Phi( <m> ) . . . . . . . . . . . . . . . . . . Euler's totient function
##
InstallMethod( Phi,
 "value of Euler's totient function of an integer",
 true, [ IsInt ], 0,

 function ( m )

 local phi, p;

 # make <m> it nonnegative, handle trivial cases
 if m < 0 then m := -m; fi;
 if m < Length(PrimeResiduesCache) then
 return Length(PrimeResiduesCache[m+1]);
 fi;

 # compute $phi$
 phi := m;
 for p in PrimeDivisors(m) do
 phi := phi / p * (p-1);
 od;

 # return the result
 return phi;
 end );

#############################################################################
##
#M Lambda( <m> ) . . . . . . . . . . . . . . . . . . . . Carmichael function
##
InstallMethod( Lambda,
 "exponent of the group of prime residues modulo an integer",
 true, [ IsInt ], 0,

 function ( m )

 local lambda, p, q, k;

 # make <m> it nonnegative, handle trivial cases
 if m < 0 then m := -m; fi;
 if m < Length(PrimeResiduesCache) then
 return Length(PrimeResiduesCache[m+1]);
 fi;

 # loop over all prime factors $p$ of $m$
 lambda := 1;
 for p in PrimeDivisors(m) do

 # compute $p^e$ and $k = (p-1) p^(e-1)$
 q := p; k := p-1;
 while m mod (q * p) = 0 do q := q * p; k := k * p; od;

 # multiples of 8 are special
 if q mod 8 = 0 then k := k / 2; fi;

 # combine with the value known so far
 lambda := LcmInt( lambda, k );
 od;

 return lambda;
 end );

#############################################################################
##
#F OrderMod( <n>, <m>[, <bound>] ) . . . multiplicative order of an integer
##
#N 23-Apr-91 martin improve 'OrderMod' similar to 'IsPrimitiveRootMod'
##
InstallGlobalFunction( OrderMod, function ( n, m, bound... )
 local x, o, d;

 # check the arguments and reduce $n$ into the range $0..m-1$
 if m <= 0 then Error("<m> must be positive"); fi;
 if n < 0 then n := n mod m + m; fi;
 if m <= n then n := n mod m; fi;

 # return 0 if $m$ is not coprime to $n$
 if GcdInt(m,n) <> 1 then
 o := 0;

 # compute the order simply by iterated multiplying, $x= n^o$ mod $m$
 elif m < 100 then
 x := n; o := 1;
 while x > 1 do
 x := x * n mod m; o := o + 1;
 od;

 # otherwise try the divisors of $\lambda(m)$ and their divisors, etc.
 else
 if Length( bound ) = 1 then
 # We know a multiple of the desired order.
 o := bound[1];
 else
 # The default a priori known multiple is 'Lambda( m )'.
 o := Lambda( m );
 fi;
 for d in PrimeDivisors( o ) do
 while o mod d = 0 and PowerModInt(n,o/d,m) = 1 do
 o := o / d;
 od;
 od;

 fi;

 return o;
end );

#############################################################################
##
#F IsPrimitiveRootMod( <r>, <m> ) . . . . . . . . test for a primitive root
##
InstallGlobalFunction( IsPrimitiveRootMod, function ( r, m )
 local p, facs, pows, i, pow;

 # check the arguments and reduce $r$ into the range $0..m-1$
 if m <= 0 then Error("<m> must be positive"); fi;
 if r < 0 then r := r mod m + m; fi;
 if m <= r then r := r mod m; fi;

 # handle trivial cases
 if m = 2 then return r = 1; fi;
 if m = 4 then return r = 3; fi;
 if m mod 4 = 0 then return false; fi;

 # handle even numbers by testing modulo the odd part
 if m mod 2 = 0 then
 if r mod 2 = 0 then return false; fi;
 m := m / 2;
 fi;

 # check that $m$ is a prime power, otherwise no primitive root exists
 p := SmallestRootInt( m );
 if not IsPrimeInt( p ) then
 return false;
 fi;

 # check that $r^((p-1)/2) \<> 1$ mod $p$ using the Jacobi symbol
 if Jacobi( r, p ) <> -1 then
 return false;
 fi;

 # compute $pows_i := r^{{p-1}/\prod_{k=2}^{i}{facs_k}}$ ($facs_1 = 2$)
 facs := PrimeDivisors( p-1 );
 pows := [];
 pows[Length(facs)] := PowerModInt( r, 2*(p-1)/Product(facs), p );
 for i in Reversed( [2..Length(facs)-1] ) do
 pows[i] := PowerModInt( pows[i+1], facs[i+1], p );
 od;

 # check $1 \<> r^{{p-1}/{facs_i}} = pows_i^{\prod_{k=2}^{i-1}{facs_k}}$
 pow := 1;
 for i in [2..Length(facs)] do
 if PowerModInt( pows[i], pow, p ) = 1 then
 return false;
 fi;
 pow := pow * facs[i];
 od;

 # if $m$ is a prime we are done
 if p = m then
 return true;
 fi;

 # for prime powers $n$ we have to test that $r$ is not a $p$-th root
 return PowerModInt( r, p-1, p^2 ) <> 1;
end );

#############################################################################
##
#F PrimitiveRootMod( <m> ) . . . . . . . . primitive root modulo an integer
##
InstallGlobalFunction( PrimitiveRootMod, function ( arg )
 local root, m, p, start, mm;

 # get and check the arguments
 if Length(arg) = 1 then
 m := arg[1]; start := 1;
 if m <= 0 then Error("<m> must be positive"); fi;
 elif Length(arg) = 2 then
 m := arg[1]; start := arg[2];
 if m <= 0 then Error("<m> must be positive"); fi;
 else
 Error("usage: PrimitiveRootMod( <m>[, <start>] )");
 fi;

 # handle the trivial cases
 if m = 2 and start = 1 then return 1; fi;
 if m = 2 then return fail; fi;
 if m = 4 and start <= 3 then return 3; fi;
 if m mod 4 = 0 then return fail; fi;

 # handle even numbers
 if m mod 2 = 0 then
 mm := 2;
 m := m / 2;
 else
 mm := 1;
 fi;

 # check that $m$ is a prime power otherwise no primitive root exists
 p := SmallestRootInt( m );
 if not IsPrimeInt( p ) then
 return fail;
 fi;

 # run through all candidates for a primitive root
 root := start+1;
 while root <= mm*m-1 do
 if (mm = 1 or root mod 2 = 1)
 and IsPrimitiveRootMod( root, p )
 and (p = m or PowerModInt( root, p-1, p^2 ) <> 1)
 then
 return root;
 fi;
 root := root + 1;
 od;

 # no primitive root found
 return fail;
end );

#############################################################################
##
#F GeneratorsPrimeResidues( <n> ) . . . . . . generators of the Galois group
##
InstallGlobalFunction( GeneratorsPrimeResidues, function( n )
 local factors, # collected list of prime factors of `n'
 primes, # list of prime divisors of `n'
 exponents, # exponents of the entries in `primes'
 generators, # generator(s) of the prime part `ppart' of residues
 ppart, # one prime part of `n'
 rest, # `n / ppart'
 i, # loop over the positions in `factors'
 gcd; # one coefficient in the output of `Gcdex'

 if n = 1 then
 return rec(
 primes := [],
 exponents := [],
 generators := []
 );
 fi;

 factors:= Collected( Factors(Integers, n ) );

 primes := [];
 exponents := [];
 generators := [];

 # For each prime part `ppart',
 # the generator must be congruent to a primitive root modulo `ppart',
 # and congruent 1 modulo the rest `N/ppart'.

 for i in [ 1 .. Length( factors ) ] do

 primes[i] := factors[i][1];
 exponents[i] := factors[i][2];
 ppart := primes[i] ^ exponents[i];
 rest := n / ppart;

 if primes[i] = 2 then

 gcd:= Gcdex( ppart, rest ).coeff2 * rest;
 if ppart mod 8 = 0 then
 # Choose the generators `**' and `*5'.
 generators[i]:= [ ( -2 * gcd + 1 ) mod n, # generator `**'
 ( 4 * gcd + 1 ) mod n ]; # generator `*5'
 else
 # Choose the generator `**'.
 generators[i]:= ( -2 * gcd + 1 ) mod n;
 fi;

 else
 generators[i] := ( ( PrimitiveRootMod( ppart ) - 1 )
 * Gcdex( ppart, rest ).coeff2 * rest + 1 ) mod n;
 fi;

 od;

 return rec(
 primes := primes,
 exponents := exponents,
 generators := generators
 );
end );

#############################################################################
##
#F Jacobi( <n>, <m> ) . . . . . . . . . . . . . . . . . . . . Jacobi symbol
##
##
InstallGlobalFunction( Jacobi, JACOBI_INT );

#############################################################################
##
#F Legendre( <n>, <m> ) . . . . . . . . . . . . . . . . . . Legendre symbol
##
InstallGlobalFunction( Legendre, function ( n, m )
 local p, q, o;

 # check the arguments and reduce $n$ into the range $0..m-1$
 if m <= 0 then Error("<m> must be positive"); fi;
 if n < 0 then n := n mod m + m; fi;
 if m <= n then n := n mod m; fi;

 # handle small values
 if n in [0,1,4,9,16] then return 1; fi;
 if m < 6 then return -1; fi;

 # check that $n$ is a quadratic residue modulo every prime power $q$
 for p in PrimeDivisors( m ) do

 # find prime power $q$ and reduce $n$
 q := p; while m mod (q * p) = 0 do q := q * p; od;
 o := n mod q;

 # the largest power of $p$ that divides $o$ must be even
 while o > 0 and o mod p = 0 do
 if o mod p^2 <> 0 then return -1; fi;
 o := o / p^2;
 od;

 # check that $m$ is a quadratic residue modulo $p$ with Jacobi
 if p = 2 then
 if o > 0 and o mod Minimum(8,q) <> 1 then return -1; fi;
 else
 if Jacobi( o, p ) = -1 then return -1; fi;
 fi;

 od;

 # else $n$ is a quadratic residue modulo $m$
 return 1;
end );

#############################################################################
##
#F RootMod( <n>, <m> ) . . . . . . . . . . . . . . . root modulo an integer
#F RootMod( <n>, <k>, <m> ) . . . . . . . . . . . . root modulo an integer
##
## Carried over from GAP3: This code requires that k is a prime.
BindGlobal( "RootModPrime", function ( n, k, p )
 local r, # <k>th root of <n> mod , result
 kk, # largest power of <k> dividing -1
 s, # cofactor for <r>
 ss, # power of <s>
 t, # <kk>th root of unity mod 
 tt, # power of <t>
 i; # loop variable

 # reduce $n$ into the range $0..p-1$
 Info( InfoNumtheor, 1, "RootModPrime(", n, ",", k, ",", p, ")" );
 n := n mod p;

 # If n is non-zero, then the code below requires that p is a Fermat
 # pseudoprime with respect to n, i.e. that $n^(p-1) mod p = 1$ holds.
 # This is of course automatically true if $p$ is a prime, but for efficiency
 # reasons we actually use this with pseudo primes.
 if n <> 0 and PowerModInt( n, p-1, p ) <> 1 then
 Error( " is not a Fermat pseudoprime with respect to <n>. Please report this error to the GAP team" );
 fi;

 # handle $p = 2$
 if p = 2 then
 r := n;

 # handle $n = 0$ (only possibility that <n> and are not rel. prime)
 elif n = 0 then
 r := 0;

 # it's easy if $k$ is invertible mod $p-1 = \phi(p)$
 elif GcdInt( p-1, k ) = 1 then
 Info( InfoNumtheor, 2, " <r> = <n>^", 1/k mod (p-1), " mod " );
 r := PowerModInt( n, 1/k mod (p-1), p );

 # check that $n$ has a $k$th root (Euler's criterium)
 elif PowerModInt( n, (p-1)/k, p ) = 1 then

 # $p-1 = x kk$, $x$ mod $k <> 0$
 kk := 1; while (p-1)/kk mod k = 0 do kk := kk*k; od;
 Info( InfoNumtheor, 2,
 " ", p, "-1 = <x> * ", kk, ", <x> mod ", k, " <> 1" );

 # find $r$ up to a $kk$-th root of 1, i.e., $n s = r^k, s^{kk/k} = 1$
 r := PowerModInt( n, 1/k mod ((p-1)/kk), p );
 s := PowerModInt( r, k, p ) / n mod p;
 Info( InfoNumtheor, 2,
 " <n>*", s, "=", r, "^", k, ", ", s, "^", kk/k, "=1" );

 # find a generator $t$ of the subgroup of $kk$-th roots of 1,
 # i.e., $t^{kk/k} <> 1, t^{kk} = 1$, therefore $s = (t^l)^k$
 i:=2; t:=PowerModInt(i,(p-1)/kk,p); tt:=PowerModInt(t,kk/k,p);
 while tt=1 do
 i:=i+1; t:=PowerModInt(i,(p-1)/kk,p); tt:=PowerModInt(t,kk/k,p);
 od;
 Info( InfoNumtheor, 2,
 " ", t, "^", kk/k, " <> 1, ", t, "^", kk, " = 1" );

 # $n s = r^k, s^{kk/k} = 1, t^{kk/k} <> 1, t^{kk} = 1$
 while kk <> k do
 Info( InfoNumtheor, 2,
 " <n>*", s, "=", r, "^", k, ", ", s, "^", kk/k, "=1" );
 kk := kk/k;
 i := t;
 t := PowerModInt( t, k, p );
 ss := PowerModInt( s, kk/k, p );
 while ss <> 1 do
 r := r * i mod p;
 s := s * t mod p;
 ss := ss * tt mod p;
 od;
 od;
 Info( InfoNumtheor, 2, " <n>*1=", r, "^", k, ", 1^1=1" );

 # otherwise $n$ has no root
 else
 r := fail;

 fi;

 # return the root $r$
 Info( InfoNumtheor, 1, "RootModPrime returns ", r );
 return r;
end );

DeclareGlobalName( "RootModPrimePower" );
BindGlobal( "RootModPrimePower", function ( n, k, p, l )
 local r, # <k>th root of <n> mod ^<l>, result
 s, # <k>th root of <n> mod smaller power
 t; # temporary variable

 # delegate prime case
 Info( InfoNumtheor, 1,
 "RootModPrimePower(", n, ",", k, ",", p, "^", l, ")" );
 if l = 1 then
 r := RootModPrime( n, k, p );

 # special case
 elif n mod p^l = 0 then
 r := 0;

 # if $n$ is a multiple of $p^k$ return $p (\sqrt[k]{n/p^k} mod p^l/p^k)$
 elif n mod p^k = 0 then
 s := RootModPrimePower( n/p^k, k, p, l-k );
 if s <> fail then
 r := s * p;
 else
 r := fail;
 fi;

 # if $n$ is a multiple of $p$ but not of $p^k$ then no root exists
 elif n mod p = 0 then
 r := fail;

 # handle the case that the root may not lift
 elif k = p then

 Info( InfoNumtheor, 3, "k=p case" );

 # compute the root mod $p^{l/2}$, or $p^{l/2+1}$ if 32 divides $p^l$
 if 2 < p or l < 5 then
 s := RootModPrimePower( n, k, p, QuoInt(l+1,2) );
 else
 s := RootModPrimePower( n, k, p, QuoInt(l+3,2) );
 fi;

 if s=fail then
 r:=fail;
 else
 # lift the root to $p^l$, use higher precision
 Info( InfoNumtheor, 2, " lift root with Newton / Hensel" );
 t := PowerModInt( s, k-1, p^(l+1) );
 r := (s + (n - t * s) / (k * t)) mod p^l;
 if PowerModInt(r,k,p^l) <> n mod p^l then
 r := fail;
 fi;
 fi;

 # otherwise lift the root with Newton / Hensel
 else

 # compute the root mod $p^{l/2}$, or $p^{l/2+1}$ if 32 divides $p^l$
 if 2 < p or l < 5 then
 s := RootModPrimePower( n, k, p, QuoInt(l+1,2) );
 else
 s := RootModPrimePower( n, k, p, QuoInt(l+3,2) );
 fi;
 Info( InfoNumtheor, 3, "lift case s=",s );

 if s=fail then
 r:=fail;
 else
 # lift the root to $p^l$
 Info( InfoNumtheor, 2, " lift root with Newton / Hensel" );
 t := PowerModInt( s, k-1, p^l );
 r := (s + (n - t * s) / (k * t)) mod p^l;
 fi;

 fi;

 # return the root $r$
 Info( InfoNumtheor, 1, "RootModPrimePower returns ", r );
 return r;
end );

InstallGlobalFunction( RootMod, function ( arg )
 local n, # <n>, first argument
 k, # <k>, optional second argument
 m, # <m>, third argument
 p, # prime divisor of <m>
 q, # power of 
 l, # <q> = ^<l>
 qq, # product of prime powers dividing <m>
 ii, # inverse of <qq> mod <q>
 r, # <k>th root of <n> mod <qq>
 s, # <k>th root of <n> mod <q>
 f, # factors
 i; # loop

 # get the arguments
 if Length(arg) = 2 then n := arg[1]; k := 2; m := arg[2];
 elif Length(arg) = 3 then n := arg[1]; k := arg[2]; m := arg[3];
 else Error("usage: RootMod( <n>, <m> ) or RootMod( <n>, <k>, <m> )");
 fi;
 Info( InfoNumtheor, 1, "RootMod(", n, ",", k, ",", m, ")" );

 # check the arguments and reduce $n$ into the range $0..m-1$
 if m <= 0 then Error("<m> must be positive"); fi;
 n := n mod m;

 if not IsPrime(k) then
 # try over factors of k
 f:=Factors(k);
 l:=n;
 for i in f do
 l:=RootMod(l,i,m);
 if l=fail then
 Info( InfoNumtheor, 2, "must try multiple roots");
 # it failed. This might have been because of taking the wrong root
 # do again with all roots
 l:=RootsMod(n,k,m);
 if Length(l)=0 then
 return fail;
 else
 return l[1];
 fi;

 fi;

 od;
 return l;
 fi;

 # combine the root modulo every prime power $p^l$
 r := 0; qq := 1;
 for p in PrimeDivisors( m : UseProbabilisticPrimalityTest ) do

 # find prime power $q = p^l$
 q := p; l := 1;
 while m mod (q * p) = 0 do q := q * p; l := l + 1; od;

 # compute the root mod $p^l$
 s := RootModPrimePower( n, k, p, l );
 if s = fail then
 Info( InfoNumtheor, 1, "RootMod returns 'fail'" );
 return fail;
 fi;

 # combine $r$ (the root mod $qq$) with $s$ (the root mod $p^l$)
 ii := 1/qq mod q;
 r := r + qq * ((s - r)*ii mod q);
 qq := qq * q;

 od;

 # return the root $rr$
 Info( InfoNumtheor, 1, "RootMod returns ", r );
 return r;
end );

#############################################################################
##
#F RootsMod( <n>, <k>, <m> ) . . . . . . . . . . . . roots modulo an integer
##
BindGlobal( "RootsModPrime", function ( n, k, p )
 local rr, # <k>th roots of <n> mod , result
 r, # one particular <k>th root of <n> mod 
 kk, # largest power of <k> dividing -1
 s, # cofactor for <r>
 ss, # power of <s>
 t, # <kk>th root of unity mod 
 tt, # power of <t>
 i; # loop variable

 # reduce $n$ into the range $0..p-1$
 Info( InfoNumtheor, 1, "RootsModPrime(", n, ",", k, ",", p, ")" );
 n := n mod p;

 # If n is non-zero, then the code below requires that p is a Fermat
 # pseudoprime with respect to n, i.e. that $n^(p-1) mod p = 1$ holds.
 # This is of course automatically true if $p$ is a prime, but for efficiency
 # reasons we actually use this with pseudo primes.
 if n <> 0 and PowerModInt( n, p-1, p ) <> 1 then
 Error( " is not a Fermat pseudoprime with respect to <n>. Please report this error to the GAP team" );
 fi;

 # handle $p = 2$
 if p = 2 then
 rr := [ n ];

 # handle $n = 0$ (only possibility that <n> and are not rel. prime)
 elif n = 0 then
 rr := [ 0 ];

 # it's easy if $k$ is invertible mod $p-1 = \phi(p)$
 elif GcdInt( p-1, k ) = 1 then
 Info( InfoNumtheor, 2, " <r> = <n>^", 1/k mod (p-1), " mod " );
 rr := [ PowerModInt( n, 1/k mod (p-1), p ) ];

 # check that $n$ has a $k$th root (Euler's criterium)
 elif PowerModInt( n, (p-1)/k, p ) = 1 then

 # $p-1 = x kk$, $x$ mod $k <> 0$
 kk := 1; while (p-1)/kk mod k = 0 do kk := kk*k; od;
 Info( InfoNumtheor, 2,
 " ", p, "-1 = <x> * ", kk, ", <x> mod ", k, " <> 1" );

 # find $r$ up to a $kk$-th root of 1, i.e., $n s = r^k, s^{kk/k} = 1$
 r := PowerModInt( n, 1/k mod ((p-1)/kk), p );
 s := PowerModInt( r, k, p ) / n mod p;
 Info( InfoNumtheor, 2,
 " <n>*", s, "=", r, "^", k, ", ", s, "^", kk/k, "=1" );

 # find a generator $t$ of the subgroup of $kk$-th roots of 1,
 # i.e., $t^{kk/k} <> 1, t^{kk} = 1$, therefore $s = (t^l)^k$
 i:=2; t:=PowerModInt(i,(p-1)/kk,p); tt:=PowerModInt(t,kk/k,p);
 while tt=1 do
 i:=i+1; t:=PowerModInt(i,(p-1)/kk,p); tt:=PowerModInt(t,kk/k,p);
 od;
 Info( InfoNumtheor, 2,
 " ", t, "^", kk/k, " <> 1, ", t, "^", kk, " = 1" );

 # $n s = r^k, s^{kk/k} = 1, t^{kk/k} <> 1, t^{kk} = 1$
 while kk <> k do
 Info( InfoNumtheor, 2,
 " <n>*", s, "=", r, "^", k, ", ", s, "^", kk/k, "=1" );
 kk := kk/k;
 i := t;
 t := PowerModInt( t, k, p );
 ss := PowerModInt( s, kk/k, p );
 while ss <> 1 do
 r := r * i mod p;
 s := s * t mod p;
 ss := ss * tt mod p;
 od;
 od;
 Info( InfoNumtheor, 2, " <n>*1=", r, "^", k, ", 1^1=1" );

 # combine $r$ (a particular root) with the powers of $t$
 rr := [ r ];
 for i in [2..k] do
 r := r * t mod p;
 AddSet( rr, r );
 od;

 # otherwise $n$ has no root
 else
 rr := [];

 fi;

 # return the roots $rr$
 Info( InfoNumtheor, 1, "RootsModPrime returns ", rr );
 return rr;
end );

DeclareGlobalName( "RootsModPrimePower" );
BindGlobal( "RootsModPrimePower", function ( n, k, p, l )
 local rr, # <k>th roots of <n> mod ^<l>, result
 r, # one element of <rr>
 ss, # <k>th roots of <n> mod smaller power
 s, # one element of <ss>
 t; # temporary variable

 # delegate prime case
 Info( InfoNumtheor, 1,
 "RootsModPrimePower(", n, ",", k, ",", p, "^", l, ")" );
 if l = 1 then
 rr := RootsModPrime( n, k, p );

 # special case
 elif n mod p^l = 0 then
 t := QuoInt( l-1, k ) + 1;
 rr := [ 0 .. p^(l-t)-1 ] * p^t;

 # if $n$ is a multiple of $p^k$ return $p (\sqrt[k]{n/p^k} mod p^l/p^k)$
 elif n mod p^k = 0 then
 ss := RootsModPrimePower( n/p^k, k, p, l-k );
 rr := [];
 for s in ss do
 for t in [ 0 .. p^(k-1)-1 ] do
 AddSet( rr, s * p + t * p^(l-k+1) );
 od;
 od;

 # if $n$ is a multiple of $p$ but not of $p^k$ then no root exists
 elif n mod p = 0 then
 rr := [];

 # handle the case that the roots split
 elif k = p then

 Info( InfoNumtheor, 3, "k=p case" );

 # compute the root mod $p^{l/2}$, or $p^{l/2+1}$ if 32 divides $p^l$
 if 2 < p or l < 5 then
 ss := RootsModPrimePower( n, k, p, QuoInt(l+1,2) );
 else
 ss := RootsModPrimePower( n, k, p, QuoInt(l+3,2) );
 fi;

 # lift the roots to $p^l$, use higher precision
 rr := [];
 for s in ss do
 Info( InfoNumtheor, 2, " lift root with Newton / Hensel" );
 t := PowerModInt( s, k-1, p^(l+1) );
 r := (s + (n - t * s) / (k * t)) mod p^l;
 if PowerModInt(r,k,p^l) = n mod p^l then
 for t in [0..k-1]*p^(l-1)+1 do
 AddSet( rr, r * t mod p^l );
 od;
 fi;
 od;

 # otherwise lift the roots with Newton / Hensel
 else

 # compute the root mod $p^{l/2}$, or $p^{l/2+1}$ if 32 divides $p^l$
 if 2 < p or l < 5 then
 ss := RootsModPrimePower( n, k, p, QuoInt(l+1,2) );
 else
 ss := RootsModPrimePower( n, k, p, QuoInt(l+3,2) );
 fi;

 # lift the roots to $p^l$
 rr := [];
 for s in ss do
 Info( InfoNumtheor, 2, " lift root with Newton / Hensel" );
 t := PowerModInt( s, k-1, p^l );
 r := (s + (n - t * s) / (k * t)) mod p^l;
 AddSet( rr, r );
 od;

 fi;

 # return the roots $rr$
 Info( InfoNumtheor, 1, "RootsModPrimePower returns ", rr );
 return rr;
end );

InstallGlobalFunction( RootsMod, function ( arg )
 local n, # <n>, first argument
 k, # <k>, optional second argument
 m, # <m>, third argument
 p, # prime divisor of <m>
 q, # power of 
 l, # <q> = ^<l>
 f, # factors
 qq, # product of prime powers dividing <m>
 ii, # inverse of <qq> mod <q>
 rr, # <k>th roots of <n> mod <qq>
 r, # one element of <rr>
 ss, # <k>th roots of <n> mod <q>
 s, # one element of <ss>
 tt; # temporary variable

 # get the arguments
 if Length(arg) = 2 then n := arg[1]; k := 2; m := arg[2];
 elif Length(arg) = 3 then n := arg[1]; k := arg[2]; m := arg[3];
 else Error("usage: RootsMod( <n>, <m> ) or RootsMod( <n>, <k>, <m> )");
 fi;
 Info( InfoNumtheor, 1, "RootsMod(", n, ",", k, ",", m, ")" );

 # check the arguments and reduce $n$ into the range $0..m-1$
 if m <= 0 then Error("<m> must be positive"); fi;
 n := n mod m;

 if not IsPrime(k) then
 # try over factors of k
 f:=Factors(k);
 l:=[n];
 for ii in f do
 l:=Concatenation(List(l,x->RootsMod(x,ii,m)));
 od;
 return l;
 fi;

 # combine the roots modulo every prime power $p^l$
 rr := [0]; qq := 1;
 for p in PrimeDivisors( m : UseProbabilisticPrimalityTest ) do

 # find prime power $q = p^l$
 q := p; l := 1;
 while m mod (q * p) = 0 do q := q * p; l := l + 1; od;

 # compute the roots mod $p^l$
 ss := RootsModPrimePower( n, k, p, l );

 # combine $rr$ (the roots mod $qq$) with $ss$ (the roots mod $p^l$)
 tt := [];
 ii := 1/qq mod q;
 for r in rr do
 for s in ss do
 Add( tt, r + qq * ((s-r)*ii mod q) );
 od;
 od;
 rr := tt;
 qq := qq * q;

 od;

 # return the roots $rr$
 Info( InfoNumtheor, 1, "RootsMod returns ", rr );
 return Set( rr );
end );

#############################################################################
##
#F RootsUnityMod( <m> ) . . . . . . . . . roots of unity modulo an integer
#F RootsUnityMod( <k>, <m> ) . . . . . . . roots of unity modulo an integer
##
BindGlobal( "RootsUnityModPrime", function ( k, p )
 local rr, # <k>th roots of 1 mod , result
 r, # <k>th root of unity mod 
 t, # <k>th root of unity mod 
 i; # loop variable

 # reduce $n$ into the range $0..p-1$
 Info( InfoNumtheor, 1, "RootsUnityModPrime(", k, ",", p, ")" );

 # handle $p = 2$
 if p = 2 then
 rr := [ 1 ];

 # it's easy if $k$ is invertible mod $p-1 = \phi(p)$
 elif GcdInt( p-1, k ) = 1 then
 rr := [ 1 ];

 # check that $n$ has a $k$th root (Euler's criterium)
 else

 # find a generator $t$ of the subgroup of $k$-th roots of 1.
 i:=2; t:=PowerModInt(i,(p-1)/k,p);
 while t=1 do
 i:=i+1; t:=PowerModInt(i,(p-1)/k,p);
 od;

 # combine $r$ (a particular root) with the powers of $t$
 r := 1;
 rr := [ 1 ];
 for i in [2..k] do
 r := r * t mod p;
 AddSet( rr, r );
 od;

 fi;

 # return the roots $rr$
 Info( InfoNumtheor, 1, "RootsUnityModPrime returns ", rr );
 return rr;
end );

DeclareGlobalName( "RootsUnityModPrimePower" );
BindGlobal( "RootsUnityModPrimePower", function ( k, p, l )
 local rr, # <k>th roots of <n> mod ^<l>, result
 r, # one element of <rr>
 ss, # <k>th roots of <n> mod smaller power
 s, # one element of <ss>
 t; # temporary variable

 # delegate prime case
 Info( InfoNumtheor, 1,
 "RootsUnityModPrimePower(", k, ",", p, "^", l, ")" );
 if l = 1 then
 rr := RootsUnityModPrime( k, p );

 # if $k$ is invertible mod $\phi(p^l)$ then there is only one root
 elif GcdInt(k,(p-1)*p) = 1 then
 rr := [ 1 ];

 # if $p = k = 2$
 elif p = 2 and k = 2 then
 rr := Set( [ 1, p^(l-1)-1, p^(l-1)+1, p^l-1 ] );

 # if $p = k$
 elif p = k then
 rr := [0..k-1]*p^(l-1)+1;

 # special case to speed up things a little bit
 elif k = 2 then
 rr := [ 1, p^l-1 ];

 # otherwise lift the roots with Newton / Hensel
 else

 # compute the root mod $p^{l/2}$
 ss := RootsUnityModPrimePower( k, p, QuoInt(l+1,2) );

 # lift the roots to $p^l$
 rr := [];
 for s in ss do
 Info( InfoNumtheor, 2, " lift root with Newton / Hensel" );
 t := PowerModInt( s, k-1, p^l );
 r := (s + (1 - t * s) / (k * t)) mod p^l;
 AddSet( rr, r );
 od;

 fi;

 # return the roots $rr$
 Info( InfoNumtheor, 1, "RootsUnityModPrimePower returns ", rr );
 return rr;
end );

InstallGlobalFunction( RootsUnityMod, function ( arg )
 local k, # <k>, optional first argument
 m, # <m>, second argument
 p, # prime divisor of <m>
 q, # power of 
 l, # <q> = ^<l>
 qq, # product of prime powers dividing <m>
 ii, # inverse of <qq> mod <q>
 rr, # <k>th roots of <n> mod <qq>
 r, # one element of <rr>
 ss, # <k>th roots of <n> mod <q>
 s, # one element of <ss>
 tt; # temporary variable

 # get the arguments
 if Length(arg) = 1 then k := 2; m := arg[1];
 elif Length(arg) = 2 then k := arg[1]; m := arg[2];
 else Error("usage: RootsUnityMod( <m> ) or RootsUnityMod( <k>, <m> )");
 fi;
 Info( InfoNumtheor, 1, "RootsUnityMod(", k, ",", m, ")" );

 # combine the roots modulo every prime power $p^l$
 rr := [0]; qq := 1;
 for p in PrimeDivisors( m ) do

 # find prime power $q = p^l$
 q := p; l := 1;
 while m mod (q * p) = 0 do q := q * p; l := l + 1; od;

 # compute the roots mod $p^l$
 ss := RootsUnityModPrimePower( k, p, l );

 # combine $rr$ (the roots mod $qq$) with $ss$ (the roots mod $p^l$)
 tt := [];
 ii := 1/qq mod q;
 for r in rr do
 for s in ss do
 Add( tt, r + qq * ((s-r)*ii mod q) );
 od;
 od;
 rr := tt;
 qq := qq * q;

 od;

 # return the roots $rr$
 Info( InfoNumtheor, 1, "RootsUnityMod returns ", rr );
 return Set( rr );
end );

#############################################################################
##
#F LogMod( <n>, <r>, <m> ) . . . . . . discrete logarithm modulo an integer
##
InstallGlobalFunction( LogModShanks,function(b,a,n)
local ai, m, m2, am, l, g, c, p, i;
 b:=b mod n;
 a:=a mod n;
 ai:=Gcdex(a,n);
 if ai.gcd=1 then
 # coprime case -- use Shanks's method
 # don't make the list longer than 5 million entries
 m:=Minimum(RootInt(n,2)+1,5*10^6);
 m2:=QuoInt(n,m)+1; # remaining part
 # calculate a^m mod n once
 am:=PowerMod(a,m,n);
 l:=[0..m-1];
 g:=[];
 c:=1;
 # create powers of a^m
 for i in l do
 Add(g,c);
 c:=c*am mod n;
 od;
 # dort the list (we'll potentially have to search often)
 SortParallel(g,l);
 c:=b;
 ai:=ai.coeff1;
 for i in [0..m2] do
 p:=PositionSorted(g,c);
 # positionsorted gives position to insert -- so we have to check
 if p<=m and g[p]=c then
 return l[p]*m+i;
 fi;
 c:=c*ai mod n;
 od;
 return fail;
 else
 Error("not coprime");
 fi;
end);

# Pollard Rho method for Index.
# Implemented by Sean Gage and AH

BindGlobal("LogModRhoIterate",function(n,g,p)
local p3, zp3, q, x, xd, a, ad, b, bd, m, r;
 p3:=QuoInt(p,3);
 zp3:=QuoInt(2*p,3);
 q := p-1;
 x := 1;
 xd := 1;
 a := 0;
 ad := 0;
 b := 0;
 bd := 0;
 repeat
 if x < p3 then
 x := (x * n) mod p;
 a := (a + 1) mod q;
 elif x < zp3 then
 x := (x * x) mod p;
 a := (a*2) mod q;
 b := (b*2) mod q;
 else
 x := (x * g) mod p;
 b := (b + 1) mod q;
 fi;
 if xd <p3 then
 xd := (xd * n) mod p;
 ad := (ad + 1) mod q;
 elif xd < zp3 then
 xd := (xd * xd) mod p;
 ad := (ad*2) mod q;
 bd := (bd*2) mod q;
 else
 xd := (xd * g) mod p;
 bd := (bd + 1) mod q;
 fi;
 if xd < p3 then
 xd := (xd * n) mod p;
 ad := (ad + 1) mod q;
 elif xd < zp3 then
 xd := (xd * xd) mod p;
 ad := (ad*2) mod q;
 bd := (bd*2) mod q;
 else
 xd := (xd * g) mod p;
 bd := (bd + 1) mod q;
 fi;
 until x=xd;

 m := (a-ad) mod q;
 r := (bd-b) mod q;
 return [m,r];
end);

InstallGlobalFunction(DoLogModRho,function(q,r,ord,f,p)
local fact, s, t, Q, R, MN, M, N, rep, d, k, theta, Qp,o,i;
 Info(InfoNumtheor,1,"DoLogModRho(",q,",",r,",",ord,",",p,")");
 fact:=[];
 s:=ord;
 for i in f do
 t:=s/i;
 if IsInt(t) then
 s:=t;
 Add(fact,i);
 fi;
 od;

 if Length(fact)>1 then
 d:=ord;
 while (d=ord) and Length(fact)>0 do
 s:=Remove(fact);
 t:=ord/s;
 Q:=PowerMod(q,s,p);
 R:=PowerMod(r,s,p);
 # iterate
 MN:=LogModRhoIterate(Q,R,p);
 M:=MN[1];
 N:=MN[2];
 rep:=GcdRepresentation(ord,s*M);
 d:=rep[1]*ord+rep[2]*s*M;
 od;
 if d<ord then
 k:=(rep[2]*s*N/d);
 if Gcd(DenominatorRat(k),ord)<>1 then
 return fail; # can't invert (can't happen if not primitive root)
 fi;
 k:=k mod ord;
 theta:=PowerMod(r,ord/d,p);
 Qp:=q/PowerMod(r,k,p) mod p;
 i:=DoLogModRho(Qp,theta,d,f,p);
 if i=fail then return i;fi; # bail out
 o:=(k+i*(ord/d)) mod ord;
 Assert(1,PowerMod(r,o,p)=q);
 return o;
 fi;
 fi;
 # naive case, iterate
 MN:=LogModRhoIterate(q,r,p);
 M:=MN[1];
 N:=MN[2];
 rep:=GcdRepresentation(ord,M);
 d:=rep[1]*ord+rep[2]*M;
 k:=(rep[2]*N/d);
 if Gcd(DenominatorRat(k),ord)<>1 then
 return fail; # can't invert (can't happen if not primitive root)
 fi;
 k:=k mod ord;
 theta:=PowerMod(r,ord/d,p);
 Qp:=q/PowerMod(r,k,p) mod p;
 for i in [1..d] do
 if Qp=1 then
 Assert(1,PowerMod(r,k,p)=q);
 return k;
 fi;
 k:=(k+ord/d) mod ord;
 Qp:=Qp/theta mod p;
 od;
 # process failed (because r was not a primitive root)
 return fail;
end);

InstallGlobalFunction( LogMod,function(b,a,n)
local c, p,f,l;

 b:=b mod n;
 a:=a mod n;
 if IsPrime(n) and Gcd(a,n)=1 then
 # use rho method
 f:=Factors(Integers,n-1:quiet); # Quick factorization, don't stop if its too hard
 l:=DoLogModRho(b,a,n-1,f,n);
 if l<>fail then
 return l;
 fi;
 fi;
 if Gcd(a,n)=1 then
 return LogModShanks(b,a,n);
 else
 # not coprime -- use old method
 c := 1;
 p := 0;
 while c <> b do
 c := (c * a) mod n;
 p := p + 1;
 if p = n then
 return fail;
 fi;
 od;
 return p;
 fi;
end);

#############################################################################
##
#M Sigma( <n> ) . . . . . . . . . . . . . . . sum of divisors of an integer
##
InstallMethod( Sigma,
 "sum of divisors of an integer",
 true, [ IsInt ], 0,

 function( n )

 local sigma, p, q, k;

 # make <n> it nonnegative, handle trivial cases
 if n < 0 then n := -n; fi;
 if n = 0 then Error("Sigma: <n> must not be 0"); fi;
 if n <= Length(DivisorsIntCache) then
 return Sum(DivisorsIntCache[n]);
 fi;

 # loop over all prime $p$ factors of $n$
 sigma := 1;
 for p in PrimeDivisors(n) do

 # compute $p^e$ and $k = 1+p+p^2+..p^e$
 q := p; k := 1 + p;
 while n mod (q * p) = 0 do q := q * p; k := k + q; od;

 # combine with the value found so far
 sigma := sigma * k;
 od;

 return sigma;
 end );

#############################################################################
##
#M Tau( <n> ) . . . . . . . . . . . . . . number of divisors of an integer
##
InstallMethod( Tau,
 "number of divisors of an integer",
 true, [ IsInt ], 0,

 function( n )

 local tau, p, q, k;

 # make <n> it nonnegative, handle trivial cases
 if n < 0 then n := -n; fi;
 if n = 0 then Error("Tau: <n> must not be 0"); fi;
 if n <= Length(DivisorsIntCache) then
 return Length(DivisorsIntCache[n]);
 fi;

 # loop over all prime factors $p$ of $n$
 tau := 1;
 for p in PrimeDivisors(n) do

 # compute $p^e$ and $k = e+1$
 q := p; k := 2;
 while n mod (q * p) = 0 do q := q * p; k := k + 1; od;

 # combine with the value found so far
 tau := tau * k;
 od;

 return tau;
 end );

#############################################################################
##
#F MoebiusMu( <n> ) . . . . . . . . . . . . . . Moebius inversion function
##
InstallGlobalFunction( MoebiusMu, function ( n )
 local factors;

 if n < 0 then n := -n; fi;
 if n = 0 then Error("MoebiusMu: <n> must be nonzero"); fi;
 if n = 1 then return 1; fi;

 factors := Factors(Integers, n );
 if factors <> Set( factors ) then return 0; fi;
 return (-1) ^ Length(factors);
end );

#############################################################################
##
#F TwoSquares( <n> ) . . . . . repres. of an integer as a sum of two squares
##
InstallGlobalFunction( TwoSquares, function ( n )
 local c, d, p, q, l, x, y;

 # check arguments and handle special cases
 if n < 0 then Error("<n> must be positive");
 elif n = 0 then return [ 0, 0 ];
 elif n = 1 then return [ 0, 1 ];
 fi;

 # write $n = c^2 d$, where $c$ has only prime factors $2$ and $4k+3$,
 # and $d$ has at most one $2$ and otherwise only prime factors $4k+1$.
 c := 1; d := 1;
 for p in PrimeDivisors( n ) do
 q := p; l := 1;
 while n mod (q * p) = 0 do q := q * p; l := l + 1; od;
 if p = 2 and l mod 2 = 0 then
 c := c * 2 ^ (l/2);
 elif p = 2 and l mod 2 = 1 then
 c := c * 2 ^ ((l-1)/2);
 d := d * 2;
 elif p mod 4 = 1 then
 d := d * q;
 elif p mod 4 = 3 and l mod 2 = 0 then
 c := c * p ^ (l/2);
 else # p mod 4 = 3 and l mod 2 = 1
 return fail;
 fi;
 od;

 # handle special cases
 if d = 1 then return [ 0, c ];
 elif d = 2 then return [ c, c ];
 fi;

 # compute a square root $x$ of $-1$ mod $d$, which must exist since it
 # exists modulo all prime powers that divide $d$
 x := RootMod( -1, d );

 # and now the Euclidean Algorithm strikes again
 y := d;
 while d < y^2 do
 p := x;
 x := y mod x;
 y := p;
 od;

 # return the representation
 return [ c * x, c * y ];
end );

InstallGlobalFunction(PValuation,function(n,p)
 if not IsInt(p) or not IsRat(n) or p = 0 then
 Error("wrong parameters");
 fi;
 if n = 0 then
 return infinity;
 elif IsInt(n) then
 return PVALUATION_INT(n,p);
 fi;
 return PVALUATION_INT(NumeratorRat(n),p) - PVALUATION_INT(DenominatorRat(n),p);
end);

Messung V0.5 in Prozent

¤ Dauer der Verarbeitung: 0.42 Sekunden (vorverarbeitet am 2026-10-13) ¤

Wurzel

Suchen

PVS Prover

Isabelle Prover

NIST Cobol Testsuite

Cephes Mathematical Library

Vienna Development Method

Haftungshinweis

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.