Discovering Number Theory

Jeffrey J. Holt
University of Virginia

John W. Jones
Arizona State University

Copyright © 2001 by W. H. Freeman and Company

Divisibility and Factorization

Initializations

\(Off[General::spell]; 
Off[General::spell1]; \)

Divisors and GCDs

In the Prelab section of this chapter, you had the opportunity to compute the divisors of some integers by hand.  Mathematica has a built-in command called Divisors that will automatically compute all of the positive divisors for a given integer.  Here it is in action:

Divisors[12345]

The output is the set of all the positive divisors of 12345, given in the form of a Mathematica list.  Lists in Mathematica are contained in curly braces, { }.  Mathematica only lists the positive divisors because it is easy to compute all of the divisors from these.  We would just have to take these numbers and their negations to get all of the divisors of 12345.

In the Prelab section, an illustration was given that showed how to compute gcd(45, 60).  The first step is to compute the positive divisors of both 45 and 60.  Here's the Mathematica version:

Divisors[45]Divisors[60]

The next step is to pick out the largest divisor that occurs in both lists.  A handy way to do this is by using the Mathematica Intersection command.  This command takes two (or more) lists as input, and produces a new sorted list of all elements common to all of the input lists.  In our case, we will get a list of all positive divisors of both 45 and 60.  Here it is:

Intersection[Divisors[45], Divisors[60]]

With no work at all it is now clear that gcd(45, 60) = 15.  The above example is cool, but after all, you could do that one by hand.  Let's try some larger choices to really see the power of Mathematica:

\(a = 223092870; 
\(b = 6227020800; 

    Intersection[Divisors[a], Divisors[b]]\)\)

Imagine trying to do this one by hand!  By the way, what is gcd(a, b)?  Then use Mathematica to compute the positive divisors of gcd(a, b).  (You have to replace ??? with gcd(a, b) in the next cell before executing it.)

Divisors[???]

Research Question 1

In the preceding computations, we have found the set of common positive divisors of a = 223092870 and b  =  6227020800, and have found the set of positive divisors of gcd(a, b).  Compare these two sets.  Change the values of a and b, and repeat the experiment: compute the set of common positive divisors of a and b, compute the set of positive divisors of gcd(a, b), and compare the two sets.  Once you think you see a pattern, state a conjecture and try to prove it!

Hints:  Since your conjecture will involve sets, when working on the proof you should keep in mind the suggestions concerning proofs involving sets in the "Tips for Writing Proofs" section of the Course Guide.  Also, you may find it helpful to use the results in the "GCDs and Factorization" section below when working on your proof.

Divisors and Factorization

How do we find the set of positive divisors for an integer?  If you happen to have Mathematica at your disposal, you could just type in the Divisors command, as you did above.  However, it is important to understand how this works, especially if you ever stray more than a few feet from your computer.

The most basic method for computing divisors is exhaustive trial division.  If we want to find the positive divisors for an integer n, we just take the integers 1, 2, 3, . . . , n, divide n by each, and those that divide evenly make up the set of positive divisors for n.  Here is a function which does exactly that.

\(\(getdivs1[n_]\  := \ 
      Module[{divlist}, 
	divlist = {\ }; 
	
        Do[If[n/m\  == \ Floor[n/m], 
		AppendTo[divlist, m]], {m, 1, n}]; 
        
	Return[divlist]]; \)\)

The test for whether  m is a divisor of n in this function is n/m == Floor[n/m].  The == tests whether two quantities are equal.  On the basis of what you learned about the Floor command in the "Introduction to Using Mathematica", can you see why this works?  Let's try our new function:

getdivs1[12345]

We can compare it with the built-in Divisors command:

Divisors[12345]

Our method seems to work and is rather simple, but it is also quite inefficient.

One way to improve upon the above procedure is by applying the following observation: if m is a divisor of n, then k = n/m is also a divisor of n, because mk=n.  Thus, the positive divisors can be organized into pairs of the form (m, n/m), where m < n/m.  (The one exception to this is if n is a perfect square and m = (Sqrt[n]), in which case  m =  n/m.)  For example, if n = 100, then the positive divisors of n are given by

