## Computing Sums

Suppose we want to compute the sum of a sequence of numbers $x_0$, $x_1$, $x_2$, $x_3$, $\dots$, $x_n$. There are at least two approaches:

1. Compute the entire sequence, store it as a list $[x_0,x_1,x_2,\dots,x_n]$ and then use the built-in function `sum`.
2. Initialize a variable with value 0 (and name it `result` for example), create and add each element in the sequence to `result` one at a time.

The advantage of the second approach is that we don't need to store all the values at once. For example, here are two ways to write a function which computes the sum of squares.

For the first approach, use a list comprehension:


```python
def sum_of_squares_1(N):
    "Compute the sum of squares 1**2 + 2**2 + ... + N**2."
    return sum([n**2 for n in range(1,N + 1)])
```


```python
sum_of_squares_1(4)
```




    30



For the second approach, use a `for` loop with the initialize-and-update construction:


```python
def sum_of_squares_2(N):
    "Compute the sum of squares 1**2 + 2**2 + ... + N**2."
    # Initialize the output value to 0
    result = 0
    for n in range(1,N + 1):
        # Update the result by adding the next term
        result = result + n**2
    return result
```


```python
sum_of_squares_2(4)
```




    30



Again, both methods yield the same result however the second uses less memory!

## Computing Products

There is no built-in function to compute products of sequences therefore we'll use an initialize-and-update construction similar to the example above for computing sums.

Write a function called `factorial` which takes a positive integer $N$ and return the factorial $N!$.


```python
def factorial(N):
    "Compute N! = N(N-1) ... (2)(1) for N >= 1."
    # Initialize the output variable to 1
    product = 1
    for n in range(2,N + 1):
        # Update the output variable
        product = product * n
    return product
```

Let's test our function for input values for which we know the result:


```python
factorial(2)
```




    2




```python
factorial(5)
```




    120



We can use our function to approximate $e$ using the Taylor series for $e^x$:

$$
e^x = \sum_{k=0}^{\infty} \frac{x^k}{k!}
$$

For example, let's compute the 100th partial sum of the series with $x=1$:


```python
sum([1/factorial(k) for k in range(0,101)])
```




    2.7182818284590455



## Searching for Solutions

