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 [ ]:
?mod
In [ ]:
mod(1.78,2π),mod(1.78+2π,2π)  # both return approximately the same number

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 [ ]:
r=-3.5
θ=5.1

z=r*exp(im*θ)

This is the same as:

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

which is the same as:

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

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 [ ]:
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 [ ]:
r=-3.5
θ=5.1

RadialComplex(r,θ)

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 /?

In [ ]:

(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 [ ]:
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    
    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 [ ]:
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.)

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

In [ ]:
RadialComplex(2.,2.)^5

Excercise 2

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

In [ ]:
function *(a::Rat,b::Rat)
    ## TODO: return a Rat that represents a*b
end

(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 [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

Override show

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 [ ]:
import Base.show

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

RadialComplex(2.0,0.1)

Exercise 3

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

In [ ]:

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 the inbuilt function copy to create a copy of a that is completely independent:

In [ ]:
import Base.copy

function copy(z::RadialComplex)
    RadialComplex(z.r,z.θ) 
end

We can compare copy versus =:

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

a.r=2.

b  # b.r has changed when a.r is changed
In [ ]:
c  # c.r has not changed when a.r is changed

Excercise 4

Implement copy(r::Rat).

In [ ]:

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 [ ]:
import Base.Complex128  

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

Complex128(RadialComplex(2.,0.1))

Excercise 5

Override the inbuilt function Float64(r::Rat) to return a Float64 approximation to r. Check that Float64(Rat(5,2)) returns 2.5.

In [ ]:

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) Complete the 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 [ ]:
function expn(z::Rat,n::Integer)
    #TODO: sum n terms of the Taylor series
end

(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. Also, e is a special constant just like pi.

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

(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 [ ]:

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

In [ ]: