NumberTheory`NumberTheoryFunctions`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. This product of primes contains no squared factors.
Out[2]= |  |
The square number  divides  .
Out[3]= |  |
SquareFreeQ can handle large integers.
Out[4]= |  |
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  .
Out[5]= |  |
Even for large numbers, the next prime can be computed rather quickly.
Out[6]= |  |
This gives the largest prime less than  .
Out[7]= |  |
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  .
Out[8]= |  |
Operations involving prime factors. Here is the list of prime factors of  .
Out[9]= |  |
We can find the least prime factor directly.
Out[10]= |  |
PrimeFactorList operates on rational numbers.
Out[11]= |  |
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.
Out[12]= |  |
We can verify this by noting that the list of prime factors has a single element.
Out[13]= |  |
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  .
Out[14]= |  |
This confirms the result.
Out[15]= |  |
For longer lists the routine is still quite fast.
Out[16]= |  |
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  .
Out[17]= |  |
This verifies the result.
Out[18]= |  |
This returns all integers less than 11 which satisfy the relation.
Out[19]= |  |
If  does not have a square root modulo  , SqrtMod[d, n] will remain unevaluated. SqrtModList will return an empty list in this case.
PowerMod::"root": "The equation x^\!\(2\) = \!\(3\) (mod \!\(5\)) has no integer solutions."
Out[20]= |  |
This checks that 3 is not a square modulo 5.
Out[21]= |  |
Even for large modulus, the square root can be computed fairly quickly.
Out[22]= |  |
SqrtMod[d, n] also works for composite  .
Out[23]= |  |
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  .
Out[24]= |  |
This confirms that it does generate the group.
Out[25]= |  |
Here is a primitive root of a prime power.
Out[26]= |  |
Here is a primitive root of twice a prime power.
Out[27]= |  |
If the argument is composite and not a prime power or twice a prime power, the function does not evaluate.
Out[28]= |  |
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  .
Out[29]= |  |
This verifies the result.
Out[30]= |  |
The case  is essentially the same as FactorInteger[n, GaussianIntegers -> True]
Out[31]= |  |
Here is a fairly large composite number.
Out[32]= |  |
This computes its quadratic representation using  .
Out[33]= |  |
This verifies the result.
Out[34]= |  |
Even for large numbers, you can get a quadratic representation fairly quickly.
Out[35]= |  |
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  .
Out[36]= |  |
Here we verify that the input integer is a fundamental discriminant.
Out[37]= |  |
The class number is just the number of inequivalent quadratic forms with the given discriminant, if the discriminant is fundamental.
Out[38]= |  |
ClassNumber can accept positive discriminants, unlike ClassList.
Out[39]= |  |
Gauss conjectured and H. Stark proved in 1968 that this is the last negative discriminant of class number one.
Out[40]= |  |
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  .
Out[41]= |  |
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.
Out[42]= |  |
This checks that the representations are valid.
Out[43]= |  |
Here is the ordered list of the same representations.
Out[44]= |  |
The asymptotic average value of  is  .
Out[45]= |  |
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. 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  .
Out[47]= |  |
If SumOfFactors[n] equals  , then  is said to be perfect.
Out[48]= |  |
Here is the aliquot sequence starting at 95. Note that it is cyclic, terminating at  .
Out[49]= |  |
This returns the cyclic portion directly.
Out[50]= |  |
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.
Out[51]= |  |
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.
Out[52]= |  |
Unitary indicates only unitary divisors are to be added. A nontrivial factor  of  is unitary if  and  are coprime.
Out[53]= |  |
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  .
Out[54]= |  |
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  .
Out[55]= |  |
Exponential divisors use the same notation as Infinitary, but in this case either  or  divides  .
Out[56]= |  |
The ModifiedExponential case is similar to Exponential, except that here the condition for  to be a divisor is that  divides  and  divides  .
Out[57]= |  |
Powers of  become perfect in this case, since the sum of the nontrivial divisors of  is  .
Out[58]= |  |
| 
