Sue Robinson - Fall 1995

Final Project
Phi Number System

Topic

In the October 1995 issue of Mathematics Magazine there is an article by Cecil Rousseau about the phi number system.  In the article is a proof of a theorem, and within the proof is an algorithm.  For my topic I tried to write that algorithm in Mathematica.

Theorem
For every positive integer n, there is a corresponding finite sequence of distinct integers k1, k2, ... km such that n = phi^k1 + phi^k2 + ... + phi^km, where phi = (1 + Sqrt 5)/2, and k1>k2>k3>...>km.

Algorithm
The algorithm in the proof is actually finding the distinct k's so that the positive integer n can be expressed as a finite series of powers of phi.

Phi Number System

The phi number system has the golden ratio, phi = (1  + Sqrt 5)/2, as its base.  Any positive real number can be expressed to base phi using the digits 0 and 1.
For the base representations to be unique there cannot be two consecutive powers of phi.  If ther are two consecutive powers of phi they can be combined to give a different power, and so the representation is not unique.

Proof that two consecutive powers of phi can be expressed as a different power of phi

phi^2 = ( (1+ Sqrt 5)/2)^2
          = (1 + 2 Sqrt 5 + 5)/4
          = 3/2 + (Sqrt 5)/2
          = 1 + ((1 + Sqrt 5)/2)
phi^2 = phi + 1

phi^3 + phi^4 = phi^3 (1 + phi),     but (1 + phi) = phi^2
                       = phi^3 phi^2
phi^3 + phi^4 = phi^5

Algorithm

The algorithm part of the proof is actually finding the powers of phi (i.e. k1, k2, ... km) for some positive integer n.  The algorithm works in the following way:

start with sigma = n, a positive integer
find k1 such that phi^k1 <= (sigma) < phi^(k1+1)
find k2 such that phi^k2 <= (sigma - phi^k1) < phi^(k2+1)
find k3 such that phi^k3 <= (sigma - phi^k1 - phi^k2) < phi^(k3+1)
|
|
find km such that phi^km <= (sigma - phi^k1 - phi^k2 - ... - phi^(km-1)) < 		phi^(km+1)
until sigma - phi^k1 - phi^k2 - ... - phi^km = 0

Then n = phi^k1 + phi^k2 + phi^k3 + ... + phi^km, where the k's are distinct and k1>k2>k3>...>km

My Approach to the Algorithm

My implementation of the algorithm is to find the first k such that phi^k <= n < phi ^(k+1) by starting at k = 1, and increasing the value of k until I find the phi^k that is closest to n.

Since I know that k1>k2>k3>... >km, I decreased the value of k by 1 until I found k2 such that phi^k2 <= (n - phi^k1) < phi^(k2+1).  I continued to work backwards through values of k until I found km such that n = phi^k1 + phi^k2 +... + phi^km.

Problems

I encountered some problems trying to implement this algorithm becausse Mathematica did not do what I expected it to do in some circumstances.  For example, it did not evaluate expressions of the form 3 + Sqrt[5] < 0, and it did not evaluate expressions with square roots in the denominator.  My first attempt to solve the problem was to use the numeric function N[ ], but accuracy was lost in the approximations, and the program did not terminate properly.  See the following example:

phi = (1 + Sqrt[5])/2;

N[6-phi^3-phi^1-phi^-4] == 0

N[6-phi^3-phi^1-phi^-4]

Clear[phi]

Mathematica did not consider -4.33681 10-19 to be equal to zero, and did not terminate at this point.  In fact, the powers of phi that I got for the positive integer 6 using numeric approximations were {3, 1, -5, -7, -9, -11, -13, -15, -17, -19, -21, -22}.  I realized that this answer was incorrect because two consecutive digits were 1 (i.e phi^-21 and phi^-22), which does not give a unique answer.  I traced the problem back to k = -4 and discovered the problem described above.  Phi^-22 + phi^-21 can be expressed as phi^-20, then phi^-20 + phi^-19 can be expressed as phi^-18 etc. until phi^-6 + phi^-5 = phi^-4.  There are then no two consecutive powers of phi.  And {-4, 1, 3} is a distinct set of powers of phi for the positive integer 6.