Divisors[100]

For this value of n, we have the pairs (1, 100),  (2, 50), (4, 25), and (5, 20), and (because 100 is a perfect square) the single divisor (Sqrt[100]) = 10.  Note that for each pair of positive divisors, the smaller divisor is less than (Sqrt[n]).  This is also the case in general: if m | n and m < n/m, then m < (Sqrt[n]).  (The proof of this assertion is left as a homework exercise.)  Thus, all positive divisors can be found by first using exhaustive trial division by all integers up to (Sqrt[n]), and then generating the pairs to identify the remaining positive divisors.

For example, suppose that n = 152.  Then we have

N[Sqrt[152]]

Dividing each of 1, 2, 3, . . . , 12 into 152 reveals the divisors 1, 2, 4, and 8.  The remaining positive divisors of 152 are then given by

152/1  =  152152/2  =  76152/4  =  38152/8  =  19

Here's a check:

Divisors[152]

Here is a function which uses this second method for finding all positive divisors of an integer n:

\(\(getdivs2[n_]\  := \ 
      Module[{divlist}, 
	divlist = {\ }; 
	
        Do[If[n/m\  == \ Floor[n/m], 
		AppendTo[divlist, {m, n/m}]], 

          	{m, 1, Sqrt[n]}]; 
	Return[Union[Flatten[divlist]]]]; \)\)

Here it is in action.

getdivs2[12345]

How much time do we save by using the second method instead of the first?  We can measure how long it takes each method to compute the divisors of 10000 using the Timing command, as shown:

Timing[getdivs1[10000]]

Timing[getdivs2[10000]]

The second method is much faster.  The reason for the difference in speed is clear: getdivs1 has to test n possible integers, while getdivs2 checks only (Sqrt[n]) integers.  So the running time for getdivs1 should be roughly proportional to n whereas the running time for getdivs2 should be roughly proportional to (Sqrt[n]).  (There are other things going on inside the computer, so one should not expect timing estimates like this to be exact.  However, they can provide good ballpark estimates.)

Experiment a little with each function.  What do you think will happen to the running times if we multiply n by 10 for each function?  Then check the result by changing 10000 in the commands above.

Now try to answer the following question.

Exercise 1

On the basis of the timings above, extrapolate the amount of time (in years) required to find the positive divisors for (10^18) using getdivs1 and getdivs2.  For simplicity, you may assume that all years have 365 days.

A problem related to finding the set of divisors for an integer is that of finding the prime-power factorization of an integer.  Indeed, we shall see later in the chapter how to use the factorization of an integer to easily compute all of the divisors for the integer.

A number of different methods have been developed for finding the factorization of an integer.  The method that we shall describe here is fairly simple, and borrows from the ideas used above to find divisors; that is, it uses trial division.

Recall that finding the prime factorization of an integer requires us to express an integer n as

