# MTHS 3600 - HW22

Date: __4/14/14______________

The following function, build_F, returns a two-dimensional list, F, with N rows and M columns.  The entry at F[i][j] is the tuple (i,j+1). Evaluate the following block to define the function and to use it to build a 10-by-10 table. Note that the entries in this table can be interpreted as all possible non-negative fractions whose numerator, $n$, satisfies $0 \leq n < N$ and whose denominator, $d$, satisfies  $0 < d \leq M$.

def build_F(N,M): F = [] for i in range(N): row = [] for j in range(M): row.append((i,j+1)) F.append(row) return F N = 10 M = 10 F = build_F(N,M) for i in range(N): print F[i]
 [(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10)] [(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9), (1, 10)] [(2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10)] [(3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10)] [(4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (4, 6), (4, 7), (4, 8), (4, 9), (4, 10)] [(5, 1), (5, 2), (5, 3), (5, 4), (5, 5), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10)] [(6, 1), (6, 2), (6, 3), (6, 4), (6, 5), (6, 6), (6, 7), (6, 8), (6, 9), (6, 10)] [(7, 1), (7, 2), (7, 3), (7, 4), (7, 5), (7, 6), (7, 7), (7, 8), (7, 9), (7, 10)] [(8, 1), (8, 2), (8, 3), (8, 4), (8, 5), (8, 6), (8, 7), (8, 8), (8, 9), (8, 10)] [(9, 1), (9, 2), (9, 3), (9, 4), (9, 5), (9, 6), (9, 7), (9, 8), (9, 9), (9, 10)] [(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10)] [(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9), (1, 10)] [(2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10)] [(3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10)] [(4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (4, 6), (4, 7), (4, 8), (4, 9), (4, 10)] [(5, 1), (5, 2), (5, 3), (5, 4), (5, 5), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10)] [(6, 1), (6, 2), (6, 3), (6, 4), (6, 5), (6, 6), (6, 7), (6, 8), (6, 9), (6, 10)] [(7, 1), (7, 2), (7, 3), (7, 4), (7, 5), (7, 6), (7, 7), (7, 8), (7, 9), (7, 10)] [(8, 1), (8, 2), (8, 3), (8, 4), (8, 5), (8, 6), (8, 7), (8, 8), (8, 9), (8, 10)] [(9, 1), (9, 2), (9, 3), (9, 4), (9, 5), (9, 6), (9, 7), (9, 8), (9, 9), (9, 10)]

Using upward sloping diagonals we can write the infinite version of this two dimensional array as an infinite one dimensional list.  It would start as follows:

[(0,1), (1,1), (0,2), (2,1), (1,2), (0,3), (3,1), . . . ]

In the following block the function, build_F1d, returns the initial portion of this one dimensional list.  Set N to 7 to see the same initial list.  Then set N to 55 to see the list the "fractions" in the upper half of the preceding array.

def build_F1d(n): F1d = [] if n <= 0: return F1d F1d = [(0,1)] count = 1 for diag in range(1,n): i = diag j = 1 while 0 <= i: F1d.append((i,j)) count += 1 if count == n: return F1d i -= 1 j += 1 return F1d N = 7 F1 = build_F1d(N) print F1
 [(0, 1), (1, 1), (0, 2), (2, 1), (1, 2), (0, 3), (3, 1)] [(0, 1), (1, 1), (0, 2), (2, 1), (1, 2), (0, 3), (3, 1)]

The previous two blocks of code define a procedure for uniquely identifying every non-negative fraction with a positive integer and vice versa.  This shows that there is a one-to-one corresponce between the non-negative fractions and the positive integers.  In other words, the non-negative fractions are countably infinite.

To make this a formal bijection we really need two functions, NonNegFracToPosInt, and its inverse, PosIntToNonNegFrac. Given a non-negative fraction NonNegFracToPosInt returns the corresponding positive integer.  Given a positive interger PosIntToNonNegFrac returns the corresponding non-negative fraction.

Exercise 1:  Construct the functions, NonNegFracToPosInt and PosIntToNonNegFrac, which together define the bijection between the positive integers and the non-negative fractions.  Complete the definitions of the two functions and avoid using any loops.

Hint, the Triangular Numbers have come up several times during this semester.

def NonNegFracToPosInt(i,j): if i<0 or j<1: return None diag = i+j n = ((diag-1)*(diag))/2 + j return n def PosIntToNonNegFrac(n): if n <= 0: return None # # # return (i,j)
print NonNegFracToPosInt(0,1), PosIntToNonNegFrac(1) print NonNegFracToPosInt(0,3), PosIntToNonNegFrac(6) print NonNegFracToPosInt(2,2), PosIntToNonNegFrac(8) print NonNegFracToPosInt(9,1), PosIntToNonNegFrac(46) print NonNegFracToPosInt(0,10), PosIntToNonNegFrac(55)
 1 Traceback (click to the left of this block for traceback) ... NameError: global name 'j' is not defined 1 Traceback (most recent call last): File "", line 1, in File "_sage_input_6.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cHJpbnQgTm9uTmVnRnJhY1RvUG9zSW50KDAsMSksIFBvc0ludFRvTm9uTmVnRnJhYygxKQpwcmludCBOb25OZWdGcmFjVG9Qb3NJbnQoMCwzKSwgUG9zSW50VG9Ob25OZWdGcmFjKDYpCnByaW50IE5vbk5lZ0ZyYWNUb1Bvc0ludCgyLDIpLCBQb3NJbnRUb05vbk5lZ0ZyYWMoOCkKcHJpbnQgTm9uTmVnRnJhY1RvUG9zSW50KDksMSksIFBvc0ludFRvTm9uTmVnRnJhYyg0NikKcHJpbnQgTm9uTmVnRnJhY1RvUG9zSW50KDAsMTApLCBQb3NJbnRUb05vbk5lZ0ZyYWMoNTUp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpKqB2Nv/___code___.py", line 3, in print NonNegFracToPosInt(_sage_const_0 ,_sage_const_1 ), PosIntToNonNegFrac(_sage_const_1 ) File "/tmp/tmp42rPHX/___code___.py", line 16, in PosIntToNonNegFrac return (i,j) NameError: global name 'j' is not defined

We say that two fractions $a/b$ and $c/d$ are equivalent if and only if $a\cdot d = c\cdot b$. It is straightforward to verify that this requirement defines an equivalence relation on the positive "fractions" as described above, namely all pairs of positive integers where the second integer is always nonzero.

Exercise 2.  On a separate sheet of paper verify that this proposed relationship is an equivalence relationship because it satisfes the following three defining properties.

