Lab 2: Creating a rational type

Please download this file to the lab machine and open it in Jupyter, and complete the exercises in the fields provided. Remember that Jupter can be run by evaluating inside Julia

using IJulia
notebook()

Hint: If you call

@async notebook()

it will run Jupyter in the background, and you can continue to use the Julia session from the terminal.

Creating RadialComplex and Rat types

This lab asks you to design and implement a new composite type Rat to represent rational numbers: that is, numbers of the form $${p \over q}$$ where $p$ and $q$ are integers. We will also investigate why we use floating point arithmetic in place of rational numbers.

To explain how to make a new type, I'll construct a type RadialComplex that represents complex numbers as $z=re^{i\theta}$. I'll then ask you to replicate the definitions for defining a type Rat to represent rational functions.

We begin with the definition of RadialComplex. We want the type to have two fields: r and θ.

We will want to enforce that r and theta satisfy abs(r) ≥ 0 and 0≤θ<2π. The latter is doable by using the command mod(p,q), which will return a number r satisfying 0≤r<q so that p=r+k*q for some integer k:

In [7]:
?mod
search: 
Out[7]:
mod(x, y)

Modulus after division, returning in the range [0,y), if y is positive, or (y,0] if y is negative.

mod modf mod1 module mod2pi Module module_name module_parent chmod

In [8]:
mod(1.78,2π),mod(1.78+2π,2π)  # both return approximately the same number
Out[8]:
(1.78,1.7799999999999994)

We can flip the sign of r to be positive. Note that $-1=e^{i\pi}$. So if $r<0$ we can multiply both terms by $-1$ to get $$re^{i\theta}= (-r)*(-1)*e^{i\theta} = (-r)*e^{i(\theta+\pi)}.$$ Here is an example:

In [9]:
r=-3.5
θ=5.1

z=r*exp(im*θ)
Out[9]:
-1.3229220994954307 + 3.2403513881470634im

This is the same as:

In [10]:
z=(-r)*exp(im*(θ+π))
Out[10]:
-1.3229220994954305 + 3.2403513881470634im

which is the same as:

In [11]:
z=(-r)*exp(im*mod(θ+π,2π))
Out[11]:
-1.3229220994954312 + 3.2403513881470634im

Types with constructors

We are now ready to make the type RadialComplex. If we put a function call called RadialComplex inbetween the type and the end with two arguments, this is called a constructor, and will always be called before the type is created. Calling new with two arguments inside the constructor will create the type.

In [6]:
type RadialComplex
    r::Float64
    θ::Float64
    
    function RadialComplex(r::Float64,θ::Float64)
        if r<0
            # flip the sign of r and add π to θ
            r=-r
            θ=θ+π
        end
        θ=mod(θ,2π)   # ensures θ is in [0,2π)
        new(r,θ)      # creates a new RadialComplex with the updated r and θ
    end
end

We can try it out with the r and θ used above:

In [6]:
r=-3.5
θ=5.1

RadialComplex(r,θ)
Out[6]:
RadialComplex(3.5,1.9584073464102065)

Exercise 1

(a) Type ?gcd and ?div to find out what the commands gcd and div do. Compare 3/6 to div(3,6). What is the difference between div and /?

(b) Complete the type Rat started below. Rat has two fields: p represents the numerator and q represents the denomenator. Both are Int32.

Inside the constructor, ensure that you remove the greater common divisor of p and q before creating the type. Check that Rat(6,4) returns Rat(3,2).

In [15]:
type Rat
    p::Int32
    q::Int32
    function Rat(p,q)
        # TODO: div p and q by their greatest common divisor, and call new to create Rat

        d=gcd(p,q)
        new(div(p,d),div(q,d))        
    end
end

WARNING: You can only create a type once. You may need to restart Julia if you make a mistake in the definition of your type.

Override *

We can override inbuilt commands like *, +, etc. for our special types. For example, for RadialComplex we want to implement the following: $$(r_1 e^{i\theta_1}) *(r_2 e^{i\theta_2})= (r_1r_2) e^{i(\theta_1+\theta_2)}.$$ In other words, we should multiply the radii and add the angles. This can be accomplished as follows:

In [16]:
import Base.*  # This allows us to override the inbuilt * function

function *(a::RadialComplex,b::RadialComplex)
    RadialComplex(a.r*b.r,a.θ+b.θ) 
end

RadialComplex(1.,2.)*RadialComplex(3.,4.)
Out[16]:
RadialComplex(3.0,6.0)

Note that Julia automatically gives us higher order powers, which default to calling * repeatedly:

In [19]:
RadialComplex(2.,2.)^5
Out[19]:
RadialComplex(32.0,3.7168146928204138)

Excercise 2

(a) Implement * for Rat by completing the following function:

In [20]:
function *(a::Rat,b::Rat)
    ## TODO: return a Rat that represents a*b
    Rat(a.p*b.p,a.q*b.q)
end
Out[20]:
* (generic function with 140 methods)

(b) Implement + and - for Rat. Don't forget to import Base.+ and Base.- first, and make sure you do the operations for rationals correctly by making the denominators match.

In [21]:
import Base.+, Base.-

function +(a::Rat,b::Rat)
    ## TODO: return a Rat that represents a*b
    Rat(a.p*b.q + b.p*a.q,a.q*b.q)
end

function -(a::Rat,b::Rat)
    ## TODO: return a Rat that represents a*b
    Rat(a.p*b.q - b.p*a.q,a.q*b.q)
end
Out[21]:
- (generic function with 205 methods)

(c) inv(x) is an inbuilt function that returns $1/x$ for numbers and $x^{-1}$ for matrices. Override Base.inv(r::Rat) to return 1/r, represented as a Rat. Check that inv(Rat(2,3)) returns Rat(3,2).

