This is documentation for Mathematica 5, which was
based on an earlier version of the Wolfram Language.
View current documentation (Version 11.2)

Documentation / Mathematica / Add-ons & Links / Standard Packages / NumberTheory /


This package contains a variety of functions that are useful for number theory applications. For more information on the mathematics related to these functions, see references such as G. H. Hardy and E. M. Wright, An Introduction to the Theory of Numbers, Oxford University Press, 1988, D. E. Knuth, Seminumerical Algorithms, Addison-Wesley, 1981, E. Grosswald, Representations of Integers as Sums of Squares, Springer, 1985, and H. Cohen, A Course in Computational Algebraic Number Theory, Springer-Verlag, 1993.

Testing for a squared factor.

SquareFreeQ[n] checks to see if has a square prime factor. This is done by computing MoebiusMu[n] and seeing if the result is zero; if it is, then is not squarefree, otherwise it is. Computing MoebiusMu[n] involves finding the smallest prime factor of . If has a small prime factor (less than or equal to ), this is very fast. Otherwise, FactorInteger is used to find .

This loads the package.

In[1]:= <<NumberTheory`NumberTheoryFunctions`

This product of primes contains no squared factors.

In[2]:= SquareFreeQ[2*3*5*7]


The square number divides .

In[3]:= SquareFreeQ[60]


SquareFreeQ can handle large integers.

In[4]:= SquareFreeQ[2^101 - 1]


Finding prime numbers.

NextPrime[n] finds the smallest prime such that , while PreviousPrime[n] finds the largest prime such that . For less than 20 digits, the algorithm does a direct search using PrimeQ on the odd numbers greater than . For with more than 20 digits, the algorithm builds a small sieve and first checks to see whether the candidate prime is divisible by a small prime before using PrimeQ. This seems to be slightly faster than a direct search.

This gives the next prime after .

In[5]:= NextPrime[10]


Even for large numbers, the next prime can be computed rather quickly.

In[6]:= NextPrime[10^100]//Timing


This gives the largest prime less than .

In[7]:= PreviousPrime[34]


For Random[Prime, range], a random prime is computed by first finding a random integer in that range and then computing the smallest prime , assuming is in that range, otherwise the process is repeated until a prime is found in the range. If no prime exists in the specified range, the input is returned unevaluated with an error message.

Here's a random prime between and .

In[8]:= Random[Prime, {10, 100}]


Operations involving prime factors.

Here is the list of prime factors of .

In[9]:= PrimeFactorList[713]


We can find the least prime factor directly.

In[10]:= LeastPrimeFactor[713]


PrimeFactorList operates on rational numbers.

In[11]:= PrimeFactorList[78/41]


The algorithm for PrimePowerQ involves first computing the least prime factor of and then attempting division by until either is obtained, in which case is a prime power, or until division is no longer possible, in which case is not a prime power.

Here is a number that is a power of a single prime.

In[12]:= PrimePowerQ[12167]


We can verify this by noting that the list of prime factors has a single element.

In[13]:= PrimeFactorList[12167]


Solving simultaneous congruences.

The Chinese Remainder Theorem states that a certain class of simultaneous congruences always has a solution. ChineseRemainder[, ] finds the smallest non-negative integer such that Mod[r, ] is . The solution is unique modulo the least common multiple of the elements of . The code for ChineseRemainder was contributed by Stan Wagon of Macalaster College.

This means that , , and .

In[14]:= ChineseRemainder[{0, 1, 2}, {4, 9, 121}]


This confirms the result.

In[15]:= Mod[244, {4, 9, 121}]


For longer lists the routine is still quite fast.

In[16]:= ChineseRemainder[Range[20], Prime[Range[20]]]


Finding square roots modulo n.

SqrtMod[d, n] computes the square root of . In other words, it returns the least nonnegative integer , where , when exists.

For given and , there may not exist an integer with . Clearly must be a perfect square , so to have a solution one must have that JacobiSymbol[d, n] be equal to 1. This condition is also sufficient if is a prime. The algorithm used for the case when is prime was discovered by Shanks.

Note that the square root need not be unique. SqrtMod returns the least nonnegative integer less than whose square is congruent to . SqrtModList returns all such integers that satisfy this relationship.

This finds the smallest nonnegative integer so that is equal to .

In[17]:= SqrtMod[3, 11]


This verifies the result.

In[18]:= Mod[5^2, 11]


This returns all integers less than 11 which satisfy the relation.

In[19]:= SqrtModList[3,11]


If does not have a square root modulo , SqrtMod[d, n] will remain unevaluated. SqrtModList will return an empty list in this case.

In[20]:= SqrtMod[3, 5]


This checks that 3 is not a square modulo 5.

