C++ For Mathematicians (2006) [eng]

Chapter 5

Arrays

In Chapter 3 we introduced the following problem. Let pn be the probability that two integers, chosen uniformly and independently from {1,2,...,n}, are relatively prime. What can we say about pn as n → ∞? We computed pn for various values of n by direct enumeration. As n approached 105, the program became slow and so we sought another approach.

In Chapter 4 we tried a Monte Carlo approach. Although our program enables us to estimate pn for n equal to one million (or larger), to get decent accuracy we need to run the program for too many iterations.

In this chapter we take another approach using Euler’s totient.

5.1Euler’s totient

Euler’s totient is a function ϕ defined on the positive integers. For n Z+, we define ϕ(n) to be the number of elements of {1,2,...,n} that are relatively prime to n. For example, ϕ(10) = 4 because the only integers between 1 and 10 (inclusive) that are relatively prime to 10 are 1, 3, 7, and 9.

Euler’s totient has immediate relevance to our problem. We want to count the number of pairs (a,b) with 1 ≤ a,b ≤ n and gcd(a,b) = 1. Alternatively, we can count the number of such pairs with a ≤ b. Of course, we miss the pairs with a > b, so we would just double the result and subtract one (for double-counting the pair (1,1)). Symbolically, we have

	n	"−	k=1	#
	1		n
pn =	2		1 + 2 ∑ ϕ(k) .

We can calculate ϕ(k) by considering all the integers from 1 to k and check which are relatively prime to k. This, however, would lead to an algorithm that is no different from the one in Chapter 3.

Fortunately, there are more efficient ways to calculate ϕ.

We begin our analysis of Euler’s totient with a few special cases.

If p is a prime, then all but one member of the set {1,2,..., p} are relatively prime to p. Therefore ϕ(p) = p −1.

67

68			C++ for Mathematicians
Next consider ϕ	(p2)	for a	prime p. In the set			{	1,2,..., p2	}	the only numbers
				2
that are not relatively prime to p					are the multiples of p, and there are p of those.

Therefore ϕ(p2) = p2 − p = p(p −1).

More generally, consider ϕ(pn) where p is prime and n Z+. In {1,2,..., pn} only the multiples of p are not relatively prime to pn, and there are pn−1 of those. This gives the following.

Proposition 5.1. Let p be a prime and let n be a positive integer. Then ϕ(pn) = pn − pn−1 = pn−1(p −1).

Now that we have examined ϕ for powers of primes, let us examine the special case ϕ(77). Note that for 1 ≤ n ≤ 77, we have

gcd(n,77) = 1		gcd(n,7) = 1 and		gcd(n,11) = 1.
By Proposition 3.1,
gcd(n,7) = gcd(7,n mod 7)		and	gcd(n,11) = gcd(11,n mod 11).
Let n1 = n mod 7 and n2 = n mod 11. If n is relatively prime to 77 and
x ≡ n1	(mod 7)	and	x ≡ n2	(mod 11)

then x is also relatively prime to 77. Furthermore, there is a unique such x in {1,2,...,77}; this is a consequence of the Chinese Remainder Theorem.

Theorem 5.2 (Chinese Remainder). Let a and b be relatively prime positive integers, and let c,d Z. Then the system of congruences

x ≡ c mod a x ≡ d mod b

has a unique solution in {1,2,...,ab}.

So, every x {1,2,...,77} that is relatively prime to 77 satisfies a pair of congruences of the form

x ≡ n1 (mod 7)

and

x ≡ n2 (mod 11)

where n1 is relatively prime to 7 and n2 is relatively prime to 11. Conversely, for every pair (n1,n2) with (a) 1 ≤ n1 ≤ 7, (b) 1 ≤ n2 ≤ 11, (c) gcd(n1,7) = 1, and

(d) gcd(n2,11) = 1, there is a unique x between 1 and 77 that satisfies the above congruences and is relatively prime to 77.

Therefore ϕ(77) equals the number of choices for (n1,n2), and that is precisely ϕ(7) ×ϕ(11) = 6 ×10 = 60.

Arrays

69

