User documentation
Generalities
The functions in the NumTheory
file are predominantly basic
operations from number theory. Most of the functions may be
applied to machine integers or big integers (i.e. values of
type BigInt
). Please recall that computational number theory is
not the primary remit of CoCoALib, so do not expect to find a
complete collection of operations here -- you would do better to
look at Victor Shoup's NTL (Number Theory Library), or PARI/GP,
or some other specialized library/system.
See also BigIntOps
for very basic arithmetic operations on integers,
and BigRat
for very basic arithmetic operations on rational numbers.
Examples
The Functions Available For Use
Several of these functions give errors if they are handed unsuitable values:
unless otherwise indicated below the error is of type ERR::BadArg
.
All functions expecting a modulus will throw an error if the modulus is
less than 2 (or an unsigned long
value too large to fit into a long
).
GCD, LCM, etc.
The main functions available are:
gcd(m,n)
computes the non-negative gcd ofm
andn
. If both args are machine integers, the result is of typelong
(or error if it does not fit); otherwise result is of typeBigInt
.IsCoprime(m,n)
returnstrue
iffabs(gcd(m,n)) == 1
ExtGcd(a,b,m,n)
computes the non-negative gcd ofm
andn
; also setsa
andb
so thatgcd = a*m+b*n
. Ifm
andn
are machine integers thena
andb
must be of type (signed)long
. Ifm
andn
are of typeBigInt
thena
andb
must also be of typeBigInt
. The cofactorsa
andb
satisfyabs(a) <= abs(n)/2g
andabs(b) <= abs(m)/2g
whereg
is the gcd (inequalities are strict if possible). Error ifm=n=0
.InvMod(r,m)
computes the least positive inverse ofr
modulom
; throws error (ERR::DivByZero
) if the inverse does not exist. Throws error (ERR::BadModulus
) ifm < 2
(or too big forlong
). Result is of typelong
ifm
is a machine integer; otherwise result is of typeBigInt
.InvMod(r,m, RtnZeroOnError)
same asInvMod(r,m)
except that it returns 0 if the inverse does not exist;RtnZeroOnError
comes from an enum.InvModNoArgCheck(r,m)
computes the least positive inverse ofr
modulom
; ASSUMES0 <= r < m
and2 <= m <= MaxLong
; result is along
. Throws errorERR::DivByZero
ifgcd(r,m)
is not 1.lcm(m,n)
computes the non-negative lcm ofm
andn
. If both args are machine integers, the result is of typelong
; otherwise result is of typeBigInt
. Gives errorERR::ArgTooBig
if the lcm of two machine integers is too large to fit into along
.
There is a class called CoprimeFactorBasis_BigInt
for computing a coprime
factor basis of a set of integers:
CoprimeFactorBasis_BigInt()
default ctor; base is initially empty.CFB.myAddInfo(n)
use also the integern
when computing the factor base.CFB.myAddInfo(v)
use also the elements ofstd::vector<long> v
orstd::vector<BigInt> v
when computing the factor base.FactorBase(CFB)
returns the factor base obtained so far (asvector<BigInt>
).
Prime generation and tests
These functions are in NumTheory-prime
. The functions
NextPrime
, PrevPrime
and RandomSmallPrime
each produce a
result of type SmallPrime
(essentially a long
which is known
to be prime).
eratosthenes(n)
buildvector<bool>
sieve of Eratosthenes up ton
; entryk
corresponds to integer2*k+1
; max valid index isn/2
EratosthenesRange(lo, hi)
buildvector<bool>
sieve of Eratosthenes fromlo
up tohi
; iflo
is odd, it is replaced bylo+1
; similarly forhi
. In returned vector entryk
corresponds to integerlo+2*k
; max valid index is(hi-lo)/2
IsPrime(n)
tests the positive numbern
for primality (may be very slow for larger numbers). Gives error ifn <= 0
.IsProbPrime(n)
tests the positive numbern
for primality (fairly fast for large numbers, but in very rare cases may falsely declare a number to be prime). Gives error ifn <= 0
.IsProbPrime(n,iters)
tests the positive numbern
for primality; performsiters
iterations of the Miller-Rabin test (default value is 25). Gives error ifn <= 0
.NextPrime(n)
andPrevPrime(n)
compute next or previous positive prime (fitting into a machinelong
).NextPrime
returns 0 if no next "small" prime exists;PrevPrime
throwsOutOfRange
if arg is less than 3. Both throwBadArg
ifn < 0
.RandomSmallPrime(N)
-- generate a (uniform) random prime from 5 up toN
; error ifN < 5
orN >= 2^31
. Result is of typeSmallPrime
(essentially along
).NextProbPrime(N)
andPrevProbPrime(N)
compute next or previous positive probable prime (usesIsProbPrime
).PrevProbPrime
throwsOutOfRange
error if0 <= N < 3
. Both throwBadArg
error ifN < 0
.NextProbPrime(N,iters)
andPrevProbPrime(N,iters)
compute next or previous positive probable prime (usesIsProbPrime
with second argiters
).PrevProbPrime
throwsOutOfRange
error if0 <= N < 3
. Both throwBadArg
error ifN < 0
.
There are also iterators for generating primes (or almost primes) in increasing order:
PrimeSeq()
the sequence of primes starting with 2.PrimeSeqForCRT()
a sequence of primes starting with some "large" value, and going upwards.FastFinitePrimeSeq()
a sequence containing all primes up to some limit (much faster thanPrimeSeq
); limit is given by mem fnmyLastPrime
.FastMostlyPrimeSeq()
a sequence containing all primes and a few composites (much faster thanPrimeSeq
).NoSmallFactorSeq()
a sequence of positive integers which have no small factors.
If pseq
is one of these iterator objects then
*pseq
gives the current prime in the sequence (as a value of typeSmallPrime
orlong
forFastMostlyPrimeSeq
andNoSmallFactorSeq
)++pseq
advances 1 step along the sequenceIsEnded(pseq)
returnstrue
if the end of the sequence has been reachedCurrPrime(pseq)
same as*pseq
(only forPrimeSeq
andPrimeSeqForCRT
)NextPrime(pseq)
advances 1 step along the sequence, and returns the new "current prime" (only forPrimeSeq
andPrimeSeqForCRT
)
Factorization
factor(n)
finds the complete factorization ofn
(WARNING may be very slow for large numbers); NB implementation incompletefactor_TrialDiv(n,limit)
finds small prime factors ofn
(up to & including the specifiedlimit
); result is afactorization
object. Gives error iflimit
is not positive or too large to fit into along
. See alsoMultiplicityOf2
inBigIntOps
.factor_PollardRho(n,niters)
attempt to find a (single) factor ofn
(not nec. prime) using at mostniters
iterations; returns "empty" factorization if no factor was found; factor found may not be prime.SumOfFactors(n,k)
compute sum ofk
-th powers of positive factors ofn
SmallestNonDivisor(n)
finds smallest (positive prime) non-divisor ofn
; ifn=0
throwsERR::NotNonZero
.IsSqFree(n)
returnstrue
ifn
is square-free, otherwisefalse
; forBigInt
result is abool3
FactorMultiplicity(b,n)
find largestk
such thatpower(b,k)
dividesn
(error ifb < 2
orn
is zero) =====Pollard Rho Sequence=====PollardRhoSeq(N, start, incr)
create a sequence starting fromstart
and with incrementincr
PRS.myAdvance(k)
advance the sequence byk
stepsGetFactor(PRS)
returns a factor ofN
(may be 1 orN
)GetNumIters(PRS)
returns number of steps/iters performed
Other Functions on Integers
EulerPhi(n)
computes Euler's totient function of the positive numbern
(i.e. the number of integers up ton
which are coprime ton
, or the degree of then
-th cyclotomic polynomial). Gives error ifn <= 0
.PrimitiveRoot(p)
computes the least positive primitive root for the positive primep
. Gives error ifp
is not a positive prime. WARNING May be very slow for largep
(because it must factorizep-1
).KroneckerSymbol(res,mod)
(test ifres
is a quadratic residue) computes the Kronecker symbol, generalization of Jacobi symbol, generalization of Legendre symbolMultiplicativeOrderMod(res,mod)
computes multiplicative order ofres
modulomod
. Throws errorERR::BadArg
ifmod < 2
orgcd(res,mod)
is not 1.PowerMod(base,exp,modulus)
computesbase
to the powerexp
modulomodulus
; result is least non-negative residue. Ifmodulus
is a machine integer then the result is of typelong
(or error if it does not fit), otherwise the result is of typeBigInt
. Gives error ifmodulus <= 1
. GivesERR::DivByZero
ifexp
is negative andbase
cannot be inverted. Ifbase
andexp
are both zero, it produces 1.BinomialRepr(N,r)
produces the repr ofN
as a sum of binomial coeffs with "denoms"r, r-1, r-2, ...
Functions on Rationals
SimplestBigRatBetween(A,B)
computes the simplest rational betweenA
andB
SimplestBinaryRatBetween(A,B)
computes the simplest binary rational betweenA
andB
; result is a rational of form N*2^k where the integer N is minimal.
Continued Fractions
Several of these functions give errors if they are handed unsuitable values:
unless otherwise indicated below the error is of type ERR::BadArg
.
Recall that any real number has an expansion as a continued fraction (e.g. see Hardy & Wright for definition and many properties). This expansion is finite for any rational number. We adopt the following conventions which guarantee that the expansion is unique:
- the last partial quotient is greater than 1 (except for the expansion of integers <= 1)
- only the very first partial quotient may be non-positive.
For example, with these conventions the expansion of -7/3 is (-3, 1, 2).
The main functions available are:
ContFracIter(q)
constructs a new continued fraction iterator objectIsEnded(CFIter)
true iff the iterator has moved past the last partial quotientIsFinal(CFIter)
true iff the iterator is at the last partial quotientquot(CFIter)
gives the current partial quotient as aBigInt
(or throwsERR::IterEnded
)*CFIter
gives the current partial quotient as aBigInt
(or throwsERR::IterEnded
)++CFIter
moves to next partial quotient (or throwsERR::IterEnded
)ContFracApproximant()
for constructing a rational from its continued fraction quotientsCFA.myAppendQuot(q)
appends the quotientq
to the continued fractionCFA.myRational()
returns the rational associated to the continued fractionCFApproximantsIter(q)
constructs a new continued fraction approximant iteratorIsEnded(CFAIter)
true iff the iterator has moved past the last "partial quotient"*CFAIter
gives the current continued fraction approximant as aBigRat
(or throwsERR::IterEnded
)++CFAIter
moves to next approximant (or throwsERR::IterEnded
)CFApprox(q,eps)
gives the simplest cont. frac. approximant toq
with relative error at mosteps
Chinese Remaindering -- Integer Reconstruction
CoCoALib offers the class CRTMill
for reconstructing an integer from
several residue-modulus pairs via Chinese Remaindering. At the moment the
moduli from distinct pairs must be coprime.
The operations available are:
CRTMill()
ctor; initially the residue is 0 and the modulus is 1CRT.myAddInfo(res,mod)
give a new residue-modulus pair to theCRTMill
(error ifmod
is not coprime to all previous moduli)CRT.myAddInfo(res,mod,CRTMill::CoprimeModulus)
give a new residue-modulus pair to theCRTMill
asserting thatmod
is coprime to all previous moduli --CRTMill::CoprimeModulus
is a constantCombinedResidue(CRT)
the combined residue with absolute value less than or equal toCombinedModulus(CRT)/2
CombinedModulus(CRT)
the product of the moduli of all pairs given to the mill
Rational Reconstruction
CoCoALib offers two heuristic methods for reconstructing rationals from residue-modulus pairs; they have the same user interface but internally one algorithm is based on continued fractions while the other uses lattice reduction. The methods are heuristic, so may (rarely) produce an incorrect result.
NOTE the heuristic requires the (combined) modulus to be a certain amount larger than strictly necessary to reconstruct the correct answer (assuming perfect bounds are known). In practice, this means that the methods always fail if the combined modulus is too small.
The constructors available are:
RatReconstructByContFrac()
ctor for continued fraction method mill log-epsilon equal to 20RatReconstructByContFrac(LogEps)
ctor for continued fraction method mill with given log-epsilon (must be at least 3)RatReconstructByLattice(SafetyFactor)
ctor for lattice method mill with givenSafetyFactor
(0 --> use default)
The operations available are:
reconstructor.myAddInfo(res,mod)
give a new residue-modulus pair to the reconstructorIsConvincing(reconstructor)
givestrue
iff the mill can produce a convincing resultReconstructedRat(reconstructor)
gives the reconstructed rational (or an error ifIsConvincing
is not true).BadMFactor(reconstructor)
gives the "bad factor" of the combined modulus.
There is also a function for deterministic rational reconstruction which requires certain bounds to be given in input. It uses the continued fraction method.
RatReconstructWithBounds(e,P,Q,res,mod)
wheree
is upper bound for number of "bad" moduli,P
andQ
are upper bounds for numerator and denominator of the rational to be reconstructed, and(res[i],mod[i])
is a residue-modulus pair with distinct moduli being coprime.
Maintainer Documentation
- Correctness of
ExtendedEuclideanAlg
is not immediately clear, because the cofactor variables could conceivably overflow -- in fact this cannot happen (at least on a binary computer): for a proof see Shoup's book A Computational Introduction to Number Theory and Algebra, in particular Theorem 4.3 and the comment immediately following it. There is just one line where a harmless "overflow" could occur -- it is commented in the code. - I have decided to make
ExtGcd
give an error if the inputs are both zero because this is an exceptional case, and so should be handled specially. I note thatmpz_exgcd
accepts this case, and returns two zero cofactors; so if we decide to accept this case, we should do the same -- this all fits in well with the (reasonable/good) principle that "zero inputs have zero cofactors". - Several functions are more complicated than you might expect because I wanted them to be correct for all possible machine integer inputs (e.g. including the most negativelong
value). - In some cases the function which does all the work is implemented as a file
local function operating on
unsigned long
values: the function should normally be used only via the "dispatch" functions whose args are of typeMachineInt
orBigInt
.
- The fns for primes (testing and generating) are in the file
NumTheory-prime
. - The impl of
eratosthenes
is fairly straightforward given that I chose to represent only odd numbers in the table: thek
-th entry corresponds to the integer2*k+1
. Overflow cannot occur: recall that the table size is at most half the biggestlong
. I'm hoping that C++11 will avoid the cost of copying the result upon returning. Anyway, I think the code is a decent compromise between readability, speed and space efficiency. The impl ofEratosthenesRange
is similar but the table covers just the given range (only odd numbers are represented; index 0 corresponds to smallest odd integer greater than or equal to the start of the range).
- The "prime sequence" classes are a bit messier than I'd like.
It was delicate getting correct the switch-over from one technique to the
other (in those classes where 2 techniques were used). The limit of 23
for
NoSmallFactorSeq
is somewhat arbitrary. Not sure the code is 32-bit safe.
- The continued fraction functions are all pretty simple. The only tricky
part is that the "end" of the
ContFracIter
is represented by bothmyFrac
andmyQuot
being zero. This means that a newly created iterator for zero is already ended. CFApproximantsIter
delegates most of the work toContFracIter
.
Bugs, Shortcomings, etc.
- Several functions return
long
values when perhapsunsigned long
would possibly be better choice (since it offers a greater range, and in the case ofgcd
it would permit the fn to return a result always, rather than report "overflow"). The choice of return type was dictated by the coding conventions, which were in turn dictated by the risks of nasty surprises to unwary users unfamiliar with the foibles of unsigned values in C++. NextPrime
has dodgy semantics: what happens when the end of the iterator is reached? In fact, all the non-finite "prime seq" iterators do not handle end-of-iterator properly!- Should there also be procedural forms of functions which return
BigInt
values? (e.g. gcd, lcm, InvMod, PowerMod, and so on). (2016-06-27 this will probably become irrelevant when using "move" semantics in C++11). - Certain implementations of
PowerMod
should be improved (e.g. to usePowerModSmallModulus
whenever possible). Is behaviour for 0^0 correct? KroneckerSymbol
I have chosen to make available justKroneckerSymbol
rather than the more widely-knownLegendreSymbol
because GMP makes Kronecker available, and it is always defined; whereasLegendreSymbol
would have to check that its 2nd arg is a prime (which would dominate the cost of the call)LucasTest
should produce a certificate, and be made publicly accessible.- How should the cont frac iterators be printed out???
ContFracIter
could be rather more efficient for rationals having very large numerator and denominator. One way would be to compute with num and den divided by the same large factor (probably a power of 2), and taking care to monitor how accurate these "scaled" num and den are. I'll wait until there is a real need before implementing (as I expect it will turn out a bit messy).CFApproximantsIter::operator++()
should be made more efficient.
Main changes
2022
- Feb (v0.99720):
SmoothFactor
has been renamedfactor_TrialDiv
- added
factor_PollardRho
-