```python
%matplotlib inline
```


```python
from Crypto.Util import number
from sympy.ntheory import factorint
from turtle import *
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import numpy.polynomial.polynomial as p
import time
import random
import hashlib
import sys
```

# Basic Algebra Exercise
## Functions, Polynomials, Complex Numbers. Applications of Abstract Algebra

### Problem 1. Polynomial Interpolation
We know that if we have a set of $n$ data points with coordinates $(x_1; y_1), (x_2; y_2), \dots, (x_n; y_n)$, we can try to figure out what function may have generated these points.

Please note that **our assumptions about the data** will lead us to choosing one function over another. This means that our results are as good as our data and assumptions. Therefore, it's extremely important that we write down our assumptions (which sometimes can be difficult as we sometimes don't realize we're making them). It will be better for our readers if they know what those assumptions and models are.

In this case, we'll state two assumptions:
1. The points in our dataset are generated by a polynomial function
2. The points are very precise, there is absolutely no error in them. This means that the function should pass **through every point**

This method is called *polynomial interpolation* (*"polynomial"* captures assumption 1 and *"interpolation"* captures assumption 2).

It can be proved (look at [Wikipedia](https://en.wikipedia.org/wiki/Polynomial_interpolation) for example) that if we have $n$ data points, there is only one polynomial of degree $n-1$ which passes through them. In "math speak": "the vector spaces of $n$ points and polynomials of degree $n-1$ are isomorphic (there exists a bijection mapping one to the other)".

There are a lot of ways to do interpolation. We can also write the function ourselves if we want but this requires quite a lot more knowledge than we already covered in this course. So we'll use a function which does this for us. `numpy.polyfit()` is one such function. It accepts three main parameters (there are others as well, but they are optional): a list of $x$ coordinates, a list of $y$ coordinates, and a polynomial degree.

Let's say we have these points:
```python
points = np.array([(0, 0), (1, 0.8), (2, 0.9), (3, 0.1), (4, -0.8), (5, -1.0)])
```

First, we need to "extract" the coordinates:
```python
x = points[:, 0]
y = points[:, 1]
```

Then, we need to calculate the interpolating polynomial. For the degree, we'll set $n-1$:
```python
coefficients = np.polyfit(x, y, len(points) - 1)
poly = np.poly1d(coefficients)
```

After that, we need to plot the function. To do this, we'll create a range of $x$ values and evaluate the polynomial at each value:
```python
plot_x = np.linspace(np.min(x), np.max(x), 1000)
plot_y = poly(plot_x)
```

Finally, we need to plot the result. We'll plot both the fitting polynomial curve (using `plt.plot()`) and the points (using `plt.scatter`). It's also nice to have different colors to make the line stand out from the points.
```python
plt.plot(plot_x, plot_y, c = "green")
plt.scatter(x, y)
plt.xlabel("x")
plt.ylabel("y")
plt.show()
```
Don't forget to label the axes!

Your task now is to **wrap the code in a function**. It should accept a list of points, the polynomial degree, min and max value of $x$ used for plotting. **Be extremely careful to ensure that the function uses its parameters!**

We'll use this function to try some other cases.


```python
def center_spines(ax, centerx=0, centery=0):
    """Centers the axis spines at <centerx, centery> on the axis 'ax' """

    # Set the axis's spines to be centered at the given point
    # (Setting all 4 spines so that the tick marks go in both directions)
    ax.spines['left'].set_position(('data', centerx))
    ax.spines['bottom'].set_position(('data', centery))
    ax.spines['right'].set_position(('data', centerx - 1))
    ax.spines['top'].set_position(('data', centery - 1))

    # Hide the line (but not ticks) for "extra" spines
    for side in ['right', 'top']:
        ax.spines[side].set_color('none')

    # On both the x and y axes
    for axis, center in zip([ax.xaxis, ax.yaxis], [centerx, centery]):
        # Hide the ticklabels at <centerx, centery>
        formatter = CenteredFormatter()
        formatter.center = center
        axis.set_major_formatter(formatter)

    # Add offset ticklabels at <centerx, centery> using annotation
    # (Should probably make these update when the plot is redrawn...)
    xlabel, ylabel = map(formatter.format_data, [centerx, centery])
    ax.annotate('%s' % xlabel, (centerx, centery),
            xytext=(-2, -4), textcoords='offset points',
            ha='right', va='top')

    ax.xaxis.set_ticks_position('bottom')
    ax.yaxis.set_ticks_position('left')

# Custom formatter
class CenteredFormatter(mpl.ticker.ScalarFormatter):
    """Acts exactly like the default Scalar Formatter, but yields an empty
    label for ticks at "center"."""
    center = 0
    def __call__(self, value, pos=None):
        if value == self.center:
            return ''
        else:
            return mpl.ticker.ScalarFormatter.__call__(self, value, pos)
```


```python
def get_x_y_points(points):
    x = points[:, 0]
    y = points[:, 1]
    return x, y

def get_min_max_points(points):
    min_points = np.min(points[:, 0])
    max_points = np.max(points[:, 0])
    return min_points, max_points

def get_degrees(points):
    return len(points) - 1

def interpolate_polynomial(points, degree, min_x, max_x):
    """
    Interpolates a polynomial of the specified degree through the given points and plots it
    points - a list of points (x, y) to plot
    degree - the polynomial degree
    min_x, max_x - range of x values used to plot the interpolating polynomial
    """
    x, y = get_x_y_points(points)
    # x = points[:, 0]
    # y = points[:, 1]
    coefficients = np.polyfit(x, y, degree)
    poly = np.poly1d(coefficients)
    plot_x = np.linspace(min_x, max_x, 5000)
    plot_y = poly(plot_x)
    return x, y, plot_x, plot_y
    
def get_plot_graph(x, y, plot_x, plot_y, plot_title):
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.plot(plot_x, plot_y, c = "r")
    plt.scatter(x, y, c = "g")
    # plt.xticks(np.arange(-2, 6, 0.5))
    # plt.yticks(np.arange(-2, 2, 0.5))
    plt.xlabel("x", loc = 'right', rotation = 0)
    plt.ylabel("y", loc = 'top', rotation = 0)
    plt.title(plot_title, pad=20)
    # ax.set_aspect('equal')
    center_spines(ax)
    return plt.show()

```


```python
points = np.array([(0, 0), (1, 0.8), (2, 0.9), (3, 0.1), (4, -0.8), (5, -1.0)])
min_points, max_points = get_min_max_points(points)
degree = get_degrees(points)
x, y, plot_x, plot_y = interpolate_polynomial(points, degree, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 4 - interpolation")
```

We see this is a very nice fit. This is expected, of course. Let's try to expand our view a little. Let's try to plot other values of $x$, further than the original ones. This is **extrapolation**.


```python
x, y, plot_x, plot_y = interpolate_polynomial(points, degree, -5, 10)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 4 - extrapolation")
```

Hmmm... it seems our polynomial goes a little wild outside the original range. This is to show how **extrapolation can be quite dangerous**.

Let's try a lower polynomial degree now. We used 4, how about 3, 2 and 1?
**Note:** We can add titles to every plot so that we know what exactly we're doing. Te title may be passed as an additional parameter to our function.


```python
x, y, plot_x, plot_y = interpolate_polynomial(points, 3, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 3")
x, y, plot_x, plot_y = interpolate_polynomial(points, 2, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 2")
x, y, plot_x, plot_y = interpolate_polynomial(points, 1, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 1")
```

We see the fitting curves (or line in the last case) struggle more and more and they don't pass through every point. This breaks our assumptions but it can be very useful.

Okay, one more thing. How about increasing the degree? Let's try 5, 7 and 10. Python might complain a little, just ignore it, everything is fine... sort of :).


```python
x, y, plot_x, plot_y = interpolate_polynomial(points, 5, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 5")
x, y, plot_x, plot_y = interpolate_polynomial(points, 7, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 7")
x, y, plot_x, plot_y = interpolate_polynomial(points, 10, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 10")
```

Those graphs look pretty much the same. But that's the point exactly. I'm being quite sneaky here. Let's try to expand our view once again and see what our results really look like.


```python
x, y, plot_x, plot_y = interpolate_polynomial(points, 5, -10, 10)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 5 - extrapolation")
x, y, plot_x, plot_y = interpolate_polynomial(points, 7, -10, 10)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 5 - extrapolation")
x, y, plot_x, plot_y = interpolate_polynomial(points, 10, -10, 10)
get_plot_graph(x, y, plot_x, plot_y, "Polynomial of degree 5 - extrapolation")
```

Now we see there are very wild differences. Even though the first two plots look quite similar, look at the $y$ values - they're quite different.

So, these are the dangers of interpolation. Use a too high degree, and you get "the polynomial wiggle". These are all meant to represent **the same** data points but they look insanely different. Here's one more comparison.


```python
x, y, plot_x, plot_y = interpolate_polynomial(points, len(points) - 1, -2, 7)
get_plot_graph(x, y, plot_x, plot_y, "The polynomial wiggle 1")
x, y, plot_x, plot_y = interpolate_polynomial(points, len(points) + 1, -2, 7)
get_plot_graph(x, y, plot_x, plot_y, "The polynomial wiggle 2")
```

Now we can see what big difference even a small change in degree can make. This is why we have to choose our interpolating functions very carefully. Generally, a lower degree means a simpler function, which is to be preferred. See [Occam's razor](https://en.wikipedia.org/wiki/Occam%27s_razor).

And also, **we need to be very careful about our assumptions**.


```python
points = np.array([(-5, 0.03846), (-4, 0.05882), (-3, 0.1), (-2, 0.2), (-1, 0.5), (0, 1), (1, 0.5), (2, 0.2), (3, 0.1), (4, 0.05882), (5, 0.03846)])
min_points, max_points = get_min_max_points(points)
degree = get_degrees(points)
x, y, plot_x, plot_y = interpolate_polynomial(points, degree, min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polinomial with lower degree")
```

This one definitely looks strange. Even stranger, if we remove the outermost points... ($x = \pm 5$), we get this


```python
points = np.array([(-4, 0.05882), (-3, 0.1), (-2, 0.2), (-1, 0.5), (0, 1), (1, 0.5), (2, 0.2), (3, 0.1), (4, 0.05882)])
min_points, max_points = get_min_max_points(points)
x, y, plot_x, plot_y = interpolate_polynomial(points, len(points - 1), min_points, max_points)
get_plot_graph(x, y, plot_x, plot_y, "Polinomial with lower degree")
```

This is because the generating function is not a polynomial. It's actually:
$$ y = \frac{1}{1 + x^2} $$

Plot the polynomial interpolation and the real generating function **on the same plot**. You may need to modify the original plotting function or just copy its contents.


```python
def get_generating_function_data(plot_x):
    return 1/(1 + plot_x ** 2)
    
def get_plot_graph(x, y, plot_x, ys, plot_title, labels, colors):
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    counter = 0
    for plot_y in ys:
        plt.plot(plot_x, plot_y, label = labels[counter], c = colors[counter])
        counter += 1
    plt.scatter(x, y, c = "g")
    # plt.xticks(np.arange(-2, 6, 0.5))
    # plt.yticks(np.arange(-2, 2, 0.5))
    plt.xlabel("x", loc = 'right', rotation = 0)
    plt.ylabel("y", loc = 'top', rotation = 0)
    plt.title(plot_title, pad=20)
    plt.legend(loc='upper right')
    # ax.set_aspect('equal')
    center_spines(ax)
    return plt.show()
```


```python
points = np.array([(0, 0), (1, 0.8), (2, 0.9), (3, 0.1), (4, -0.8), (5, -1.0)])
#points = np.array([(-4, 0.05882), (-3, 0.1), (-2, 0.2), (-1, 0.5), (0, 1), (1, 0.5), (2, 0.2), (3, 0.1), (4, 0.05882)])
min_points, max_points = get_min_max_points(points)
degree = get_degrees(points)
x, y, plot_x, plot_y = interpolate_polynomial(points, degree, min_points, max_points)
y_gen = get_generating_function_data(plot_x)
ys = [plot_y, y_gen]
title = "Polynomial Interpolation and generating function"
labels = ['Polinomial', r'$y = \frac{1}{1 + x^2}$']
colors = ['r', 'lightblue']
get_plot_graph(x, y, plot_x, ys, title, labels, colors)
```

### Problem 2. Complex Numbers as Vectors
We saw that a complex number $z = a + bi$ is equivalent to (and therefore can be represented as) the ordered tuple $(a; b)$, which can be plotted in a 2D space. So, complex numbers and 2D points are equivalent. What is more, we can draw a vector from the origin of the coordinate plane to our point. This is called a point's **radius-vector**.

Let's try plotting complex numbers as radius vectors. Don't forget to label the real and imaginary axes. Also, move the axes to the origin. Hint: These are called "spines"; you'll need to move 2 of them to the origin and remove the other 2 completely. Hint 2: You already did this in the previous lab.

We can use `plt.quiver()` to plot the vector. It can behave a bit strangely, so we'll need to set the scale of the vectors to be the same as the scale on the graph axes:
```python
plt.quiver(0, 0, z.real, z.imag, angles = "xy", scale_units = "xy", scale = 1)
```

Other than that, the main parameters are: $x_{begin}$, $y_{begin}$, $x_{length}$, $y_{length}$ in that order.

Now, set the aspect ratio of the axes to be equal. Also, add grid lines. Set the axis numbers (called ticks) to be something like `range(-3, 4)` for now.
```python
plt.xticks(range(-3, 4))
plt.yticks(range(-3, 4))

```

If you wish to, you can be a bit more clever with the tick marks. Find the minimal and maximal $x$ and $y$ values and set the ticks according to them. It's a good practice not to jam the plot too much, so leave a little bit of space. That is, if the actual x-range is $[-2; 2]$, set the plotting to be $[-2.5; 2.5]$ for example. Otherwise, the vector heads (arrows) will be "jammed" into a corner or side of the plot.


```python
def plot_complex_number(z):
    """
    Plots the complex number z as a radius vector in the 2D space
    """
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.quiver(0, 0, z.real, z.imag, angles = "xy", scale_units = "xy", scale = 1, color = 'r')
    plt.xticks(range(-3, 4))
    plt.yticks(range(-3, 4))
    plt.xlabel("Re", loc = 'right')
    plt.ylabel("Imag", loc = 'top', rotation=0)
    plt.title("Complex Number", pad=20)
    ax.set_aspect('equal')
    center_spines(ax)
    return plt.show()

plot_complex_number(2 + 3j)
```

How about many numbers? We'll need to get a little bit more creative. First, we need to create a 2D array, each element of which will be a 4-element array: `[0, 0, z.real, z.imag]`. Next, `plt.quiver()` can accept a range of values. Look at [this StackOverflow post](https://stackoverflow.com/questions/12265234/how-to-plot-2d-math-vectors-with-matplotlib) for details and adapt your code.


```python
def plot_complex_numbers(numbers, colors):
    """
    Plots the given complex numbers as radius vectors in the 2D space
    """
    zeros = [0 for number in numbers]
    real_parts = [number.real for number in numbers]
    imaginary_parts = [number.imag for number in numbers]
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.quiver(zeros, zeros, real_parts, imaginary_parts, angles = "xy", scale_units = "xy", scale = 1, color = colors)
    plt.xticks(range(-3, 5))
    plt.yticks(range(-3, 5))
    plt.xlabel("Re", loc = 'right')
    plt.ylabel("Imag", loc = 'top', rotation=0)
    plt.title("Complex Number", pad=20)
    ax.set_aspect('equal')
    center_spines(ax)
    return plt.show()
plot_complex_numbers([2 + 3j, -2 - 1j, -3, 2j], ["green", "red", "blue", "orange"])
```

Now let's see what the operations look like. Let's add two numbers and plot the result.


```python
z1 = 2 + 3j
z2 = 1 - 1j
plot_complex_numbers([z1, z2, z1 + z2], ["red", "blue", "green"])
```

We can see that adding the complex numbers is equivalent to adding vectors (remember the "parallelogram rule"). As special cases, let's try adding pure real and pure imaginary numbers:


```python
z1 = 2 + 3j
z2 = 2 + 0j
plot_complex_numbers([z1, z2, z1 + z2], ["red", "blue", "green"])
```


```python
def plot_complex_numbers(numbers, colors):
    """
    Plots the given complex numbers as radius vectors in the 2D space
    """
    zeros = [0 for number in numbers]
    real_parts = [number.real for number in numbers]
    imaginary_parts = [number.imag for number in numbers]
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.quiver(zeros, zeros, real_parts, imaginary_parts, angles = "xy", scale_units = "xy", scale = 1, color = colors)
    plt.xticks(range(-3, 5))
    plt.yticks(range(-3, 6))
    plt.xlabel("Re", loc = 'right')
    plt.ylabel("Imag", loc = 'top', rotation=0)
    plt.title("Complex Number", pad=20)
    ax.set_aspect('equal')
    center_spines(ax)
    return plt.show()
z1 = 2 + 3j
z2 = 0 + 2j
plot_complex_numbers([z1, z2, z1 + z2], ["red", "blue", "green"])
```

How about multiplication? First we know that multiplying by 1 gives us the same vector and mulpiplying by -1 gives us the reversed version of the same vector. How about multiplication by $\pm i$?


```python
def plot_complex_numbers(numbers, colors):
    """
    Plots the given complex numbers as radius vectors in the 2D space
    """
    zeros = [0 for number in numbers]
    real_parts = [number.real for number in numbers]
    imaginary_parts = [number.imag for number in numbers]
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.quiver(zeros, zeros, real_parts, imaginary_parts, angles = "xy", scale_units = "xy", scale = 1, color = colors)
    plt.xticks(range(-3, 5))
    plt.yticks(range(-3, 5))
    plt.xlabel("Re", loc = 'right')
    plt.ylabel("Imag", loc = 'top', rotation=0)
    plt.title("Complex Number multiplication $\pm i$", pad=20)
    ax.set_aspect('equal')
    center_spines(ax)
    return plt.show()
z = 2 + 3j
plot_complex_numbers([z, z * 1], ["red", "blue"])
plot_complex_numbers([z, z * -1], ["red", "blue"])
plot_complex_numbers([z, z * 1j], ["red", "blue"])
plot_complex_numbers([z, z * -1j], ["red", "blue"])
```

So, multiplication by $i$ is equivalent to 90-degree rotation. We can actually see the following equivalence relationships between multiplying numbers and rotation about the origin:

| Real | Imaginary | Result rotation |
|------|-----------|-----------------|
| 1    | 0         | $0^\circ$       |
| 0    | 1         | $90^\circ$      |
| -1   | 0         | $180^\circ$     |
| 0    | -1        | $270^\circ$     |

Once again, we see the power of abstraction and algebra in practice. We know that complex numbers and 2D vectors are equivalent. Now we see something more: addition and multiplication are equivalent to translation (movement) and rotation!

Let's test the multiplication some more. We can see the resulting vector is the sum of the original vectors, but *scaled and rotated*:


```python
def plot_complex_numbers(numbers, colors):
    """
    Plots the given complex numbers as radius vectors in the 2D space
    """
    zeros = [0 for number in numbers]
    real_parts = [number.real for number in numbers]
    imaginary_parts = [number.imag for number in numbers]
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.quiver(zeros, zeros, real_parts, imaginary_parts, angles = "xy", scale_units = "xy", scale = 1, color = colors)
    plt.xticks(range(-3, 10))
    plt.yticks(range(-3, 5))
    plt.xlabel("Re", loc = 'right')
    plt.ylabel("Imag", loc = 'top', rotation=0)
    plt.title("Complex Number translation (movement) and rotation", pad=20)
    ax.set_aspect('equal')
    center_spines(ax)
    return plt.show()
z1 = 2 + 3j
z2 = 1 - 2j
plot_complex_numbers([z1, z2, z1 * z2], ["red", "blue", "green"])
```

### Problem 3. Recursion and Fractals


> "To understand recursion, you first need to understand recursion."

There are three main parts to a recursive function:
1. Bottom - when the recursion should finish
2. Operation - some meaningful thing to do
3. Recursive call - calling the same function
4. Clean-up - returning all data to its previous state (this reverses the effect of the operation)

Let's do one of the most famous recursion examples. And I'm not talking about Fibonacci here. Let's draw a tree using recursive functions.

The figure we're going to draw is called a **fractal**. It's self-similar, which means that if you zoom in on a part of it, it will look the same. You can see fractals everywhere in nature, with broccoli being one of the prime examples. Have a look:



First, we need to specify the recursive part. In order to draw a tree, we need to draw a line of a given length (which will be the current branch), and then draw two more lines to the left and right. By "left" and "right", we should mean "rotation by a specified angle".

So, this is how to draw a branch: draw a line and prepare to draw two more branches to the left and right. This is going to be our recursive call.

To make things prettier, more natural-looking (and have a natural end to our recursion), let's draw each "sub-branch" a little shorter. If the branch becomes too short, it won't have "child branches". This will be the bottom of our recursion.

There's one more important part of recursion, and this is **"the clean-up"**. After we did something in the recursive calls, it's very important to return the state of everything as it was **before** we did anything. In this case, after we draw a branch, we go back to our starting position. 

Let's first import the most import-ant (no pun intended...) Python drawing library: `turtle`! In order to make things easier, we'll import all methods directly.
```python
from turtle import *
```

You can look up the docs about turtle if you're more interested. The basic things we're going to use are going forward and backward by a specified number of pixels and turning left and right by a specified angle (in degrees).

Let's now define our recursive function:
```python
def draw_branch(branch_length, angle):
    if branch_length > 5:
        forward(branch_length)
        right(angle)
        draw_branch(branch_length - 15, angle)
        left(2 * angle)
        draw_branch(branch_length - 15, angle)
        right(angle)
        backward(branch_length)
```

And let's call it:
```python
draw_branch(100, 20)
```

We need to start the tree not at the middle, but toward the bottom of the screen, so we need to make a few more adjustments. We can wrap the setup in another function and call it. Let's start one trunk length below the center (the trunk length is the length of the longest line).

**Note:** It's important to call `done()` after the drawing is finished. If you miss it, the window with the turtle drawing will freeze and throw an exception.

```python
def draw_tree(trunk_length, angle):
    speed("fastest")
    left(90)
    up()
    backward(trunk_length)
    down()
    draw_branch(trunk_length, angle)
    done()
```

Note that the graphics will show in a separate window. Also note that sometimes you might get bugs. If you do, go to Kernel > Restart.


```python
def draw_branch(branch_length, angle):
    if branch_length > 5:
        forward(branch_length)
        right(angle)
        draw_branch(branch_length - 15, angle)
        left(2 * angle)
        draw_branch(branch_length - 15, angle)
        right(angle)
        backward(branch_length)
        

def draw_tree(trunk_length, angle):
    color('brown')
    shape('turtle')
    speed("fastest")
    left(90)
    up()
    backward(trunk_length)
    down()
    draw_branch(trunk_length, angle)
    done()
```


```python
draw_branch(100, 20)
```


```python
draw_tree(100, 20)
```

Experiment with different lengths and angles. Especially interesting angles are $30^\circ$, $45^\circ$, $60^\circ$ and $90^\circ$.


```python
draw_tree(100, 30)
```


```python
draw_tree(100, 45)
```


```python
draw_tree(100, 90)
```

Now modify the original function a little. Draw the lines with different thickness. Provide the trunk thickness at the initial call. Similar to how branches go shorter, they should also go thinner.


```python
def draw_branch(branch_length, angle):
    size = int(branch_length / 10)
    pensize(size)
    if branch_length > 5:
        forward(branch_length)
        right(angle)
        draw_branch(branch_length - 15, angle)
        left(2 * angle)
        draw_branch(branch_length - 15, angle)
        right(angle)
        backward(branch_length)
```


```python
draw_tree(100, 90)
```

#### * Optional problem
Try to draw another kind of fractal graphic using recursion and the `turtle` library. Two very popular examples are the "Koch snowflake" and the "Sierpinski triangle". You can also modify the original tree algorithm to create more natural-looking trees. You can, for example, play with angles, number of branches, lengths, and widths. The Internet has a lot of ideas about this :). Hint: Look up **"L-systems"**.

### Problem 4. Run-length Encoding
One application of algebra and basic math can be **compression**. This is a way to save data in less space than it originally takes. The most basic form of compression is called [run-length encoding](https://en.wikipedia.org/wiki/Run-length_encoding).

Write a function that encodes a given text. Write another one that decodes.

We can see that RLE is not very useful in the general case. But it can be extremely useful if we have very few symbols. An example of this can be DNA and protein sequences. DNA code, for example, has only 4 characters.

Test your encoding and decoding functions on a DNA sequence (you can look up some on the Internet). Measure how much your data is compressed relative to the original.


```python
def encode(text):
    """
    Returns the run-length encoded version of the text
    (numbers after symbols, length = 1 is skipped)
    """
    if len(text) == 1:
        return text
    else:
        counter = 1
        result = ''
        for i in range(len(text)):
            # Check if text[i] is last element of the input text
            if i == len(text) - 1:
                if text[i] != text[i-1]:
                    result += text[i]
                else: 
                    result += text[i] + str(counter) 
            # Check if current and next elements are equal
            elif text[i] == text[i+1]:
                counter += 1
            # Check if current and next elements are not equal
            else:
                if counter == 1:
                    result += text[i]
                else:
                    result += text[i] + str(counter)
                    counter = 1
        return result
            

def decode(text):
    """
    Decodes the text using run-length encoding
    """
    # Check if the length of the string is equal to 1
    if len(text) == 1:
        return text
    else:
        result = ''
        for i in range(len(text)):
            # Check if current element is digit - just skip it
            if text[i].isdigit():
                pass
            # Check if this is the last element and add it to the result if it's not a digit
            elif i == len(text) - 1:
                if not text[i].isdigit():
                    result += text[i] 
            # Check if the current element is not digit and next element is digit - add element multiplied by digit to the result
            elif text[i + 1].isdigit():
                result += text[i] * int(text[i + 1])
            # If current and next elements are not digits - add current element to the result
            else:
                result += text[i]
        return result
```


```python
# Tests
# Test that the functions work on their own
assert encode("AABCCCDEEEE") == "A2BC3DE4"
assert decode("A2BC3DE4") == "AABCCCDEEEE"

# Test that the functions really invert each other
assert decode(encode("AABCCCDEEEE")) == "AABCCCDEEEE"
assert encode(decode("A2BC3DE4")) == "A2BC3DE4"
```

### * Problem 5. Function Invertibility and Cryptography
As we already saw, some functions are able to be inverted. That is, if we know the output, we can see what input generated it directly. This is true if the function is **one-to-one correspondence** (bijection).

However, not all functions are created the same. Some functions are easy to compute but their inverses are extremely difficult. A very important example is **number factorization**. It's relatively easy (computationally) to multiply numbers but factoring them is quite difficult. Let's run an experiment.

We'll need a function to generate random n-bit numbers. One such number can be found in the `Crypto` package
```python
from Crypto.Util import number
random_integer = number.getRandomNBitInteger(n_bits)
```

We could, of course, write our factorization by hand but we'll use `sympy`
```python
from sympy.ntheory import factorint
factorint(1032969399047817906432668079951) # {3: 2, 79: 1, 36779: 1, 7776252885493: 1, 5079811103: 1}
```

This function returns a `dict` where the keys are the factors, and the values - how many times they should be multiplied.

We'll also need a tool to accurately measure performance. Have a look at [this one](https://docs.python.org/3/library/time.html#time.time) for example.

Specity a sequence of bit lengths, in increasing order. For example, you might choose something like `[10, 20, 25, 30, 32, 33, 35, 38, 40]`. Depending on your computer's abilities you can go as high as you want. For each bit length, generate a number. See how much time it takes to factor it. Then see how much time it takes to multiply the factors. Be careful how you measure these. You shouldn't include the number generation (or any other external functions) in your timing.

In order to have better accuracy, don't do this once per bit length. Do it, for example, five times, and average the results.

Plot all multiplication and factorization times as a function of the number of bits. You should see that factorization is much, much slower. If you don't see this, just try larger numbers :D.


```python
def multiply_factorizing_time_calculations(bit_lengths):
    random_nums = []
    m_times = []
    f_times = []
    for bit in bit_lengths:
        random_nums.append(number.getRandomNBitInteger(bit))
    result = 1
    print(random_nums)
    for num in random_nums:
        f_start_time = time.time()
        factorizing_data = factorint(num)
        f_end_time = time.time()
        f_duration = f_end_time - f_start_time
        f_times.append(f_duration)
        
        m_start_time = time.time()
        for num, power in factorizing_data.items():
            result *= num ** power
        m_end_time = time.time()
        m_duration = m_end_time - m_start_time
        m_times.append(m_duration)
    return m_times, f_times

def plot_multiply_factorizing_times(bit_lengths, m_times, f_times):
    ax = plt.gca()
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.xlabel("Bits", loc = 'right')
    plt.ylabel("Time, s", loc = 'top', rotation=0)
    plt.title("Multiply and Factorizing Times", pad=20)
    #ax.set_aspect('equal')
    center_spines(ax)
    
    plt.plot(bit_lengths, f_times, label = 'Factorizing')
    plt.plot(bit_lengths, m_times, label = 'Multiply')
    plt.legend(loc='upper left')
    return plt.show()
```


```python
bit_lengths = [10, 20, 25, 30, 32, 33, 35, 38, 40, 100, 125, 150]
m_times, f_times = multiply_factorizing_time_calculations(bit_lengths)
plot_multiply_factorizing_times(bit_lengths, m_times, f_times)
```

### * Problem 6. Diffie - Hellman Simulation
As we already saw, there are functions which are very easy to compute in the "forward" direction but really difficult (computationally) to invert (that is, determine the input from the output). There is a special case: the function may have a hidden "trap door". If you know where that door is, you can invert the function easily. This statement is at the core of modern cryptography.

Look up **Diffie - Hellman key exchange** (here's a [video](https://www.youtube.com/watch?v=cM4mNVUBtHk) on that but feel free to use anything else you might find useful).

Simulate the algorithm you just saw. Generate large enough numbers so the difference is noticeable (say, factoring takes 10-15 seconds). Simulate both participants in the key exchange. Simulate an eavesdropper.

First, make sure after both participants run the algotihm, they have *the same key* (they generate the same number).

Second, see how long it takes for them to exchange keys.

Third, see how long it takes the eavesdropper to arrive at the correct shared secret.

You should be able to see **the power of cryptography**. In this case, it's not that the function is irreversible. It can be reversed, but it takes a really long time (and with more bits, we're talking billions of years). However, if you know something else (this is called a **trap door**), the function becomes relatively easy to invert.


```python
g=9
p=1001

a=random.randint(5, 10)

b=random.randint(10,20)

A = (g**a) % p
B = (g**b) % p

print('g: ',g,' (a shared value), n: ',p, ' (a prime number)')

print('\nAlice calculates:')
print('a (Alice random): ',a)
print('Alice value (A): ',A,' (g^a) mod p')

print('\nBob calculates:')
print('b (Bob random): ',b)
print('Bob value (B): ',B,' (g^b) mod p')

print('\nAlice calculates:')
keyA=(B**a) % p
print('Key: ',keyA,' (B^a) mod p')
print('Key: ',hashlib.sha256(str(keyA).encode()).hexdigest())

print('\nBob calculates:')
keyB=(A**b) % p
print('Key: ',keyB,' (A^b) mod p')
print('Key: ',hashlib.sha256(str(keyB).encode()).hexdigest())
```

    g:  9  (a shared value), n:  1001  (a prime number)
    
    Alice calculates:
    a (Alice random):  6
    Alice value (A):  911  (g^a) mod p
    
    Bob calculates:
    b (Bob random):  16
    Bob value (B):  9  (g^b) mod p
    
    Alice calculates:
    Key:  911  (B^a) mod p
    Key:  a5ccb1c538e34663a658b1be28b16455ee5285efb10e6f1d4caba1f69ec9782b
    
    Bob calculates:
    Key:  911  (A^b) mod p
    Key:  a5ccb1c538e34663a658b1be28b16455ee5285efb10e6f1d4caba1f69ec9782b


### ** Problem 7. The Galois Field in Cryptography
Research about the uses of the Galois field. What are its properties? How can it be used in cryptography? Write a simple cryptosystem based on the field.

You can use the following questions to facilitate your research:
* What is a field?
* What is GF(2)? Why is it an algebraic field?
* What is perfect secrecy? How does it relate to the participants in the conversation, and to the outside eavesdropper?
* What is symmetrical encryption?
* How to encrypt one-bit messages?
* How to extend the one-bit encryption system to many buts?
* Why is the system decryptable? How do the participants decrypt the encrypted messages?
* Why isn't the eavesdropper able to decrypt?
* What is a one-time pad?
    * How does the one-time pad achieve perfect secrecy?
* What happens if we try to use a one-time pad many times?
    * Provide an example where you break the "many-time pad" security
* What are some current enterprise-grade applications of encryption over GF(2)?
* Implement a cryptosystem based on GF(2). Show correctness on various test cases

### ** Problem 8. Huffman Compression Algorithm
Examine and implement the **Huffman algorithm** for compressing data. It's based on information theory and probiability theory. Document your findings and provide your implementation.

This algorithm is used for **lossless compression**: compressing data without loss of quality. You can use the following checklist:

* What is the difference betwenn lossless and lossy compression?
* When can we get away with lossy compression?
* What is entropy?
* How are Huffman trees constructed?
    * Provide a few examples
* How can we get back the uncompressed data from the Huffman tree?
* How and where are Huffman trees stored?
* Implement the algorithm. Add any other formulas / assumptions / etc. you might need.
* Test the algorithm. A good meaure would be percentage compression: $$\frac{\text{compressed}}{\text{uncompressed}} * 100\%$$
* How well does Huffman's algorithm perform compared to other compression algorithms (e.g. LZ77)?

a069299145334f15512289863622c452ad88c006
ipynb
Jupyter Notebook
basic_algebra/Basic Algebra Exercise.ipynb
PetkoAndreev/Python-Math-Concepts-For-Developers
756639ff35bc81e1050f850dd30e76948da290d9
MIT
2021-05-27T07:57:24.000Z
2021-05-27T07:57:24.000Z
basic_algebra/Basic Algebra Exercise.ipynb
PetkoAndreev/Python-Math-Concepts-For-Developers
756639ff35bc81e1050f850dd30e76948da290d9
MIT
basic_algebra/Basic Algebra Exercise.ipynb
PetkoAndreev/Python-Math-Concepts-For-Developers
756639ff35bc81e1050f850dd30e76948da290d9
MIT
Qwen/Qwen-72B
1. YES
2. YES