A careful reading of the analysis of ϕ(77) reveals that the only facts we used about 7 and 11 is that they are relatively prime. Thus the argument can be rewritten to give a proof of the following.

Proposition 5.3. If a and b are relatively prime, then ϕ(ab) = ϕ(a)ϕ(b).

From these propositions we derive the following formula for ϕ(n).

Theorem 5.4. Let n be a positive integer and let p1, p2,..., pt be the distinct prime divisors of n. Then

ϕ(n) = n	1 − p1	1 − p2			···	1 − pt	.
	1			1		1
Proof. We factor n into primes as
		e1	e2		et
	n = p1		p2	··· pt

where the pj are distinct primes and the ej are positive integers. By Propositions 5.3 and 5.1 we have

ϕ(n) = ϕ(pe11 )ϕ(pe22 ) ···ϕ(ptet )

= pe11−1(p1 −1) pe21−1(p2 −1) ··· ptet −1(pt −1)

= p1e1		1 − p1		p2e2	1 − p2		··· ptet		1 − pt
			1			1			1
= n	1 − p1		1 − p2			···	1 − pt		.
		1			1		1

Thus, if we know the prime factors of n, we can calculate ϕ(n). This results steers the discussion for the rest of this chapter. We begin by developing algorithms to factor long integers. We create a factor procedure that produces the full prime factorization of an integer; we use an array to hold the result.

5.2Array fundamentals

The first goal of this chapter is to create a procedure that takes an integer and gives us a list of the prime factors of that integer. For example, if the integer is 120, then the list of prime factors is (2,2,2,3,5).

There are several ways to hold lists in C++. The most primitive method is to use an array. An array is simply a list of values of a given type. To declare an array of, say, long integers, we use a statement such as this:

int vals[10];

70	C++ for Mathematicians

This declares val to be an array of ten1 integers. The elements of the array are indexed starting from zero. That is, the ten elements of val are these:

val[0] val[1] val[2] val[3] val[4] val[5] val[6] val[7] val[8] val[9]

Each of these behaves exactly as does a long variable. One can have statements such as val[3]=2*val[0]+val[5]; or ++val[7];. Elements of an array are accessed using square brackets, not parentheses. If an array is declared to have n elements, the first is always indexed by subscript 0 and the last by subscript n − 1. Element n of such an array is not defined.

To print out all the members of an array, one can use a for loop:

for (int k=0; k<10; k++) { cout << val[k] << endl;

}

However, the following does not work: cout << val << endl;. This does not print the values of the val array, even though it is legal C++. Furthermore, arithmetic operations on elements of an array must be specified element by element. If you wish to increase every element of an array by one, you need a statement such as this: for(int k=0; k<10; k++) val[k]++;. Unfortunately, the statement val++; is not illegal, but it does not do what you might want or expect. If alpha and beta are two arrays, the expression alpha==beta does not determine if corresponding elements in the arrays are equal, but the compiler might not complain because it is a legal C++ expression.

Worse yet, if val is an array of ten elements, then only val[0] through val[9] are legitimate array elements. Surprisingly, using val[10] or val[17] is not prevented, but severe problems may result from accessing elements beyond the normal bounds of an array.

The C++ array is the most computationally efficient mechanism for storing a list, but it is not the only manner. We explore several alternatives in Chapter 8.

It is helpful to have a basic understanding the underlying mechanism by which arrays work. To do that, we need to mention the concept of a pointer. In this book we scrupulously avoid dealing with pointers; they are confusing and easily lead to subtle programming errors. However, we do need to be at least vaguely aware that C++ uses pointers. What is a pointer? A pointer is a variable whose value is a location in the computer’s memory. Each byte of a computer’s memory is numbered, and a pointer holds such a number. In the case of arrays, the name of the array is actually a pointer to the location of the first element of the array (e.g., the memory location that holds val[0]). So a statement of the form val++; changes the pointer val so that it now points to a different location in memory. (In fact, it would now point to the location of val[1] but this is more than we want or need to know right now.) A statement of the form cout << val << endl; prints out the location