I wrote two functions called PhiPower and PositiveValue to overcome the problems of evaluation and correct termination.

PhiPower

When I used Mathematica to evaluate phi^m where m is a negative integer, Mathematica left the expression in a 2^m/(1+Sqrt[5])^m format, see the following example:

phi = (1+Sqrt[5])/2;

phi^(-2)

This format created problems with evaluations because Mathematica did not evaluate the denominator.  It left the denominator in a form containing square roots.  See the following examples:

phi + phi^-1

Expand[%]

Together[%]

I determined the expression for phi^-1 that did not have square roots in the denominator.

phi = (1 + Sqrt 5)/2

phi^-1 = 1/phi = 2/(1 + Sqrt 5)
           = 2 (1 - Sqrt 5)/(1 + Sqrt 5)(1 - Sqrt 5)
           = 2 (1- Sqrt 5)/ (1-5)
phi^-1 = (Sqrt 5 - 1)/2

In the following example the name phinegative refers to phi^-1.

phinegative = (Sqrt[5]-1)/2;

phi + phinegative

Expand[%]

Clear[phi]

I then wrote a function called PhiPower which would handle the negative powers of phi differently to the positive powers.

PhiPower[n_] := Which[n>0, ((Sqrt[5]+1)/2)^n,
					  n==0, 1,
					  n<0, ((Sqrt[5]-1)/2)^-n]

Examples using PhiPower

PhiPower[2]

PhiPower[0]

PhiPower[-4]

PositiveValue

I still had the problem that Mathematica does not evaluate expressions like
6 + 5/2 -17 (Sqrt 5) /2 > 0, and I needed to check these expressions when finding the powers of phi.  The following are some examples of Mathematica's evaluation of > tests.

Examples:

6>0

-2>0

6+5/2-17 Sqrt[5]/2>0

I used pattern matching to separate the parts of an expression of the form c + d Sqrt 5 to find c and d.  I then used these values of c and d to determine whether c + d Sqrt 5 > 0.  If c + d Sqrt 5 is positive then I wanted my function to return the value "True".  There are many different cases to consider when determining if the value is positive:

1)  For c > 0 and d > 0, it is obvious that c + d Sqrt 5 > 0, so expression is True.
 
2)  For c < 0 and d < 0, it is obvious that c + d Sqrt 5 < 0, so expression is False.

3)  For c > 0 and d < 0,  for c + d Sqrt 5 > 0
                                                   d Sqrt 5 > -c
                                                     Sqrt 5 < -c/d     since d<0
                                                            5 < (-c)^2/d^2
                                                         c^2 > 5 d^2      for c + d Sqrt 5 > 0, so expression is True.
                                                     
 4)  For c < 0 and d > 0,  for c + d Sqrt 5 > 0
                                                    d Sqrt 5 > -c
                                                       Sqrt 5 > -c/d     since d>0
                                                              5 > c^2/d^2
                                                          c^2 < 5 d^2      for c + d Sqrt 5 > 0, so expression is True.

PositiveValue[c_+d_ Sqrt[5]] := 
		Which[c>0 && d>0, True,
			  c<0 && d<0, False,
			  c>0 && d<0, c^2>5*d^2,
			  c<0 && d>0, c^2<5*d^2]

The expression 6+5/2-17 Sqrt[5]/2 > 0 can now be checked using the PositiveValue function.

PositiveValue[6+5/2-17 Sqrt[5]/2]

However, this PositiveValue function only works when the number given to it is of the form c + d Sqrt 5.  If c or d is zero and therefore omitted, the pattern does not match.  Further definitions of PositiveValue are needed to deal with these cases, and they are given below.

PositiveValue[c_/; IntegerQ[c]] :=
		Which[c>0, True,
			  c<0, False]

PositiveValue[d_ Sqrt[5]] :=
		Which[d>0, True,
			  d<0, False]

PositiveValue[0] := True

Examples using PositiveValue:

PositiveValue[3+6 Sqrt[5]]

PositiveValue[-1+3 Sqrt[5]]

PositiveValue[-1-7 Sqrt[5]]