In [23]:
import Base.inv

inv(a::Rat)=Rat(a.q,a.p)

inv(Rat(2,3))
Out[23]:
Rat(3,2)

(d) Implement / for two Rats, by calling inv and *. The output should be a Rat.

In [24]:
import Base./

/(a::Rat,b::Rat)=a*inv(b)
Out[24]:
/ (generic function with 49 methods)

(e) Check that (Rat(3,4)*Rat(5,6)-Rat(3,5))/Rat(1,40) returns Rat(1,1)

In [25]:
(Rat(3,4)*Rat(5,6)-Rat(3,5))/Rat(1,40)
Out[25]:
Rat(1,1)

Override show

Base.show is an inbuilt function that tells Julia how to display a type. For example, it dictates that Float64(1) should be displayed as 1.0.

Now that RadialComplex is defined, we want to override Base.show to display RadialComplex(r,θ) in a pretty format: r*exp(θ*im). This function takes in a special inbuilt type IO which represents where we want to print to. Passing it as the first argument to print causes the output to be printed to IJulia as desired:

In [26]:
import Base.show

show(io::IO,z::RadialComplex)=print(io,string(z.r)*"exp("*string(z.θ)*"im)")

RadialComplex(2.0,0.1)
Out[26]:
2.0exp(0.1im)

Exercise 3

Override show(io::IO,r::Rat) to display a rational function Rat(p,q) as p/q.

In [28]:
show(io::IO,r::Rat)=print(io,string(r.p)*"/"*string(r.q))
Out[28]:
3/4

Override copy

Remember from lectures that, for composite types, if we use = as in

a=RadialComplex(1.,2.)
b=a

the fields of a point to the same location in memory as the fields in b. Thus modifying a.r will modify b.r and vice-versa.

Therefore, it is useful to override Base.copy to create a copy of a that is completely independent:

In [29]:
import Base.copy

function copy(z::RadialComplex)
    RadialComplex(z.r,z.θ) 
end
Out[29]:
copy (generic function with 33 methods)

We can compare copy versus =:

In [30]:
a=RadialComplex(1.,2.)
b=a
c=copy(a)

a.r=2.

b  # b.r has changed when a.r is changed
Out[30]:
2.0exp(2.0im)
In [31]:
c  # c.r has not changed when a.r is changed
Out[31]:
1.0exp(2.0im)

Excercise 4

Implement copy(r::Rat).

In [32]:
copy(r::Rat)=Rat(r.p,r.q)
Out[32]:
copy (generic function with 34 methods)

Converting to inbuilt number types

We want to be able to convert between our new types and inbuilt number types. For example, the following routine converts a RadialComplex to a regular Complex128:

In [33]:
import Base.Complex128  

Complex128(z::RadialComplex)=z.r*exp(im*z.θ)

Complex128(RadialComplex(2.,0.1))
Out[33]:
1.9900083305560516 + 0.1996668332936563im

Excercise 5

Override Base.Float64(r::Rat) to return a Float64 approximation to r.

In [34]:
import Base.Float64

Float64(r::Rat)=r.p/r.q
Out[34]:
Float64
In [35]:
Float64(Rat(5,6))
Out[35]:
0.8333333333333334

Rat versus Float64

Unlike Float64, Rat has no error for algebraic operations. Why do we use floating point representations then? The answer is that Rat is highly susceptible to overflow. The following excercise demonstrates the problem.

Recall the Taylor expansion of $e^z$: $$ e^z=1 + z+{z^2 \over 2} + {z^3 \over 3!} + {z^4 \over 4!} + \cdots $$

Excercise 6

(a) Write a function expn(z::Rat,n::Integer) that finds a rational approximation to exp(z) by summing the first n terms of the Taylor series. Hint: use the inbuilt function factorial.

In [36]:
function expn(z::Rat,n::Integer)
    r=Rat(1,1)
    for k=2:n
        r+=z^(k-1)*Rat(1,factorial(k-1))
    end
    r
end
Out[36]:
expn (generic function with 1 method)

(b) What is the accuracy of expn(Rat(1,1),9), which should approximate e? Remember, you can use your previous routine Float64(r::Rat) to convert to a Float64.

In [41]:
Float64(expn(Rat(1,1),9))-e
Out[41]:
-3.0586177750535626e-6

(c) What is the accuracy of expn(Rat(1,1),10)? What has gone wrong?

In [43]:
Float64(expn(Rat(1,1),10))-e
Out[43]:
-2.0784897278989707

(d) Write expn(z::Float64,n::Integer) that uses floating point arithmetic to sum the first n terms of the Taylor series.

In [44]:
function expn(z::Float64,n::Integer)
    r=1.0
    for k=2:n
        r+=z^(k-1)/factorial(k-1)
    end
    r
end
Out[44]:
expn (generic function with 2 methods)

(d) Does expn(1.0,10) suffer the same inaccuracies as expn(Rat(1,1),10)? Why not?

In [46]:
expn(1.0,10)-e
Out[46]:
-3.0288585284310443e-7

(e) What is the error of expn(1.,20)?

In [47]:
expn(1.0,20)-e
Out[47]:
4.440892098500626e-16

(f) (Advanced) If you used factorial, then this approach will fail for expn(1.,30). Rewrite the expn to build up 1/k! by hand using floating point multiplication.

In [50]:
function expn(z::Float64,n::Integer)
    r=1.0
    fac=1.0
    for k=2:n
        r+=z^(k-1)*fac
        fac/=k
    end
    r
end
Out[50]:
expn (generic function with 2 methods)

(g) What happens to the accuracy of expn(1.0,n) as n becomes large?
Can you explain why the accuracy is limited?

In [57]:
expn(1.0,100)-e
Out[57]:
4.440892098500626e-16