1Of course, we can declare an array to hold hundreds or thousands of elements. However, some compilers place a limit on the maximum size of an array. There is a simple way to exceed that limit and this is explained in the footnote accompanying Program 5.8 on page 82.

Arrays

71

number where the first element (with subscript zero) of val is housed. Try it! You should get a result that looks something like this: 0xbffffa20. The leading 0x is C++’s way of indicating that the number that follows is in base-16. The remainder is the number in standard base-16 (where a is a digit equal to ten, b is a digit equal to eleven, and so on).

When you declare an array with a statement such as int val[10]; the computer sets aside a block of memory to hold ten ints. When you refer to, say, val[3], the computer figures out where in memory this quantity sits by adding 3*sizeof(int) to the base address val.

You may be thinking: Why do I need to know all this stuff?!? Mostly, you do not need to worry about this. The important things you need to know are these.

•Array elements are indexed starting from 0.

•Operations cannot be performed on arrays as a whole; you need to operate on the elements of arrays individually.

•The name of the array is a valid C++ entity. Once declared, the only information that the name of the array carries is the location in memory of the start of the array and the type of elements held in the array.

The third point is important when we write procedures that have array arguments. In the next section we develop a procedure called factor that takes two arguments. The first is the integer we want to factor. The second is an array to hold the prime factors.

5.3A procedure to factor integers

In this section we develop a procedure to factor integers. The input to the procedure is a long integer, n. The procedure gives us two results: an array holding the prime factors and an integer telling us how many prime factors (multiplicities counted) we found. For example, if the number to be factored is 20, the procedure finds that the factors are (2,2,5) and that there are three prime factors.

We use C++’s return statement to report on the number of factors found. However, the return mechanism does not work for arrays. Instead, we give the procedure an array in the argument list.

The declaration of the factor procedure is this:

long factor(long n, long* flist);

The first argument, n, is the number we want to factor. The second argument, flist, is an array to hold the answer (the factors of n). The type of flist is long*. The star indicates that flist is an array. (Technically, the star signifies that flist is a pointer.) The long to the left of the word factor signifies that factor returns a long integer—the number of prime factors found.

72	C++ for Mathematicians

The technicality that flist is a pointer is relevant to us only in this regard. Although C++ uses call by value as its standard way to pass arguments to procedures, the array flist is not copied to the factor procedure. Instead, the location (i.e., pointer) of the array is sent. Therefore the factor procedure is capable of modifying the elements of flist in the procedure that called flist.

The factor procedure has no mechanism to ensure that the array flist contains enough elements to hold the primes found. It is the responsibility of the user to be sure that the array is large enough to hold the answers. If the size of your computer’s long integers is 4 bytes, then these integers are less than 232. Therefore, no long can have more than 32 factors. Perhaps you have a computer in which long integers are 64 bits; in that case, if flist has at least 64 elements, no problem can arise. This means the procedure that calls factor should look something like this:

long prime_factors[200]; long nfactors;

nfactors = factor(60, prime_factors);

After this code runs nfactors holds 4 and the array prime_factors holds the values 2, 2, 3, and 5 in positions 0 through 3, respectively. Positions 4 through 199 are unaffected by this call to factor.

It is time to design the factoring algorithm. We begin by handling exceptional cases. What if the user sends a negative number or zero to be factored? For negative values, we can simply replace the argument by its absolute value. Asking to factor zero is asking for trouble. Returning a value of 0 is not quite right, because that is what we would return when we are asked to factor 1. Returning a positive value is also misleading. We settle for the unhappy choice of returning the value −1 in case the user requests a prime factorization of 0.

If we are asked to factor 1, we simply report there are no prime factors (i.e., return 0;).

So we suppose that n ≥ 2. The method we use to factor n is this. We check if n is divisible by 2. If so, we record 2 in the first element (index zero) of flist and replace n by n/2. We keep track of where we last wrote into the flist array with a variable we name idx. Initially idx equals zero. Once we record a prime in flist[idx] we then increase idx by 1. We keep doing this until n is no longer divisible by 2.

At that point, we check n for divisibility by 3. We continue dividing until no factors of 3 remain. With each factor of 3 that we find, we record a 3 in flist[idx], advance idx by one, and replace n by n/3.