1. reflexive:  $a/b \sim a/b$.
${\raise0.7ex\hbox{a} \!\mathord{\left/ {\vphantom {a b}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{b}} \sim {\raise0.7ex\hbox{c} \!\mathord{\left/ {\vphantom {c d}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{d}}$
2. symmetric:  $a/b \sim c/d$  implies that  $c/d \sim a/b$.
3. transitive: if  $a/b \sim c/d$  and  $c/d \sim p/q$  then  $a/b \sim p/q$.

Exercise 3.  On the printout from the first block, use pencil(s), pen(s), or highlighter(s), to circle or highlight each set of equivalent fractions containing more than 3 fractions. Hint, the first row is one such set, but many sets have less than 4 members (some only 1) and should not be hightlighted. Also write a short program which prints the same table, but prints the pair of integers in the conventional fraction form $p/q$. This version of the table should quickly confirm your hand drawn results.

def build_F_rational(N,M): F = [] for i in range(N): fraction = [] for j in range(M): fraction.append(i/(j+1)) F.append(fraction) return F N = 10 M = 10 F = build_F_rational(N,M) for i in range(N): for j in range(M): print F[i][j],'\t', print
 0 0 0 0 0 0 0 0 0 0 1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 2 1 2/3 1/2 2/5 1/3 2/7 1/4 2/9 1/5 3 3/2 1 3/4 3/5 1/2 3/7 3/8 1/3 3/10 4 2 4/3 1 4/5 2/3 4/7 1/2 4/9 2/5 5 5/2 5/3 5/4 1 5/6 5/7 5/8 5/9 1/2 6 3 2 3/2 6/5 1 6/7 3/4 2/3 3/5 7 7/2 7/3 7/4 7/5 7/6 1 7/8 7/9 7/10 8 4 8/3 2 8/5 4/3 8/7 1 8/9 4/5 9 9/2 3 9/4 9/5 3/2 9/7 9/8 1 9/10  0 0 0 0 0 0 0 0 0 0 1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 2 1 2/3 1/2 2/5 1/3 2/7 1/4 2/9 1/5 3 3/2 1 3/4 3/5 1/2 3/7 3/8 1/3 3/10 4 2 4/3 1 4/5 2/3 4/7 1/2 4/9 2/5 5 5/2 5/3 5/4 1 5/6 5/7 5/8 5/9 1/2 6 3 2 3/2 6/5 1 6/7 3/4 2/3 3/5 7 7/2 7/3 7/4 7/5 7/6 1 7/8 7/9 7/10 8 4 8/3 2 8/5 4/3 8/7 1 8/9 4/5 9 9/2 3 9/4 9/5 3/2 9/7 9/8 1 9/10 

When a pair of integers is explicitly written in the conventional fraction form, $p/q$, Sage automatically "reduces the fraction to lowest terms," and this fraction becomes the conventional representative for the set of equivalent fractions.  However the other equivalent fractions are absolutely essential for the arithmetic operations of addition, subtraction, multiplication, and division. So the more sophisticated mathematical view of the set of Rational Numbers, Q, is that Q is the set of equivalence classes of fractions, and that each equivalence class is conventionally represented by the fraction in the set which has been reduced to lowest terms. However, any fraction in the set can be used for arithmetic calculations. Given any member of a class, say 1234567890987654321 / 123451234598769876, how can we find the conventional representative? The answer is Euclid's algorithm.

## Euclid's Algorithm

Euclid's algorithm is an efficient method for computing the greatest common divisor (GCD) of two integers. It is named after the Greek mathematician Euclid, who described it in Books VII and X of his Elements.  The algorithm has many theoretical and practical applications. It is a key element of the RSA algorithm. It is used to solve Diophantine equations, such as finding numbers that satisfy multiple congruences (Chinese remainder theorem) or multiplicative inverses of a finite field. It can also be used in the construction of continued fractions, in the Sturm chain method for finding real roots of a polynomial, and in several modern integer factorization algorithms. Finally, it is a basic tool for proving theorems in modern number theory, such as Lagrange's four-square theorem and the fundamental theorem of arithmetic (unique factorization).

The greatest common divisor g is the largest natural number that divides both a and b without leaving a remainder. In Sage the greatest common divisor is given by the function, gcd(a,b).

Consider the problem of determining the gcd(252,105). If we reduce the larger number by subtracting the smaller one from it, the gcd does not change (since it divides both terms it must divide the difference):

gcd(252,105) = gcd(252-105,105) = gcd(147,105).

Subtract again:

gcd(147,105) = gcd(147-105,105) = gcd(42,105).

Now 105 is no longer the smaller number. However, continuing in the same way, we reduce the larger number, now 105, by subtracting the smaller one from it, the gcd is unchanged:

gcd(42,105) = gcd(42,105-42) = gcd(42,63).

The first number, 42, is still the smaller one, so again we use it to reduce the larger one:

gcd(42,63) = gcd(42,63-42) = gcd(42,21),

and finally:

gcd(42,21) = gcd(42-21,21) = gcd(21,21) = 21.

Of course, we can speed up the algorithm by dividing and using the remainder. Then the algorithm has the following steps.

252 = 2*105 + 42

105 = 2*42 + 21

42 = 2*21

Exercise 4. Complete the following function, which implements Euclid's Algorithm using the integer divide. Then test it by evaluating the subsequent block.

def euclid(a,b): if a < b: dividend = b divisor = a else: dividend = a divisor = b quotient = dividend//divisor remainder = dividend - quotient*divisor while 0 < remainder: dividend = divisor divisor = remainder quotient = dividend//divisor remainder = dividend - quotient*divisor return divisor
print euclid(13,6) print euclid(6,14) print euclid(18,6) print euclid(105,252) print euclid(1989,867)
 1 2 6 21 51 1 2 6 21 51

## The Extended Euclidean Algorithm

The following function is a very straightforward translation of the recursive pseudo-code function presented in the Wikipedia article on the Extended Euclidean Algorithm. This algorithm returns $x$ and $y$, the smallest (in absolute value) integer coefficients that satisfy Bezout's identity,

$x*a + y*b = r$

where $a$ and $b$ are given integers and $r$ is their gcd.

# The Extended Euclidean Algorithm # a and b are integers with b < a. def extended_euclid(a, b): if b == 0: return (1, 0) else: q = a//b r = a - q*b s, t = extended_euclid(b, r) return (t, s-q*t)
a,b = 252,105 x,y = extended_euclid(a,b) print x,'*',a,'+',y,'*',b,' = ',(x*a + y*b) print euclid(a,b)
 -2 * 252 + 5 * 105 = 21 21 -2 * 252 + 5 * 105 = 21 21
a,b = 1989,867 x,y = extended_euclid(a,b) print x,'*',a,'+',y,'*',b,' = ',(x*a + y*b) print euclid(a,b)
 7 * 1989 + -16 * 867 = 51 51 7 * 1989 + -16 * 867 = 51 51

If the gcd of $a$ and $b$ is 1, then we say that $a$ and $b$ are coprime (or relatively prime). When $a$ and $b$ are coprime, then Bezout's identity provides some useful information. The identity becomes

$x*a+y*b = 1$,

which means that $x$ is the multiplicative inverse of $a$ mod $b$ and that $y$ is the multiplicative inverse of $b$ mod $a$. This is illustrated in the following block.

a,b = 120,23 x,y = extended_euclid(a,b) print x,'*',a,'+',y,'*',b,' = ',(x*a + y*b) print (x*a),'mod',b,'=',(x*a)%b print (y*b),'mod',a,'=',(y*b)%a
 -9 * 120 + 47 * 23 = 1 -1080 mod 23 = 1 1081 mod 120 = 1 -9 * 120 + 47 * 23 = 1 -1080 mod 23 = 1 1081 mod 120 = 1

## The Calkin-Wilf Tree

The Calkin-Wilf Tree was introduced in a 2000 paper by Neil Calkin and Herbert Wilf. It is a binary tree with each node a fraction of the form $a/b$.  The root node is $1/1$. For every node, $a/b$, the left child is $a/(a+b)$ and the right child is $(a+b)/b$.  Using the same indexing functions that were introduced with the priority queue we can easily store the tree in a one-dimensional list.

def parent(node): return (node-1)//2 def leftchild(node): return 2*node + 1 def rightchild(node): return 2*node + 2

Exercise 5.  Complete the function, Calkin_Wilf, that given $N$, returns a list of the first $N$ nodes Calkin-Wilf tree (viewed as a complete binary tree, filled in from top-to-bottom, left-to-right).

def Calkin_Wilf(N): if N<1: return None tree =  for p in range(N): a = numerator(tree[p]) b = denominator(tree[p]) # # return tree
print Calkin_Wilf(7) print Calkin_Wilf(15)
 Traceback (click to the left of this block for traceback) ... IndexError: list index out of range Traceback (most recent call last): File "", line 1, in File "_sage_input_17.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cHJpbnQgQ2Fsa2luX1dpbGYoNykKcHJpbnQgQ2Fsa2luX1dpbGYoMTUp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpqcw5zC/___code___.py", line 3, in print Calkin_Wilf(_sage_const_7 ) File "/tmp/tmp7aptmR/___code___.py", line 8, in Calkin_Wilf a = numerator(tree[p]) IndexError: list index out of range

The following function takes the list from the Calkin_Wilf function and displays it as a tree.  Examine the tree carefully to make sure that it is correct.

def CW_tree_display(N): CW = Calkin_Wilf(N) # print CW row = 1 right_index = 2^row - 1 k = 0 while k < N: while k < right_index: print CW[k],'\t', k += 1 if N <= k: break print k = right_index row += 1 right_index = 2^row - 1 CW_tree_display(31)
 Traceback (click to the left of this block for traceback) ... IndexError: list index out of range Traceback (most recent call last): while k < N: File "", line 1, in File "/tmp/tmpwqzbKn/___code___.py", line 20, in exec compile(u'CW_tree_display(_sage_const_31 ) File "", line 1, in File "/tmp/tmpwqzbKn/___code___.py", line 4, in CW_tree_display CW = Calkin_Wilf(N) File "/tmp/tmp7aptmR/___code___.py", line 8, in Calkin_Wilf a = numerator(tree[p]) IndexError: list index out of range

Exercise 6.  Examine each row in the tree displayed above and answer the following two questions.

1. Is this a binary search tree, that is, are the nodes in order by magnitude?
2. Is there a pattern relating the denominators and numerators in each row? If so, what is the pattern?

Exercise 7.  In the block below list the first 31 numerators returned by your Calkin_Wilf routine.  Then in the next line list the first 31 denominators. Is there a pattern connecting the list of numerators to the list of denominators? If so, briefly describe it.

Exercise 8.  Complete the code for the recursive function, fusc, which is defined by the following properties:

1.  fusc(0) = 0
2.  fusc(1) = 1
3.  fusc(2n) = fusc(n)
4.  fusc(2n+1) = fusc(n) + fusc(n+1)

Then evaluate the subsequent block.  Hint, the first 6 values are: 0, 1, 1, 2, 1, 3.

Is there a connection between the sequence returned by the fusc function and the Calkin-Wilf tree?  If so, describe it.

def fusc(N): if N < 0: return None if N == 0: y = elif N == 1: y = elif N%2==0: # N is even, let N = 2*n n = y = else: n = y = return y
 Traceback (click to the left of this block for traceback) ... SyntaxError: invalid syntax Traceback (most recent call last): y = File "", line 1, in File "/tmp/tmpnOjMnS/___code___.py", line 7 y = ^ SyntaxError: invalid syntax
[fusc(k) for k in range(64)]
 Traceback (click to the left of this block for traceback) ... NameError: name 'fusc' is not defined Traceback (most recent call last): File "", line 1, in File "_sage_input_22.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("W2Z1c2MoaykgZm9yIGsgaW4gcmFuZ2UoNjQpXQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpIuCsF_/___code___.py", line 3, in exec compile(u'[fusc(k) for k in range(_sage_const_64 )] File "", line 1, in NameError: name 'fusc' is not defined

Exercise 9.  Is every node in the Calkin-Wilf tree a fraction that has been reduced to lowest terms? If so, provide the outline of a proof. If not, search for a counter-example and indicate the extent of your search and the results.

Exercise 10.  Does every positive rational number occur in the Calkin-Wilf tree? If you think so, provide the reason(s) that support this conclusion.

Based on your explorations in this and the previous exercise, what do you conclude about the claim that every rational number occurs exactly once in the Calkin-Wilf tree?