**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 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, and E. Grosswald, Representations of Integers as Sums of Squares, Springer, 1985.

Testing for a squared factor.

SquareFreeQ[n] checks to see if n has a square prime factor. It is believed that this problem is, in general, as hard as the problem of finding the prime factorization of . However, in many cases one can do much better than simply factoring and seeing if a prime factor occurs with a power equal to two or higher.

To show that is not square-free, one only needs to find a single prime factor of whose square divides . In principle then, one could perform the factorization of one prime at a time to see if the square of the prime divides . This method, while simple to describe, is rather inefficient. The actual procedure used by SquareFreeQ[

n] is somewhat more complicated to describe, but it is much more efficient computationally.

First, the greatest common divisor of and the product of the primes less than is computed. (Call this product of primes .) The next step is to compute the GCD of this result and . If this GCD is not equal to one, then has a square factor. If this GCD is equal to one, then a new quantity is defined to be n

/GCD[n,P]. If is prime (an easy thing to check), then is square-free. If is not prime, then the process continues with the factorization of (the small factors found with the GCD no longer matter since they occur to power one in the factorization of ).

The first thing to check about is whether it is a perfect power of an integer (this can be done easily). If it is not, then the elliptic curve method of factorization is used as implemented in the package NumberTheory`FactorIntegerECM`. FactorIntegerECM is used to find a single factor of . The next step is to check whether and have no common factor. If they do, then has a square factor. If they don't, then the prime factorizations of and do not intersect, and the answer can be obtained recursively by applying SquareFreeQ to both and

.

This loads the package.
In[1]:= **<<NumberTheory`NumberTheoryFunctions`**

This product of primes contains no squared factors.
In[2]:= **SquareFreeQ[2*3*5*7]**

Out[2]=

The square number divides

.
In[3]:= **SquareFreeQ[60]**

Out[3]=

This took about two and a half hours on a workstation.
In[4]:= **SquareFreeQ[2^101 - 1]**

Out[4]=

SquareFreeQ uses FactorIntegerECM as a subroutine, so it should return a value for numbers less than 40 digits.

Finding the next prime.

NextPrime[n] finds the smallest 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]**

Out[5]=

Even for large numbers, the next prime can be computed rather quickly.
In[6]:= **NextPrime[10^64]**

Out[6]=

Solving simultaneous congruences.

The Chinese Remainder Theorem states that a certain class of simultaneous congruences always has a solution. ChineseRemainderTheorem[

,

] finds the smallest integer such that Mod[

r,

] is . This integer always exists if the numbers in are pairwise relatively prime. One way to prove this fact is to use the algorithm implemented by ChineseRemainderTheorem. The program in the package was taken from Stan Wagon's book, Mathematica in Action

, W.H. Freeman, 1991, Section 8.4.

This means that , , and

.
In[7]:= **ChineseRemainderTheorem[{0, 1, 2}, {4, 9, 121}]**

Out[7]=

This confirms the result.
In[8]:= **Mod[244, {4, 9, 121}]**

Out[8]=

For longer lists the routine is still quite fast.
In[9]:= **ChineseRemainderTheorem[Range[20], Prime[Range[20]]]**

Out[9]=

Finding square roots modulo n.

SqrtMod[d,n] computes the square root of . In other words, it returns , where . The algorithm used by SqrtMod is discussed in Mathematica in Action, Section 9.6. It is interesting to note that this method does fewer divisions than Newton's method for extracting square roots, and so is faster by a constant factor.

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.

If is not prime, the algorithm that is used requires that be factored completely. In this case, for to be a perfect square, every prime factor of for which JacobiSymbol[

d,p] is must divide to an even power.

In the special case when is a prime power , SqrtMod[

d,n] is computed using a variant of the classical Newton iteration

The Newton iteration sends a square root to a solution, but does a modular division at each step.

One can compute the square root by doing modular divisions only at the first step by first computing

and then iterating

as in Newton's iteration. The answer is recovered by multiplying back by on the final step.

The method is based on the

-adic identity

where , .

For a general composite ,

is first factored into primes. The square root modulo primes and prime powers is then computed using the techniques described above. These solutions are then used with the Chinese Remainder Theorem to get the final solution.

This finds an so that is equal to

.
In[10]:= **SqrtMod[3, 11]**

Out[10]=

This verifies the result.
In[11]:= **Mod[5^2, 11]**

Out[11]=

If does not have a square root modulo , Sqrt[

d,n] will remain unevaluated.
In[12]:= **SqrtMod[3, 5]**

Out[12]=

This checks that 3 is not a square modulo 5.
In[13]:= **Mod[{0, 1, 2, 3, 4}^2, 5]**

Out[13]=

Even for large modulus, the square root can be computed fairly quickly.
In[14]:= **SqrtMod[2, 10^64 + 57]**

Out[14]=

SqrtMod[d,n] also works for composite

.
In[15]:= **SqrtMod[3, 11^3]**

Out[15]=

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 single generator if and only if is a power of a prime or twice a power of a prime. The algorithm is probabilistic but runs very fast in practice (it does not require one to factor ). This implementation was improved by some suggestions of Ferrell S. Wheeler. A variation of this implementation can be used to generate the Pratt primality certificates used in the package NumberTheory`PrimeQ`. A discussion of Pratt certificates can be found in Mathematica in Action

, Section 8.7.

Here is the primitive root of

.
In[16]:= **PrimitiveRoot[5]**

Out[16]=

This confirms that it does generate the group.
In[17]:= **Sort[Mod[2^Range[4], 5]]**

Out[17]=

Here is the primitive root of a prime power.
In[18]:= **PrimitiveRoot[1093^3]**

Out[18]=

Here is the primitive root of twice a prime power.
In[19]:= **PrimitiveRoot[2*5^5]**

Out[19]=

If the argument is composite and not a prime power or twice a prime power, the function does not evaluate.
In[20]:= **PrimitiveRoot[11*13]**

Out[20]=

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 125-124. 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:

(i) is not a perfect square , that is, JacobiSymbol[-

d,n] is .

(ii) 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 (i) 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 (ii) is given in the book D.A. Cox, Primes of the Form

, Wiley, 1989.

This gives a quadratic representation of

.
In[21]:= **QuadraticRepresentation[1, 13]**

Out[21]=

This verifies the result.
In[22]:= **(%^2) . {1, 1}**

Out[22]=

The case is essentially the same as FactorInteger[

n,GaussianIntegers->True]
In[23]:= **FactorInteger[13, GaussianIntegers -> True]**

Out[23]=

Here is a fairly large composite number.
In[24]:= **13*31*61*Prime[10^7]**

Out[24]=

This computes its quadratic representation using

.
In[25]:= **QuadraticRepresentation[3, %]**

Out[25]=

This verifies the result.
In[26]:= **(%^2) . {1, 3}**

Out[26]=

Even for large numbers, you can get a quadratic representation fairly quickly.
In[27]:= **QuadraticRepresentation[1, 10^64+57]**

Out[27]=

Using ClassList.

ClassList[d] gives a set of representatives of the equivalence classes under composition of binary quadratic forms of discriminant . It is assumed that is a negative square-free integer of the form . 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

.

Here are representatives of the three quadratic forms with discriminant

.
In[28]:= **ClassList[-23]**

Out[28]=

The class number is just the number of inequivalent quadratic forms with the given discriminant.
In[29]:= **ClassNumber[-10099]**

Out[29]=

Gauss conjectured and H. Stark proved in 1968 that this is the last negative discriminant of class number one.
In[30]:= **ClassNumber[-163]**

Out[30]=

Representing an integer as a sum of squares.

SumOfSquaresRepresentations[d,n] gives a set of representations of the integer as a sum of 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[31]:= **SumOfSquaresRepresentations[3, 100]**

Out[31]=

This checks that the representations are valid.
In[32]:= **Apply[Plus, (%^2), 2]**

Out[32]=

The asymptotic average value of is

.
In[33]:= **Sum[N[SumOfSquaresR[2, k]], {k, 200}] / 200**

Out[33]=

SumOfSquaresRepresentations and SumOfSquaresR were contributed by Stan Wagon. For more information see "The Magic of Imaginary Factoring," Mathematica in Education and Research 5:1 (1996), pages 43-47.