It would be logical to next try divisibility by 5, but C++ does not “know” that 5 is the next prime. Instead, we try divisibility by 4 which, of course, fails because we have already divided out all of n’s factors of 2.

We continue in this manner, trying a divisor until it is exhausted and then advancing to the next divisor. When do we stop? Once we have divided out all of n’s prime factors, it equals 1; that’s when we stop.

The algorithm is depicted in the flowchart in Figure 5.1.

Arrays

73

# $

%

!

"

' (

&

#

$

%

*

+

)

"

!

)

"

,

'

(

)

"

Figure 5.1: A flowchart for the factoring algorithm.

We are now ready to implement this procedure in C++. To begin, we write the header file, factor.h.

Program 5.1: Header file factor.h for the first version of the factor procedure.

1 #ifndef FACTOR_H

2#define FACTOR_H

3

4/**

5* Factor an integer n. The prime factors are saved in the second

6* argument, flist. It is the user’s responsibility to be sure that

7* flist is large enough to hold all the primes. If n is negative, we

8* factor -n instead. If n is zero, we return -1. The case n equal to

9	*	1 causes this procedure to return 0 and no primes are saved in
10	*	flist.

74	C++ for Mathematicians

11*

12* @param n the integer we wish to factor

13* @param flist an array to hold the prime factors

14* @return the number of prime factors

15*/

16long factor(long n, long* flist);

17

18 #endif

Next is the C++ source code.

Program 5.2: Source file factor.cc for the first version of the factor procedure.

1#include "factor.h"

2

3long factor(long n, long* flist) {

4

5 // If n is zero, we return -1

6if (n==0) return -1;

7

8 // If n is negative, we change it to |n|

9if (n<0) n = -n;

10

11// If n is one, we simply return 0

12if (n==1) return 0;

13

14 // At this point we know n>1

15

16int idx = 0; // index into the flist array

17int d = 2; // current divisor

18

19while (n>1) {

20while (n%d == 0) {

21flist[idx] = d;

22++idx;

23n /= d;

24}

25++d;

26}

27return idx;

28}

Lines 5–12 deal with the exceptional cases.

Line 16 sets up the variable idx. Throughout the procedure idx contains the index of the next location in the flist array where we record the prime factors. We set idx equal to zero so the first prime we record goes into the first element of flist.

Line 17 sets up the variable d. This is the current divisor we are testing.

The heart of the procedure lies in lines 19–26. We keep trying to find factors of n as long as there are factors to be found (i.e., as long as n holds a value greater than 1). On line 20 we test n for divisibility by d. If successful, we record d in the array flist at location idx, increase idx by one so it refers to the next available cell in flist, and divide out the factor of d from n. We keep doing this until n is no longer divisible by d. At that point, we increase d by one and continue.

Arrays

75

At the end, idx has been increased once for every prime factor we found. So we simply return its value at line 27.

Here is a main to test the factor procedure. This program prints out the prime factorization of all integers from 1 to 100.

Program 5.3: A main to test the factor procedure.

1 #include "factor.h"

2 #include <iostream>

3using namespace std;

4

5/**

6* A program to test the factor procedure.

7	*/
8

9int main() {

10

11 long flist[100]; // place to hold the factors

12

13for (long n=1; n<=100; n++) {

14int nfactors = factor(n,flist);

15cout << n << "\t";

16for (int k=0; k<nfactors; k++) cout << flist[k] << " ";

17cout << endl;

18}

19}

On line 11 we declare flist to hold 100 long integers. The for loop on lines
13–18 requests that we factor n for all values from 1 to 100. The call to factor is
on line 14. We save the number of factors found in a variable named nfactors and
the factors themselves populate the elements of the array flist.
Line 15 prints the current value of n and then a tab character. The sequence \t
stands for a tab; this way the list of factors ends up nicely arranged.
Line 16 prints out the factors of n separated by spaces and then line 17 starts a
new line on the screen.
The output looks like this.

1

22

33

42 2

55

62 3

77

82 2 2

93 3

102 5

1111

(many lines deleted)

955 19

962 2 2 2 2 3

9797