PositiveValue[4-2 Sqrt[5]]

PositiveValue[4-1 Sqrt[5]]

PositiveValue[4]

PositiveValue[-4]

PositiveValue[0]

PositiveValue[-1 Sqrt[5]]

First Implementation of Algorithm

My first implementation of the algorithm just found the powers of phi for a positive integer, and printed those powers as shown below.  The last line of the program checks that the series does add up to the initial positive integer; it is not necessary but shows that the answer is correct.

Clear[phi]

FindSequence[Sigma_]:=
	(
	p=1;
	m=1;
	While[PositiveValue[Expand[Sigma-PhiPower[p]]], 
				p=p+1];
	k[m]=p-1;
	NewSigma=Together[Expand[Sigma-Sum[
					PhiPower[k[i]],{i,1,m}]]];
	
	While[NewSigma =!= 0,
		(p=k[m]-1;
		While[!PositiveValue[S=Expand[NewSigma-
				PhiPower[p]]], p=p-1];
		m=m+1;
		k[m]=p;
		NewSigma=Together[S];
		)]
	Print["Positive Integer = ", Sigma]
	Print["Powers of phi = ", Table[k[i], {i, 1, m}]]
	Print[Together[Expand[Sum[PhiPower[k[i]],
			{i,1,m}]]], " = ", Sum[phi^k[i], {i,1,m}]];
	)

Examples:

FindSequence[3]

FindSequence[6]

FindSequence[33]

In the above example one of the k values was 0 and Mathematica writes x^0 as 1.  I tried to change the way that Power works, but then none of my program worked, so I decided to leave phi^0 as 1.

FindSequence[1234567890]

Enhancements to the Program

The proof of the theorem in the magazine article used the norms of Z(phi), where 
Z(phi) = {a + b phi} where a and b are integers.  It also showed that the norms are all positive integers and that:
|N(n)| > |N(n - phi^k1)| > |N(n - phi^k1 - phi^k2)| > ... >= 0.

Therefore, I decided to modify my program to produce the norms.  To produce the norms of Z(phi) I wrote two programs; ZofPhi and NormofZ.

ZofPhi

The function ZofPhi takes a number in the form c + d Sqrt 5, and returns an expression of the form a + b phi, where a and b are integers.  Once again several different definitions of ZofPhi are required for correct pattern matching whatever the values of c and d.

Given that c + d Sqrt 5 is to be represented as a + b phi,

a + b phi = c + d Sqrt 5
a + b (1/2 + (Sqrt 5)/2) = c + d Sqrt 5
a + b/2 +(b/2) Sqrt 5 = c + d Sqrt 5
b/2 = d, so b = 2 d
a + b/2 = c, so a = c - b/2, so a = c-d.

Hence the following definitions for ZofPhi:

ZofPhi[c_ + d_ Sqrt[5]]:= (c-d) + 2 d phi

ZofPhi[d_ Sqrt[5]] := 2 d phi

ZofPhi[c_/; IntegerQ[c]] := c

ZofPhi[0] := 0

Examples of XofPhi:

Clear[phi]

ZofPhi[2 + 4 Sqrt[5]]

ZofPhi[3]

ZofPhi[-6 Sqrt[5]]

ZofPhi[1/2 + 3/2 Sqrt[5]]

Notice that all the coefficients of ZofPhi are integers, as expected.

NormofZ

With Z(phi) = a + b phi, where a and b are integers, the norm of Z(phi) = a^2 + a b - b^2.  The following definitions for NormofZ allow for the various values of a and b, and calculate the norm of Z(phi).

NormofZ[a_ + b_ phi] := a^2 + a b - b^2

NormofZ[a_/; IntegerQ[a]] := a^2

NormofZ[b_ phi/; IntegerQ[b]] := - b^2

NormofZ[0] := 0

Examples of NormofZ:

NormofZ[2 + 3 phi]

NormofZ[-2 + 3 phi]

NormofZ[2 - 3 phi]

NormofZ[-2 -3 phi]

NormofZ[2]

NormofZ[-2]

NormofZ[6 phi]

NormofZ[-6 phi]

NormofZ[0]

Program that Gives Norms and K's

The following program is an enhancement of the original program and it produces the norms of Z(phi) as well as the values of k.

Clear[phi]

FindSequence[Sigma_]:=
	(
	p=1;
	m=1;
	While[PositiveValue[Expand[Sigma-PhiPower[p]]], 
				p=p+1];
	k[m]=p-1;
	NewSigma=Together[Expand[Sigma-Sum[
					PhiPower[k[i]],{i,1,m}]]];
	Norm[m]=NormofZ[ZofPhi[Expand[NewSigma]]];
	While[NewSigma =!= 0,
		(p=k[m]-1;
		While[!PositiveValue[S=Expand[NewSigma-
				PhiPower[p]]], p=p-1];
		m=m+1;
		k[m]=p;
		NewSigma=Together[S];
		Norm[m]=NormofZ[ZofPhi[Expand[NewSigma]]];
		)]
	Print["Positive Integer = ", Sigma]
	Print["Powers of phi = ", Table[k[i], {i, 1, m}]];
	Print["Norms of Z(phi) = ", Table[
								Norm[i], {i, 1, m}]];
	Print[Together[Expand[Sum[PhiPower[k[i]],
		{i,1,m}]]], " = ", Sum[phi^k[i], {i,1,m}]]
	)

Examples:

FindSequence[3]

FindSequence[33]

Notice that the norms are decreasing in value as expected.

A Recursion Solution

I also attempted a recursion solution to find the distinct k's, based on my solution to the algorithm.  First of all I wrote functions to find the first k, and a subsequent k.

FindFirstPower

This function takes a number and finds the first k such that phi^k < number < phi^(k+1).

FindFirstPower[number_] :=
	( p=1;
	While[PositiveValue[Expand[number-PhiPower[p]]],
		p+=1];
	Return[p-1]; 
	)

Examples using FindFirstPower

FindFirstPower[3]

FindFirstPower[6]

FindNextPower

This function takes a number and a power and returns the next phi power.

FindNextPower[number_,startpower_] := 
	 (p=startpower-1;
	  While[!PositiveValue[Expand
			[number-PhiPower[p]]],
			p-=1];
	  Return[p];
	 )

Examples using FindNextPower

FindNextPower[Expand[3-PhiPower[2]],2]

FindNextPower[Expand[6-PhiPower[3]],3]

RemainderList

RemainderList is a recursive function that builds a list of the powers of phi.

RemainderList[number_, k_] := 
	Module[{newk, newnumber},
	If[number===0, Return[{}],
		(newk=FindNextPower[number, k];
		 newnumber=Expand[number-PhiPower[newk]];
		 Return[Prepend[RemainderList[newnumber,newk],newk]];
		 )
		 ]
	]

Examples using RemainderList

RemainderList[Expand[3-PhiPower[2]],2]

RemainderList[Expand[6-PhiPower[3]],3]

RemainderList[Expand[33-PhiPower[7]],7]

MakePowerList

This function brings together the previous functions. It takes a number, finds the first k and the value of number-phi^k.  It then uses the recursive function RemainderList to calculate the remaining powers of phi.

MakePowerList[number_] :=
	( k = FindFirstPower[number];
	  newnumber=Expand[number-PhiPower[k]];
	  result = Prepend[RemainderList[newnumber,k],k];
	  Print["Phi powers of ",number, " = ",result];
	)

Examples:

MakePowerList[3]

MakePowerList[6]

MakePowerList[33]

MakePowerList[1234567890]

Enhancements to the Recursive Solution

I decided to include the calculations of the norms in this solution.  To do this I only had to change one of the functions – RemainderList.  The following function is the enhanced RemainderList.

RemainderList[number_, k_] := 
	Module[{newk, newnumber},
	Print["Norm = ",NormofZ[ZofPhi[Expand[number]]]];
	If[number===0, Return[{}],
		(newk=FindNextPower[number, k];
		 newnumber=Expand[number-PhiPower[newk]];
		 Return[Prepend[RemainderList[newnumber,newk],newk]];
		 )
		 ]
	]

Examples giving norms as well as k's

MakePowerList[3]

MakePowerList[6]

MakePowerList[33]