We can use `for` loops to search for integer solutions of equations. For example, suppose we would like to find all representations of a positive integer $N$ as a [sum of two squares](https://en.wikipedia.org/wiki/Sum_of_two_squares_theorem). In other words, we want to find all integer solutions $(x,y)$ of the equation:

$$
x^2 + y^2 = N
$$

Write a function called `reps_sum_squares` which takes an integer $N$ and finds all representations of $N$ as a sum of squares $x^2 + y^2 = N$ for $0 \leq x \leq y$. The function returns the representations as a list of tuples. For example, if $N = 50$ then $1^2 + 7^2 = 50$ and $5^2 + 5^2 = 50$ and the function returns the list `[(1, 7),(5, 5)]`.

Let's outline our approach before we write any code:

1. Given $x \leq y$, the largest possible value for $x$ is $\sqrt{\frac{N}{2}}$
2. For $x \leq \sqrt{\frac{N}{2}}$, the pair $(x,y)$ is a solution if $N - x^2$ is a square
3. Define a helper function called `is_square` to test if an integer is square


```python
def is_square(n):
    "Determine if the integer n is a square."
    if round(n**0.5)**2 == n:
        return True
    else:
        return False

def reps_sum_squares(N):
    '''Find all representations of N as a sum of squares x**2 + y**2 = N.

    Parameters
    ----------
    N : integer

    Returns
    -------
    reps : list of tuples of integers
        List of tuples (x,y) of positive integers such that x**2 + y**2 = N.

    Examples
    --------
    >>> reps_sum_squares(1105)
    [(4, 33), (9, 32), (12, 31), (23, 24)]
    '''
    reps = []
    if is_square(N/2):
        # If N/2 is a square, search up to x = (N/2)**0.5
        max_x = round((N/2)**0.5)
    else:
        # If N/2 is not a square, search up to x = floor((N/2)**0.5)
        max_x = int((N/2)**0.5)
    for x in range(0,max_x + 1):
        y_squared = N - x**2
        if is_square(y_squared):
            y = round(y_squared**0.5)
            # Append solution (x,y) to list of solutions
            reps.append((x,y))
    return reps
```


```python
reps_sum_squares(1105)
```




    [(4, 33), (9, 32), (12, 31), (23, 24)]



What is the smallest integer which can be expressed as the sum of squares in 5 different ways?


```python
N = 1105
num_reps = 4
while num_reps < 5:
    N = N + 1
    reps = reps_sum_squares(N)
    num_reps = len(reps)
print(N,':',reps_sum_squares(N))
```

    4225 : [(0, 65), (16, 63), (25, 60), (33, 56), (39, 52)]


## Examples

### Prime Numbers

A positive integer is [prime](https://en.wikipedia.org/wiki/Prime_number) if it is divisible only by 1 and itself. Write a function called `is_prime` which takes an input parameter `n` and returns `True` or `False` depending on whether `n` is prime or not.

Let's outline our approach before we write any code:

1. An integer $d$ divides $n$ if there is no remainder of $n$ divided by $d$.
2. Use the modulus operator `%` to compute the remainder.
3. If $d$ divides $n$ then $n = d q$ for some integer $q$ and either $d \leq \sqrt{n}$ or $q \leq \sqrt{n}$ (and not both), therefore we need only test if $d$ divides $n$ for integers $d \leq \sqrt{n}$


```python
def is_prime(n):
    "Determine whether or not n is a prime number."
    if n <= 1:
        return False
    # Test if d divides n for d <= n**0.5
    for d in range(2,round(n**0.5) + 1):
        if n % d == 0:
            # n is divisible by d and so n is not prime
            return False
    # If we exit the for loop, then n is not divisible by any d
    # and therefore n is prime
    return True
```

Let's test our function on the first 30 numbers:


```python
for n in range(0,31):
    if is_prime(n):
        print(n,'is prime!')
```

    2 is prime!
    3 is prime!
    5 is prime!
    7 is prime!
    11 is prime!
    13 is prime!
    17 is prime!
    19 is prime!
    23 is prime!
    29 is prime!


Our function works! Let's find all the primes between 20,000 and 20,100.


```python
for n in range(20000,20100):
    if is_prime(n):
        print(n,'is prime!')
```

    20011 is prime!
    20021 is prime!
    20023 is prime!
    20029 is prime!
    20047 is prime!
    20051 is prime!
    20063 is prime!
    20071 is prime!
    20089 is prime!


### Divisors

Let's write a function called `divisors` which takes a positive integer $N$ and returns the list of positive integers which divide $N$.


```python
def divisors(N):
    "Return the list of divisors of N."
    # Initialize the list of divisors (which always includes 1)
    divisor_list = [1]
    # Check division by d for d <= N/2
    for d in range(2,N // 2 + 1):
        if N % d == 0:
            divisor_list.append(d)
    # N divides itself and so we append N to the list of divisors
    divisor_list.append(N)
    return divisor_list
```

Let's test our function:


```python
divisors(10)
```




    [1, 2, 5, 10]




```python
divisors(100)
```




    [1, 2, 4, 5, 10, 20, 25, 50, 100]




```python
divisors(59)
```




    [1, 59]



### Collatz Conjecture

Let $a$ be a positive integer and consider the recursive sequence where $x_0 = a$ and

$$
x_{n+1} = \left\\{ \begin{array}{cl} x_n/2 & \text{if } x_n \text{ is even} \\\\ 3x_n+1 & \text{if } x_n \text{ is odd}  \end{array} \\right.
$$

The [Collatz conjecture](https://en.wikipedia.org/wiki/Collatz_conjecture) states that this sequence will *always* reach 1. For example, if $a = 10$ then $x_0 = 10$, $x_1 = 5$, $x_2 = 16$, $x_3 = 8$, $x_4 = 4$, $x_5 = 2$ and $x_6 = 1$.

Write a function called `collatz` which takes one input parameter `a` and returns the sequence of integers defined above and ending with the first occurrence $x_n=1$.


```python
def collatz(a):
    "Compute the Collatz sequence starting at a and ending at 1."
    # Initialize list with first value a
    sequence = [a]
    # Compute values until we reach 1
    while sequence[-1] > 1:
        # Check if the last element in the list is even
        if sequence[-1] % 2 == 0:
            # Compute and append the new value
            sequence.append(sequence[-1] // 2)
        else:
            # Compute and append the new value
            sequence.append(3*sequence[-1] + 1)
    return sequence
```

Let's test our function:


```python
print(collatz(10))
```

    [10, 5, 16, 8, 4, 2, 1]



```python
collatz(22)
```




    [22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1]



The Collatz conjecture is quite amazing. No matter where we start, the sequence always terminates at 1!


```python
a = 123456789
seq = collatz(a)
print("Collatz sequence for a =",a)
print("begins with",seq[:5])
print("ends with",seq[-5:])
print("and has",len(seq),"terms.")
```

    Collatz sequence for a = 123456789
    begins with [123456789, 370370368, 185185184, 92592592, 46296296]
    ends with [16, 8, 4, 2, 1]
    and has 178 terms.


Which $a < 1000$ produces the longest sequence?


```python
max_length = 1
a_max = 1
for a in range(1,1001):
    seq_length = len(collatz(a))
    if seq_length > max_length:
        max_length = seq_length
        a_max = a
print('Longest sequence begins with a =',a_max,'and has length',max_length)
```

    Longest sequence begins with a = 871 and has length 179


## Exercises

1. [Fermat's theorem on the sum of two squares](https://en.wikipedia.org/wiki/Fermat%27s_theorem_on_sums_of_two_squares) states that every prime number $p$ of the form $4k+1$ can be expressed as the sum of two squares. For example, $5 = 2^2 + 1^2$ and $13 = 3^2 + 2^2$. Find the smallest prime greater than $2019$ of the form $4k+1$ and write it as a sum of squares. (Hint: Use the functions `is_prime` and `reps_sum_squares` from this section.)

2. What is the smallest prime number which can be represented as a sum of squares in 2 different ways?

3. What is the smallest integer which can be represented as a sum of squares in 3 different ways?

4. Write a function called `primes_between` which takes two integer inputs $a$ and $b$ and returns the list of primes in the closed interval $[a,b]$.

5. Write a function called `primes_d_mod_N` which takes four integer inputs $a$, $b$, $d$ and $N$ and returns the list of primes in the closed interval $[a,b]$ which are congruent to $d$ mod $N$ (this means that the prime has remainder $d$ after division by $N$). This kind of list is called [primes in an arithmetic progression](https://en.wikipedia.org/wiki/Dirichlet%27s_theorem_on_arithmetic_progressions).

6. Write a function called `reciprocal_recursion` which takes three positive integers $x_0$, $x_1$ and $N$ and returns the sequence $[x_0,x_1,x_2,\dots,x_N]$ where

    $$
    x_n = \frac{1}{x_{n-1}} + \frac{1}{x_{n-2}}
    $$

7. Write a function called `root_sequence` which takes input parameters $a$ and $N$, both positive integers, and returns the $N$th term $x_N$ in the sequence:

    $$
    \begin{align}
    x_0 &= a \\\
    x_n &= 1 + \sqrt{x_{n-1}}
    \end{align}
    $$

    Does the sequence converge to different values for different starting values $a$?

8. Write a function called `fib_less_than` which takes one input $N$ and returns the list of Fibonacci numbers less than $N$.

9. Write a function called `fibonacci_primes` which takes an input parameter $N$ and returns the list of Fibonacci numbers less than $N$ which are also prime numbers.

10. Let $w(N)$ be the number of ways $N$ can be expressed as a sum of two squares $x^2 + y^2 = N$ with $1 \leq x \leq y$. Then

    $$
    \lim_{N \to \infty} \frac{1}{N} \sum_{n=1}^{N} w(n) = \frac{\pi}{8}
    $$

    Compute the left side of the formula for $N=100$ and compare the result to $\pi / 8$.

11. A list of positive integers $[a,b,c]$ (with $1 \leq a < b$) are a [Pythagorean triple](https://en.wikipedia.org/wiki/Pythagorean_triple) if $a^2 + b^2 = c^2$. Write a function called `py_triples` which takes an input parameter $N$ and returns the list of Pythagorean triples `[a,b,c]` with $c \leq N$.

a0a2c22932ad51d3467ba4d3ba92feea9eca94fa
ipynb
Jupyter Notebook
Python/1. Python Basics/Notebooks/1. The Basics.- Startups/1.0 Fundamentals Continued/1.4b Using Functions to Solve Math Problems.ipynb
okara83/Becoming-a-Data-Scientist
f09a15f7f239b96b77a2f080c403b2f3e95c9650
MIT
Python/1. Python Basics/Notebooks/1. The Basics.- Startups/1.0 Fundamentals Continued/1.4b Using Functions to Solve Math Problems.ipynb
okara83/Becoming-a-Data-Scientist
f09a15f7f239b96b77a2f080c403b2f3e95c9650
MIT
Python/1. Python Basics/Notebooks/1. The Basics.- Startups/1.0 Fundamentals Continued/1.4b Using Functions to Solve Math Problems.ipynb
okara83/Becoming-a-Data-Scientist
f09a15f7f239b96b77a2f080c403b2f3e95c9650
MIT
2022-02-09T15:41:33.000Z
2022-02-11T07:47:40.000Z
Qwen/Qwen-72B
1. YES
2. YES