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.

`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
```

Out[7]:

In [8]:

```
mod(1.78,2π),mod(1.78+2π,2π) # both return approximately the same number
```

Out[8]:

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

This is the same as:

In [10]:

```
z=(-r)*exp(im*(θ+π))
```

Out[10]:

which is the same as:

In [11]:

```
z=(-r)*exp(im*mod(θ+π,2π))
```

Out[11]:

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

*(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.*

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

Note that `Julia`

automatically gives us higher order powers, which default to calling `*`

repeatedly:

In [19]:

```
RadialComplex(2.,2.)^5
```

Out[19]:

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

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

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

*(d)* Implement `/`

for two `Rat`

s, by calling `inv`

and `*`

. The output should be a `Rat`

.

In [24]:

```
import Base./
/(a::Rat,b::Rat)=a*inv(b)
```

Out[24]:

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

`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]:

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

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

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

In [31]:

```
c # c.r has not changed when a.r is changed
```

Out[31]:

Implement `copy(r::Rat)`

.

In [32]:

```
copy(r::Rat)=Rat(r.p,r.q)
```

Out[32]:

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

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

In [35]:

```
Float64(Rat(5,6))
```

Out[35]:

`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 $$

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

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

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

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

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

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

?

In [47]:

```
expn(1.0,20)-e
```

Out[47]:

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

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