In[21]:= Mod[{0, 1, 2, 3, 4}^2, 5]


Even for large modulus, the square root can be computed fairly quickly.

In[22]:= SqrtMod[2, 10^64 + 57]//Timing


SqrtMod[d, n] also works for composite .

In[23]:= SqrtMod[3, 11^3]


Since the algorithm factors , SqrtMod[d, n] may not return a result for very large composite values of .

Computing primitive roots.

PrimitiveRoot[n] returns a generator for the group of numbers relatively prime to under multiplication . This has a generator if and only if is 2, 4, a power of an odd prime, or twice a power of an odd prime. The algorithm is deterministic and uses FactorInteger. If is a prime or prime power, the least positive primitive root will be returned.

Here is a primitive root of .

In[24]:= PrimitiveRoot[5]


This confirms that it does generate the group.

In[25]:= Sort[Mod[2^Range[4], 5]]


Here is a primitive root of a prime power.

In[26]:= PrimitiveRoot[1093^3]


Here is a primitive root of twice a prime power.

In[27]:= PrimitiveRoot[2*5^5]


If the argument is composite and not a prime power or twice a prime power, the function does not evaluate.

In[28]:= PrimitiveRoot[11*13]


PrimitiveRoot uses FactorInteger as a subroutine, so it may not return a result for very large arguments.

Finding quadratic representations.

QuadraticRepresentation[d, n] returns x, y where , if such a representation exists. Here must be a positive integer and is odd. The algorithm resembles the Euclidean algorithm, and for the case of prime was discovered by Cornacchia, see S. Wagon, "The Euclidean Algorithm Strikes Again," American Mathematical Monthly 97 (1990), pages 124-125. The generalization to all is given in the paper by K. Hardy, J. B. Muskat, and K. S. Williams, "A Deterministic Algorithm for Solving in Coprime Integers and ," Mathematics of Computation 55 (1990), pages 327-343. This algorithm uses SqrtMod[-d, n] as a subroutine and therefore requires the factorization of .

QuadraticRepresentation[d, n] may not return an answer for two reasons:

(1)  is not a perfect square , that is, JacobiSymbol[-d, n] is .

(2) The class number of the extension field is greater than one, and a prime divisor of splits into prime ideals not in the principal class.

The reason that condition (1) can imply the nonexistence of a representation is that the equation implies that , so that is a perfect square (here division is ). It follows from this that for such a representation to exist, each prime that divides and for which is not a perfect square must divide to an even power.

A complete analysis of condition (2) is given in the book D. A. Cox, Primes of the Form , Wiley, 1989.

This gives a quadratic representation of .

In[29]:= QuadraticRepresentation[1, 13]


This verifies the result.

In[30]:= (%^2) . {1, 1}


The case is essentially the same as FactorInteger[n, GaussianIntegers -> True]

In[31]:= FactorInteger[13, GaussianIntegers -> True]


Here is a fairly large composite number.

In[32]:= 13*31*61*Prime[10^7]


This computes its quadratic representation using .

In[33]:= QuadraticRepresentation[3, %]


This verifies the result.

In[34]:= (%^2) . {1, 3}


Even for large numbers, you can get a quadratic representation fairly quickly.

In[35]:= QuadraticRepresentation[1, 10^64+57]//Timing


Using ClassList.

An integer is a fundamental discriminant if it is either congruent to where is odd and square-free, or it is congruent to and is even and is square-free.

ClassList[d] gives a set of representatives of the equivalence classes under composition of binary quadratic forms of negative integer . This list will have elements for square-free and of the form or . A quadratic form is represented by the list a, b, c. The algorithm is the most straightforward one and is slow for inputs much larger in magnitude than .

ClassNumber[d] finds the class number of the quadratic number field where is a fundamental discriminant. This is equal to the length of ClassList[d] for negative . ClassNumber uses an analytic formula which is asymptotically faster than a brute force search.

Here are representatives of the three quadratic forms with discriminant .

In[36]:= ClassList[-23]


Here we verify that the input integer is a fundamental discriminant.

In[37]:= FundamentalDiscriminantQ[-10099]


The class number is just the number of inequivalent quadratic forms with the given discriminant, if the discriminant is fundamental.

In[38]:= ClassNumber[-10099]


ClassNumber can accept positive discriminants, unlike ClassList.

In[39]:= ClassNumber[65]


Gauss conjectured and H. Stark proved in 1968 that this is the last negative discriminant of class number one.

In[40]:= ClassNumber[-163]


Kronecker symbols.

Kronecker symbols are extensions of Jacobi symbols to all integers , . They are used in operations involving discriminants of quadratic fields.

This finds the Kronecker symbol .

In[41]:= KroneckerSymbol[5, 3]


Representing an integer as a sum of squares.

SumOfSquaresRepresentations[d, n] gives a set of representations of the integer as a sum of squares, while OrderedSumOfSquaresRepresentations[d, n] gives only those representations as a sum of nondecreasing nonnegative squares. SumOfSquaresR[d, n] gives , the number of representations of the integer as a sum of squares. For , SumOfSquaresR[d, n] can handle large integer values of , as long as can be factored. For other values of , SumOfSquaresR[d, n] uses recursion, thus only modestly sized values of and can be used.

Here are the representations of as a sum of squares.

In[42]:= SumOfSquaresRepresentations[3, 100]


This checks that the representations are valid.

In[43]:= Apply[Plus, (%^2), 2]


Here is the ordered list of the same representations.

In[44]:= OrderedSumOfSquaresRepresentations[3, 100]


The asymptotic average value of is .

In[45]:= Sum[N[SumOfSquaresR[2, k]], {k, 200}] / 200


SumOfSquaresRepresentations and SumOfSquaresR are modifications of algorithms contributed by Stan Wagon.

Identifying de Moivre numbers.

Solutions to the cylcotomic equation are algebraic numbers lying on the unit circle forming the vertices of a regular polygon. They are called de Moivre numbers, or roots of unity. They can be expressed in the form , where and designate which root the number corresponds to. WhichRootOfUnity returns the pair n, k.

This finds that the input is a de Moivre corresponding to the first sixth root of unity.

In[46]:= WhichRootOfUnity[(1 + I Sqrt[3])/2]

Note that WhichRootOfUnity may not recognize a root of unity if it is not simplified enough.

Functions for working with aliquot sequences.

Given an integer , we can find the positive divisors which are less than and sum them. For example, has divisors , , and , which sum to . The function SumOfFactors computes this sum directly. Now, suppose that SumOfFactors is applied recursively to the sum? The resulting sequence of numbers is called an aliquot sequence. This sequence will terminate when the sum reaches . The sequence may not terminate, but become cyclic; the function AliquotCycle returns the cyclic portion. It has been theorized that all aliquot sequences must terminate or become cyclic, though this is unproven, and some evidence suggests that it is false. The smallest number for which we do not yet know whether the sequence terminates is , as of April 2003.

Here is the sum of divisors of .

In[47]:= SumOfFactors[14]


If SumOfFactors[n] equals , then is said to be perfect.

In[48]:= SumOfFactors[6]


Here is the aliquot sequence starting at 95. Note that it is cyclic, terminating at .

In[49]:= AliquotSequence[95]


This returns the cyclic portion directly.

In[50]:= AliquotCycle[95]


If the cycle is of length greater than , it is called a amicable chain, and its elements are called sociable numbers. If the cycle is of length , the elements are called amicable numbers. is the smallest amicable number.

Here is an amicable chain of length , discovered by Poulet in 1918.

In[51]:= AliquotCycle[12496]


Options for aliquot functions.

A variety of options can be used to modify the computations in AliquotSequence and AliquotCycle. The MaxIterations and MaxTerms options can limit the size computations for the aliquot sequence, by specifying the maximum number of terms the sequence should have and the size of the terms. ShowProgress can be set to True to monitor the progress of the computation. TermIncrement can be used to add an extra increment to the result of computing the value for each term in the sequence.

The SumOfFactorsType option applies to SumOfFactors as well as the aliquot functions. It determines what definition to use for the divisors. By default, all divisors are summed. Other definitions can be applied by setting this option to Unitary, Biunitary, Exponential, ModifiedExponential, or Infinitary.

Here is the ordinary sum of divisors with default behavior.

In[52]:= SumOfFactors[12]


Unitary indicates only unitary divisors are to be added. A nontrivial factor of is unitary if and are coprime.

In[53]:= SumOfFactors[12, SumOfFactorsType -> Unitary]


Biunitary denotes nontrivial divisors for which if for every prime factor of , is not the largest power of dividing where is the largest power of dividing .

In[54]:= SumOfFactors[12, SumOfFactorsType -> Biunitary]


Here, an infinitary divisor of is one where if for every prime factor of , if is the largest power of dividing and is the largest power of dividing , then all binary digits of that are zero are also zero for .

In[55]:= SumOfFactors[12, SumOfFactorsType -> Infinitary]


Exponential divisors use the same notation as Infinitary, but in this case either or divides .

In[56]:= SumOfFactors[12, SumOfFactorsType -> Exponential]


The ModifiedExponential case is similar to Exponential, except that here the condition for to be a divisor is that divides and divides .

In[57]:= SumOfFactors[12, SumOfFactorsType -> ModifiedExponential]


Powers of become perfect in this case, since the sum of the nontrivial divisors of is .

In[58]:= AliquotSequence[16, TermIncrement -> 1]