\(TraditionalForm
    \`n\  = \ \(p\_1\%n\_1\) \(p\_2\%n\_2\) ⋯\ 
        \(p\_k\%n\_k . \)\)

The basic procedure for finding the factorization is as follows:

As an example, let us find the prime-power factorization for n = 571450.  We start by dividing n by 2:

571450/2

Although it is pretty obvious from the output, let's see if 2 will divide evenly into the above quotient:

285725/2

Since we get a fraction for output, we know that 2 won't divide in again.  Thus, we know that (2^1) is the correct power of 2 in the factorization.  Now let's move to the prime 3:

285725/3

This doesn't divide in evenly, so 3 is not included in the factorization.  Let's try 5:

285725/5

Thus 5 divides in evenly, and as we can see from the output, 5 will divide in evenly again:

57145/5

That's it for the 5s.  Hence we know that (5^2) is the correct power of 5 in the factorization.  Here's the test for the prime 7:

11429/7

No luck.  Here's the prime 11:

11429/11

Continuing in this manner, we find that 1039 is not divisible by 11, 13, 17, 19, 23, 29, or 31.  We already know that 1039 is not divisible by 2, 3, or 5 because we divided away all factors of 2, 3, and 5 already.  If we are clever, we can stop now since (Sqrt[1039]) is approximately . . .

N[Sqrt[1039]]

Because 1039 has no prime divisors less than or equal (Sqrt[1039]), it must be prime.  Thus we now have the entire factorization:

\(TraditionalForm
      \`571450 = 2\^1·5\^2·11\^1·1039\).

The method described and illustrated above will work well for fairly small integers and can be used to factor numbers on a hand calculator.  However, this technique is inefficient for larger integers, and more-advanced methods must be used in these cases.  But even the more-sophisticated factoring techniques have their limits.  The general problem of factoring large integers (upwards of 200 digits!) is a difficult one, and is currently an area of intense mathematical research.   The problem is not only interesting in its own right, but also has some important applications.  We'll see how factoring is related to cryptography in a later chapter.

GCDs and Factorization

Here's a computation we saw earlier:

\(a = 223092870; 
\(b = 6227020800; 

    Intersection[Divisors[a], Divisors[b]]\)\)

On the basis of the output, we see that gcd(a, b) = 30030.

Another way to compute gcd(a, b) is to look at the factorizations of a and b.  The Mathematica command FactorInteger computes factorizations of integers.  To illustrate, let's take a = 7920 and b = 4536.

\(a = 7920; 
\(b = 4536; \)\)FactorInteger[a]FactorInteger[b]

The output tells us that

a = 7920 = (2^4) (3^2) (5^1) (11^1)

and

b = 4536 = (2^3) (3^4) (7^1).

Now if d divides both a and b, then all of the prime divisors of d have to divide a and all of the prime divisors of d have to divide b.  The only primes that divide both a and b are 2 and 3, so any positive common divisor must have the form

d = (2^m) (3^n).

What are the largest choices of m and n we can take?  Well, on the basis of the factorizations,

if d | a then m ≤ 4, and if d | b then m ≤ 3.

So, we've got to take m ≤ 3.  Similarly,

if d | a then n ≤ 2, and if d | b then n ≤ 4.

This time, we've got to take n ≤ 2.  Thus, taking the biggest possible values for n and m, we see that the greatest common divisor is

(2^3) (3^2) = 72.

In general, suppose we have the prime-power factorizations

\(TraditionalForm
    \`a\  = \ \(p\_1\%m\_1\) \(p\_2\%m\_2\) ⋯\ p\_k\%m\_k\)

and

\(TraditionalForm
    \`b\  = \ \(p\_1\%n\_1\) \(p\_2\%n\_2\) ⋯\ 
        \(p\_k\%n\_k . \)\)

Note: It is possible that some of the (Subscript[m, 1]), . . . , (Subscript[m, k]) and (Subscript[n, 1]), . . . , (Subscript[n, k]) are equal to 0.  Writing two integers as products of prime powers over the same set of primes is a useful trick when writing proofs involving prime power factorizations.  You can use it yourself, as long as you advise the reader of what you are doing.

The greatest common divisor of a and b is

\(TraditionalForm
    \`gcd(a, \ b)\  = \ 
      \(p\_1\%\(min(m\_1, \ n\_1)\)\) p\_2\%\(min(m\_2, \ n\_2)\)\ 
        ⋯\ \ \(p\_k\%\(min(m\_k, \ n\_k)\) . \)\)

To see how this works in practice, let's look at another example.  Suppose that we have a = 165620000 and b = 65984625.  Here are the factorizations:

\(a = 165620000; 
b = 65984625; 
FactorInteger[a]\)FactorInteger[b]

To compute gcd(a, b), we just start comparing prime factors.  Since 2 is a factor of a but not b, 2 is not a factor of gcd(a, b).  Similarly, since 3 is a factor of b but not a, 3 is not a factor of gcd(a, b).  Now 5 is a factor of both a and b, and the smallest exponent on 5 is 3, so that (5^3) is a factor of gcd(a, b).  Also, 7 is a factor of both a and b, where the smallest exponent is 2, so that (7^2) is a factor of gcd(a, b).  The factors 13 and 19 are eliminated in the same manner as 2 and 3, and we arrive at

gcd(a, b) = (5^3)(7^2)= 6125.

Mathematica has a built-in command to compute greatest common divisors.  The name of this command is GCD.  We can use it to check our computation:

GCD[165620000,  65984625]

Looks good.

If you play around with Mathematica's GCD command, you will find that it computes gcds very quickly, even for huge numbers.  In fact, we can use it on numbers which are too large for Mathematica to factor:

GCD[10^300  +  1,  10^500  +  1]

How does it find gcds without factoring the numbers?  We will learn the secret in the next chapter.

Counting Divisors

By using a combination of the Divisors and Length commands, it's easy to get Mathematica to compute the number of positive divisors of a given integer.  Here's an example:

Length[Divisors[7920]]

All that this code does is compute the list of positive divisors, and then report the number of elements in the list.

In this section, we shall perform some experiments to learn how the number of positive divisors of a given integer can be determined from the factorization of the integer.  We'll do this by starting with fairly simple cases and building up to more complicated situations.  Along the way, a pattern should (hopefully!) emerge.

Let us begin with the simplest case: Suppose that n = p, where p is a prime number.  For instance, let n = 13.  Here's the positive divisors:

Divisors[13]

And here's d(n), the number of positive divisors of n:

Length[Divisors[13]]

Research Question 2

Select different primes, and repeat the above experiment: Compute the divisor list, and then compute d(n).  When you have enough data, form a conjecture that states d(n) for n = p, where p is prime.  (It probably won't take much data for you to form a conjecture!)  Then prove your conjecture.

Research Question 3

Repeat Research Question 2, but this time for n = (p^a), where p is a prime and a is a positive integer.  

Hints:  When doing computational experiments, select a small prime p and then increment the value of a by 1.  Also, note that when a = 1, you revert to the case covered by Research Question 2.  Thus, whatever your conjecture, it must be consistent with your results from the preceding Research Question.

Now suppose that n = pq, where p and q are distinct primes.  (If p = q, then we would be in the situation covered in Research Question 3.)  What is d(n) in this case?  Let's try an experiment, taking n = 3 · 7:

\(n = 3*7; 
Divisors[n]\)Length[Divisors[n]]

Research Question 4

Select different primes p and q, and repeat the above experiment: Compute the list of positive divisors for n = p q, and then compute d(n).  When you have enough data, form a conjecture that states d(n) for n = p q, where p and q are distinct primes.  

Note:  We will no longer state that you should prove your conjecture, since you should always prove (or try to prove) all conjectures made in this course.

Research Question 5

Repeat Research Question 4, but this time take n = (p^a)(q^b), where p and q are primes and a and b are positive integers.

Hints:  When doing your numerical experiments, it will make life easier if you are systematic.  (In fact, systematic experimentation is the hallmark of good scientific work.)  You might try starting with fixed values of p and q, and make simple variations in the exponents a and b.  Note also that the previous Research Questions are all special cases of this one.  Therefore, your conjecture in this case must be consistent with your previous results.

All of the earlier work has been leading up to the most general case of all.  Suppose that

\(TraditionalForm
    \`\(n\  = \ \(p\_1\%a\_1\) \(p\_2\%a\_2\) ⋯\ p\_k\%a\_k, 
    \)\)

where (Subscript[p, 1]), (Subscript[p, 2]), . . . , (Subscript[p, k]) are distinct primes and (Subscript[a, 1]), (Subscript[a, 2]), . . . , (Subscript[a, k]) are positive integers.  What is d(n) in this case?  On the basis of your earlier investigations, you may not need any additional experimentation to form a conjecture.  On the other hand, there is no harm in further experimentation, so feel free to do more if you wish.

Research Question 6

Form a conjecture that states d(n) for 

				n = (Subscript[(Subscript[p, 1]^a), 1])(Subscript[(Subscript[p, 2]^a), 2])⋯ (Subscript[(Subscript[p, k]^a), k]),
		
where (Subscript[p, 1]), (Subscript[p, 2]), . . . , (Subscript[p, k]) are distinct primes and (Subscript[a, 1]), (Subscript[a, 2]), . . . , (Subscript[a, k]) are positive integers.

It may seem that we led to this formula in an inefficient manner; why not go for the whole formula at once?  As you may have discovered along the way, it can help a great deal to work up to the general case through various special cases.  In number theory, starting with cases such as n prime, and then n a prime power, and so on, is a standard progression of study.  Keep this in mind for later investigations in the course.

An important part of proving your conjectures for Research Questions 2 through 6 is being able to list the positive divisors of an integer of the forms n = p, n = p q, n = (p^a), and so on.  The basic idea of listing these divisors was implicit in the section above where we derived a formula for the greatest common divisor of two integers in terms of their factorizations.  In the next exercise, you should formulate the statement precisely.

Exercise 2

Suppose \(TraditionalForm
      \`n\  = \ \(p\_1\%m\_1\) \(p\_2\%m\_2\) ⋯\ 
          p\_k\%m\_k\).  What are the positive divisors of n in terms of their prime-power factorizations?

The Division Algorithm

We begin this section with a statement of the Division Algorithm, which you saw at the end of the Prelab section of this chapter:

Theorem 1.2 (Division Algorithm)  Let a be an integer and b be a positive integer.  Then there exist unique integers q and r such that

a = b q + r	and	0 ≤ r < b.

Remember learning long division in grade school?  (If not, pretend that you do.)  Long division is a procedure for dividing a number a by another number b, and coming up with a quotient q and a remainder r.  The Division Algorithm is really nothing more than a guarantee that good old long division really works.  Although this result doesn't seem too profound, it is nonetheless quite handy.  For instance, it is used in proving the Fundamental Theorem of Arithmetic, and will also appear in the next chapter.  A proof of the Division Algorithm is given at the end of the "Tips for Writing Proofs" section of the Course Guide.

Now, suppose that you have a pair of integers a and b, and would like to find the corresponding q and r.  If a and b are small, then you could find q and r by trial and error.  However, suppose that a = 124389001 and b = 593.  In this case, trial and error would be extremely tedious.  What happens if we try the division using Mathematica?  Since q is the quotient, it should be close to a/b.

\(a = 124389001; 
b = 593; 
a/b\)

Maybe a decimal approximation would be more useful:

N[a/b, 10]

Now we're getting somewhere.  Since we know that a = b q + r, it follows that

(a/b) = q + (r/b).

Thus it follows that q ≤ a/b.  (Remember that 0 ≤ r < b.)  So, in our above example, it makes sense to take q = 209762, because this is the biggest integer that is less than (or equal to) a/b.  (Recall that this is the idea used repeatedly when performing long division by hand.)  With this choice of q, what is r?  That's easy enough, because we must have

a -  b q = r.

Here's the computation:

\(q = 209762; 
r = a - b\ q\)

Since r < b, we have fulfilled the requirements of the Division Algorithm, and we're done.

Let's return to the determination of q.  The procedure used in the above example will work in general: Find the largest integer that is less than or equal to a/b.  This process can be automated by using the Mathematica Floor command, which does exactly what is required.  So, we could have computed q this way:

q = Floor[a/b]

The whole process is automated in the command divalg defined below:

\(divalg[a_, b_] := 
      Module[{q, r}, 
		q\  = \ Floor[a/b]; 
		r\  = \ a\  - \ b\ q; 

        		Return[{q, r}]]\)

Note that this function is made up of two operations, separated by semicolons.  Here it is in action on our sample values of a and b:

divalg[124389001, 593]

Exercise 3

Use divalg to repeat Exercise 5 in the Prelab section of this chapter.