Let's get started with Sympy, Sicipy installation checks
#verify that scipy is installed
!pip install scipy
#verify that sympy is installed
!pip install sympy
#verify that numpy is installed
!pip install numpy
#verify that matplotlib is installed
!pip install matplotlib
Requirement already satisfied: scipy in /srv/conda/lib/python3.7/site-packages (1.2.1) Requirement already satisfied: numpy>=1.8.2 in /srv/conda/lib/python3.7/site-packages (from scipy) (1.16.2) You are using pip version 18.1, however version 19.0.3 is available. You should consider upgrading via the 'pip install --upgrade pip' command. Requirement already satisfied: sympy in /srv/conda/lib/python3.7/site-packages (1.3) Requirement already satisfied: mpmath>=0.19 in /srv/conda/lib/python3.7/site-packages (from sympy) (1.1.0) You are using pip version 18.1, however version 19.0.3 is available. You should consider upgrading via the 'pip install --upgrade pip' command. Requirement already satisfied: numpy in /srv/conda/lib/python3.7/site-packages (1.16.2) You are using pip version 18.1, however version 19.0.3 is available. You should consider upgrading via the 'pip install --upgrade pip' command. Requirement already satisfied: matplotlib in /srv/conda/lib/python3.7/site-packages (2.1.1) Requirement already satisfied: numpy>=1.7.1 in /srv/conda/lib/python3.7/site-packages (from matplotlib) (1.16.2) Requirement already satisfied: six>=1.10 in /srv/conda/lib/python3.7/site-packages (from matplotlib) (1.12.0) Requirement already satisfied: python-dateutil>=2.0 in /srv/conda/lib/python3.7/site-packages (from matplotlib) (2.8.0) Requirement already satisfied: pytz in /srv/conda/lib/python3.7/site-packages (from matplotlib) (2018.9) Requirement already satisfied: cycler>=0.10 in /srv/conda/lib/python3.7/site-packages (from matplotlib) (0.10.0) Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /srv/conda/lib/python3.7/site-packages (from matplotlib) (2.3.1) You are using pip version 18.1, however version 19.0.3 is available. You should consider upgrading via the 'pip install --upgrade pip' command.
#verify that PrimeFinders is installed
!pip install git+https://github.com/LaGuer/PrimeFinders.git
from primefinders import primefinders
primefinders.about()
'primefinders.check(n)' returns 'True' if 'n' is a prime number 'primefinders.factor(n)' returns the lowest prime factor of 'n' 'primefinders.factors(n)' returns all the prime factors of 'n' with multiplicity 'primefinders.first(n)' returns first 'n' many primes 'primefinders.upto(n)' returns all the primes less than or equal to 'n' 'primefinders.between(m,n)' returns all the primes between 'm' and 'n' 'primefinders.phi(n)' returns the Euler's phi(n) 'primefinders.lcg(n)' example of Linear Congruential Generator 'primefinders.getResult(n)' returns 'True' if 'n' getResult 'primefinders.calculate(n)' returns 1 if 'n' is prime 'primefinders.fermat(n)' returns 1 if 'n' is prime using fermat a ** (n - 1)) % n) 'primefinders.check_complex(n)' returns 1 if 'n' might be prime it said to be complex
#9410309769636119433901981144786048550681811636889
primefinders.check(81811636889)
False
primefinders.factor(81811636889)
97
primefinders.factors(812156416)
[2, 2, 2, 2, 2, 2, 2, 2, 2, 239, 6637]
primefinders.first(20)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71]
primefinders.upto(10)
[2, 3, 5, 7]
primefinders.between(1690,1712)
[1693, 1697, 1699, 1709]
primefinders.calculate(16889)
Pool filled successfully **Progression** 25% 50% 75% 100%
1
primefinders.between(1700,1800)
[1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789]
primefinders.between(1800,1900)
[1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889]
primefinders.between(1900,2000)
[1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999]
primefinders.between(2000,2100)
[2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099]
primefinders.factor(13543)
29
primefinders.phi(1801)
1800
import scipy
from scipy import constants
scipy.e
2.718281828459045
scipy.euler_gamma
0.5772156649015329
scipy.constants.m_p / scipy.constants.m_e
1836.1526737600675
Electron Compton Wavelength
# Electron Compton Wavelength is given by :
scipy.constants.hbar / (scipy.constants.c * scipy.constants.m_e)
scipy.constants.hbar
1.0545718001391127e-34
3.8615926754860414*scipy.euler_gamma
2.2289717837595644
primefinders.factors(22289717837595644)
[2, 2, 4051, 1375568861861]
primefinders.factors(10545718001391127)
[752933, 14006183819]
primefinders.factors(3698998514839191553)
[11, 1280989, 262510132607]
import sympy
from sympy import exp
sympy.exp(841)
exp(841)
scipy.pi*scipy.e
8.539734222673566
# is this number prime ?
8539734222673563**841
primefinders.getResult()
Now let's switch to some sympy example
from sympy import *
import sympy.ntheory as nt
init_printing()
nt.isprime(1733)
True
nt.isprime(8539734222673563**841)
False
nt.isprime(453591899840915007729783809943324448515944275642519571972837233909606336899794963852468361313639693847640541257130320827382582057392503295025274655349406230018127758894477562917938182041375905178819902540327566610936581907458846282589340601330150248791041946045553400702783988196496244462024173043023555660702657265785984130911118593005675198407906455520015459251346288975829859494581215488934704174197022688253615839462685436311656563909471252403283650066079310697230199791238114794817420784786652045620848211906115555849470114651960030326983885443539241490043158682185243680733405300500655909231366745228810213719756031690000834101079241174001386432783889407786148292548723971772763295482054235296488645938171785380028649324452284511593929155508141904557038141954657993486810273488273116967649462158330683522940641483355408530829183104484735130712786959169)
True
nt.nextprime(9410309769636119433901981144786048550681811636753)
(9410309769636119433901981144786048550681811636753**4)
nt.nextprime(941000030976000009636119433000090190081144786048550681811000636753**13)
import numpy as np
import matplotlib.pyplot as plt
nt.primepi(1797)
import numpy as np
import matplotlib.pyplot as plt
#%matplotlib inline
x = np.arange(0, 60)
fig, ax = plt.subplots(1, 1, figsize=(10, 12))
ax.plot(x, np.log(2)**(x/4), '--k', label='$\log(2)^(\frac{x}{4}) $ ')
ax.plot(x, np.log(2)**(x/2), '-k', label='$\log(2)^(\frac{x}{2}) $ ')
ax.plot(x, np.log(2)**(x/np.pi), 'g', label='$\log(2)^(\frac{x}{2}) $ ')
ax.plot(x, np.log(2)**(x/np.e), 'c', label='$\log(2)^(\frac{x}{2}) $ ')
#ax.plot(x, list(map(np.exp(np.pi, x)**(x/4), x)), '-k', label='$pi(x)$')
#ax.plot(x, list(map(nt.primepi, x)), '-k', label='$pi(x)$')
#ax.plot(x, log(x)/4, '-k', label='$pi(x)$')
#ax.plot(x, x / np.log(x), '--k', label='$x/log(x)$')
#ax.plot(x, x / np.log(x), '--k', label='$x/log(x)$')
#ax.plot(x, np.exp(2)**(x/4), '--k', label='$exp(2)^(x/4)$')
#ax.plot(x, np.log(2)**(x/4), '--k', label='$log(2)^(x/4)$')
ax.legend(loc=2)
<matplotlib.legend.Legend at 0x7f3ae1ebdc18>
First off we need to check if all necessary libraries are installed
Write a function that tests if a number is prime. Test it by writing out all prime numbers less than 50.
This is a "simple" solution, but not efficient.
def isprime(n):
"""
Checks to see if an integer is prime.
Parameters
----------
n : integer
Number to check
Returns
-------
isprime : Boolean
If n is prime
"""
# No number less than 2 can be prime
if n < 2:
return False
# We only need to check for divisors up to sqrt(n)
for m in range(2, int(n**0.5)+1):
if n%m == 0:
return False
# If we've got this far, there are no divisors.
return True
for n in range(50):
if isprime(n):
print("Function says that {} is prime.".format(n))
Function says that 2 is prime. Function says that 3 is prime. Function says that 5 is prime. Function says that 7 is prime. Function says that 11 is prime. Function says that 13 is prime. Function says that 17 is prime. Function says that 19 is prime. Function says that 23 is prime. Function says that 29 is prime. Function says that 31 is prime. Function says that 37 is prime. Function says that 41 is prime. Function says that 43 is prime. Function says that 47 is prime.
500 years ago some believed that the number $2^n - 1$ was prime for all primes $n$. Use your function to find the first prime $n$ for which this is not true.
We could do this many ways. This "elegant" solution says:
n = 2
while (not isprime(n)) or (isprime(2**n-1)):
n += 1
print("The first n such that 2^n-1 is not prime is {}.".format(n))
The first n such that 2^n-1 is not prime is 11.
The Mersenne primes are those that have the form $2^n-1$, where $n$ is prime. Use your previous solutions to generate all the $n < 40$ that give Mersenne primes.
for n in range(2, 41):
if isprime(n) and isprime(2**n-1):
print("n={} is such that 2^n-1 is prime.".format(n))
n=2 is such that 2^n-1 is prime. n=3 is such that 2^n-1 is prime. n=5 is such that 2^n-1 is prime. n=7 is such that 2^n-1 is prime. n=13 is such that 2^n-1 is prime. n=17 is such that 2^n-1 is prime. n=19 is such that 2^n-1 is prime. n=31 is such that 2^n-1 is prime.
Write a function to compute all prime factors of an integer $n$, including their multiplicities. Test it by printing the prime factors (without multiplicities) of $n = 17, \dots, 20$ and the multiplicities (without factors) of $n = 48$.
One effective solution is to return a dictionary, where the keys are the factors and the values are the multiplicities.
This solution uses the trick of immediately dividing $n$ by any divisor: this means we never have to check the divisor for being prime.
def prime_factors(n):
"""
Generate all the prime factors of n.
Parameters
----------
n : integer
Number to be checked
Returns
-------
factors : dict
Prime factors (keys) and multiplicities (values)
"""
factors = {}
m = 2
while m <= n:
if n%m == 0:
factors[m] = 1
n //= m
while n%m == 0:
factors[m] += 1
n //= m
m += 1
return factors
for n in range(17, 21):
print("Prime factors of {} are {}.".format(n, prime_factors(n).keys()))
print("Multiplicities of prime factors of 48 are {}.".format(prime_factors(48).values()))
Prime factors of 17 are dict_keys([17]). Prime factors of 18 are dict_keys([2, 3]). Prime factors of 19 are dict_keys([19]). Prime factors of 20 are dict_keys([2, 5]). Multiplicities of prime factors of 48 are dict_values([4, 1]).
Here we will do it directly.
def divisors(n):
"""
Generate all integer divisors of n.
Parameters
----------
n : integer
Number to be checked
Returns
-------
divs : list
All integer divisors, including 1.
"""
divs = [1]
m = 2
while m <= n/2:
if n%m == 0:
divs.append(m)
m += 1
return divs
for n in range(16, 21):
print("The divisors of {} are {}.".format(n, divisors(n)))
The divisors of 16 are [1, 2, 4, 8]. The divisors of 17 are [1]. The divisors of 18 are [1, 2, 3, 6, 9]. The divisors of 19 are [1]. The divisors of 20 are [1, 2, 4, 5, 10].
A perfect number $n$ is one where the divisors sum to $n$. For example, 6 has divisors 1, 2, and 3, which sum to 6. Use your previous solution to find all perfect numbers $n < 10,000$ (there are only four!).
We can do this much more efficiently than the code below using packages such as numpy
, but this is a "bare python" solution.
def isperfect(n):
"""
Check if a number is perfect.
Parameters
----------
n : integer
Number to check
Returns
-------
isperfect : Boolean
Whether it is perfect or not.
"""
divs = divisors(n)
sum_divs = 0
for d in divs:
sum_divs += d
return n == sum_divs
for n in range(2,10000):
if (isperfect(n)):
factors = prime_factors(n)
print("{} is perfect.\n"
"Divisors are {}.\n"
"Prime factors {} (multiplicities {}).".format(
n, divisors(n), factors.keys(), factors.values()))
6 is perfect. Divisors are [1, 2, 3]. Prime factors dict_keys([2, 3]) (multiplicities dict_values([1, 1])). 28 is perfect. Divisors are [1, 2, 4, 7, 14]. Prime factors dict_keys([2, 7]) (multiplicities dict_values([2, 1])). 496 is perfect. Divisors are [1, 2, 4, 8, 16, 31, 62, 124, 248]. Prime factors dict_keys([2, 31]) (multiplicities dict_values([4, 1])). 8128 is perfect. Divisors are [1, 2, 4, 8, 16, 32, 64, 127, 254, 508, 1016, 2032, 4064]. Prime factors dict_keys([2, 127]) (multiplicities dict_values([6, 1])).
Using your previous functions, check that all perfect numbers $n < 10,000$ can be written as $2^{k-1} \times (2^k - 1)$, where $2^k-1$ is a Mersenne prime.
In fact we did this above already:
Investigate the timeit
function in python or IPython. Use this to measure how long your function takes to check that, if $k$ on the Mersenne list then $n = 2^{k-1} \times (2^k - 1)$ is a perfect number, using your functions. Stop increasing $k$ when the time takes too long!
You could waste considerable time on this, and on optimizing the functions above to work efficiently. It is not worth it, other than to show how rapidly the computation time can grow!
%timeit isperfect(2**(3-1)*(2**3-1))
2.49 µs ± 215 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit isperfect(2**(5-1)*(2**5-1))
35.3 µs ± 4 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit isperfect(2**(7-1)*(2**7-1))
658 µs ± 64.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit isperfect(2**(13-1)*(2**13-1))
2.72 s ± 290 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
It's worth thinking about the operation counts of the various functions implemented here. The implementations are inefficient, but even in the best case you see how the number of operations (and hence computing time required) rapidly increases.
Partly taken from Newman's book, p 120.
The logistic map builds a sequence of numbers $\{ x_n \}$ using the relation
\begin{equation} x_{n+1} = r x_n \left( 1 - x_n \right), \end{equation}where $0 \le x_0 \le 1$.
Write a program that calculates the first $N$ members of the sequence, given as input $x_0$ and $r$ (and, of course, $N$).
def logistic(x0, r, N = 1000):
sequence = [x0]
xn = x0
for n in range(N):
xnew = r*xn*(1.0-xn)
sequence.append(xnew)
xn = xnew
return sequence
Fix $x_0=0.5$. Calculate the first 2,000 members of the sequence for $r=1.5$ and $r=3.5$ Plot the last 100 members of the sequence in both cases.
What does this suggest about the long-term behaviour of the sequence?
import numpy
from matplotlib import pyplot
%matplotlib inline
x0 = 0.5
N = 2000
sequence1 = logistic(x0, 1.5, N)
sequence2 = logistic(x0, 3.5, N)
pyplot.plot(sequence1[-100:], 'b-', label = r'$r=1.5$')
pyplot.plot(sequence2[-100:], 'k-', label = r'$r=3.5$')
pyplot.xlabel(r'$n$')
pyplot.ylabel(r'$x$')
pyplot.show()
This suggests that, for $r=1.5$, the sequence has settled down to a fixed point. In the $r=3.5$ case it seems to be moving between four points repeatedly.
Fix $x_0 = 0.5$. For each value of $r$ between $1$ and $4$, in steps of $0.01$, calculate the first 2,000 members of the sequence. Plot the last 1,000 members of the sequence on a plot where the $x$-axis is the value of $r$ and the $y$-axis is the values in the sequence. Do not plot lines - just plot markers (e.g., use the 'k.'
plotting style).
import numpy
from matplotlib import pyplot
%matplotlib inline
# This is the "best" way of doing it, but numpy hasn't been introduced yet
# r_values = numpy.arange(1.0, 4.0, 0.01)
r_values = []
for i in range(302):
r_values.append(1.0 + 0.01 * i)
x0 = 0.5
N = 2000
for r in r_values:
sequence = logistic(x0, r, N)
pyplot.plot(r*numpy.ones_like(sequence[1000:]), sequence[1000:], 'k.')
pyplot.xlabel(r'$r$')
pyplot.ylabel(r'$x$')
pyplot.show()
For iterative maps such as the logistic map, one of three things can occur:
Using just your plot, or new plots from this data, work out approximate values of $r$ for which there is a transition from fixed points to limit cycles, from limit cycles of a given number of values to more values, and the transition to chaos.
The first transition is at $r \approx 3$, the next at $r \approx 3.45$, the next at $r \approx 3.55$. The transition to chaos appears to happen before $r=4$, but it's not obvious exactly where.
The Mandelbrot set is also generated from a sequence, $\{ z_n \}$, using the relation
\begin{equation} z_{n+1} = z_n^2 + c, \qquad z_0 = 0. \end{equation}The members of the sequence, and the constant $c$, are all complex. The point in the complex plane at $c$ is in the Mandelbrot set only if the $|z_n| < 2$ for all members of the sequence. In reality, checking the first 100 iterations is sufficient.
Note: the python notation for a complex number $x + \text{i} y$ is x + yj
: that is, j
is used to indicate $\sqrt{-1}$. If you know the values of x
and y
then x + yj
constructs a complex number; if they are stored in variables you can use complex(x, y)
.
Write a function that checks if the point $c$ is in the Mandelbrot set.
def in_Mandelbrot(c, n_iterations = 100):
z0 = 0.0 + 0j
in_set = True
n = 0
zn = z0
while in_set and (n < n_iterations):
n += 1
znew = zn**2 + c
in_set = abs(znew) < 2.0
zn = znew
return in_set
Check the points $c=0$ and $c=\pm 2 \pm 2 \text{i}$ and ensure they do what you expect. (What should you expect?)
c_values = [0.0, 2+2j, 2-2j, -2+2j, -2-2j]
for c in c_values:
print("Is {} in the Mandelbrot set? {}.".format(c, in_Mandelbrot(c)))
Is 0.0 in the Mandelbrot set? True. Is (2+2j) in the Mandelbrot set? False. Is (2-2j) in the Mandelbrot set? False. Is (-2+2j) in the Mandelbrot set? False. Is (-2-2j) in the Mandelbrot set? False.
Write a function that, given $N$
import numpy
def grid_Mandelbrot(N):
x = numpy.linspace(-2.0, 2.0, N)
X, Y = numpy.meshgrid(x, x)
C = X + 1j*Y
grid = numpy.zeros((N, N), int)
for nx in range(N):
for ny in range(N):
grid[nx, ny] = int(in_Mandelbrot(C[nx, ny]))
return grid
Using the function imshow
from matplotlib
, plot the resulting array for a $100 \times 100$ array to make sure you see the expected shape.
from matplotlib import pyplot
%matplotlib inline
pyplot.imshow(grid_Mandelbrot(100))
<matplotlib.image.AxesImage at 0x7fe773e42fd0>
Modify your functions so that, instead of returning whether a point is inside the set or not, it returns the logarithm of the number of iterations it takes. Plot the result using imshow
again.
from math import log
def log_Mandelbrot(c, n_iterations = 100):
z0 = 0.0 + 0j
in_set = True
n = 0
zn = z0
while in_set and (n < n_iterations):
n += 1
znew = zn**2 + c
in_set = abs(znew) < 2.0
zn = znew
return log(n)
def log_grid_Mandelbrot(N):
x = numpy.linspace(-2.0, 2.0, N)
X, Y = numpy.meshgrid(x, x)
C = X + 1j*Y
grid = numpy.zeros((N, N), int)
for nx in range(N):
for ny in range(N):
grid[nx, ny] = log_Mandelbrot(C[nx, ny])
return grid
from matplotlib import pyplot
%matplotlib inline
pyplot.imshow(log_grid_Mandelbrot(100))
<matplotlib.image.AxesImage at 0x7fe773de1128>
Try some higher resolution plots, and try plotting only a section to see the structure. Note this is not a good way to get high accuracy close up images!
This is a simple example:
pyplot.imshow(log_grid_Mandelbrot(1000)[600:800,400:600])
<matplotlib.image.AxesImage at 0x7fe773d80320>
An equivalence class is a relation that groups objects in a set into related subsets. For example, if we think of the integers modulo $7$, then $1$ is in the same equivalence class as $8$ (and $15$, and $22$, and so on), and $3$ is in the same equivalence class as $10$. We use the tilde $3 \sim 10$ to denote two objects within the same equivalence class.
Here, we are going to define the positive integers programmatically from equivalent sequences.
Define a python class Eqint
. This should be
__repr__
function) to be the integer length of the sequence;__eq__
function) so that two eqint
s are equal if their sequences have same length.class Eqint(object):
def __init__(self, sequence):
self.sequence = sequence
def __repr__(self):
return str(len(self.sequence))
def __eq__(self, other):
return len(self.sequence)==len(other.sequence)
Define a zero
object from the empty list, and three one
objects, from a single object list, tuple, and string. For example
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')
Check that none of the one
objects equal the zero object, but all equal the other one
objects. Print each object to check that the representation gives the integer length.
zero = Eqint([])
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')
print("Is zero equivalent to one? {}, {}, {}".format(zero == one_list,
zero == one_tuple,
zero == one_string))
print("Is one equivalent to one? {}, {}, {}.".format(one_list == one_tuple,
one_list == one_string,
one_tuple == one_string))
print(zero)
print(one_list)
print(one_tuple)
print(one_string)
Is zero equivalent to one? False, False, False Is one equivalent to one? True, True, True. 0 1 1 1
Redefine the class by including an __add__
method that combines the two sequences. That is, if a
and b
are Eqint
s then a+b
should return an Eqint
defined from combining a
and b
s sequences.
Adding two different types of sequences (eg, a list to a tuple) does not work, so it is better to either iterate over the sequences, or to convert to a uniform type before adding.
class Eqint(object):
def __init__(self, sequence):
self.sequence = sequence
def __repr__(self):
return str(len(self.sequence))
def __eq__(self, other):
return len(self.sequence)==len(other.sequence)
def __add__(a, b):
return Eqint(tuple(a.sequence) + tuple(b.sequence))
Check your addition function by adding together all your previous Eqint
objects (which will need re-defining, as the class has been redefined). Print the resulting object to check you get 3
, and also print its internal sequence.
zero = Eqint([])
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')
sum_eqint = zero + one_list + one_tuple + one_string
print("The sum is {}.".format(sum_eqint))
print("The internal sequence is {}.".format(sum_eqint.sequence))
The sum is 3. The internal sequence is (1, 1, '1').
We will sketch a construction of the positive integers from nothing.
positive_integers
.Eqint
called zero
from the empty list. Append it to positive_integers
.Eqint
called next_integer
from the Eqint
defined by a copy of positive_integers
(ie, use Eqint(list(positive_integers))
. Append it to positive_integers
.Use this procedure to define the Eqint
equivalent to $10$. Print it, and its internal sequence, to check.
positive_integers = []
zero = Eqint([])
positive_integers.append(zero)
N = 10
for n in range(1,N+1):
positive_integers.append(Eqint(list(positive_integers)))
print("The 'final' Eqint is {}".format(positive_integers[-1]))
print("Its sequence is {}".format(positive_integers[-1].sequence))
print("That is, it contains all Eqints with length less than 10.")
The 'final' Eqint is 10 Its sequence is [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] That is, it contains all Eqints with length less than 10.
Instead of working with floating point numbers, which are not "exact", we could work with the rational numbers $\mathbb{Q}$. A rational number $q \in \mathbb{Q}$ is defined by the numerator $n$ and denominator $d$ as $q = \frac{n}{d}$, where $n$ and $d$ are coprime (ie, have no common divisor other than $1$).
Find a python function that finds the greatest common divisor (gcd
) of two numbers. Use this to write a function normal_form
that takes a numerator and divisor and returns the coprime $n$ and $d$. Test this function on $q = \frac{3}{2}$, $q = \frac{15}{3}$, and $q = \frac{20}{42}$.
def normal_form(numerator, denominator):
from fractions import gcd
factor = gcd(numerator, denominator)
return numerator//factor, denominator//factor
print(normal_form(3, 2))
print(normal_form(15, 3))
print(normal_form(20, 42))
(3, 2) (5, 1) (10, 21)
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
Define a class Rational
that uses the normal_form
function to store the rational number in the appropriate form. Define a __repr__
function that prints a string that looks like $\frac{n}{d}$ (hint: use len(str(number))
to find the number of digits of an integer). Test it on the cases above.
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
q1 = Rational(3, 2)
print(q1)
q2 = Rational(15, 3)
print(q2)
q3 = Rational(20, 42)
print(q3)
3 - 2 5 10 -- 21
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
Overload the __add__
function so that you can add two rational numbers. Test it on $\frac{1}{2} + \frac{1}{3} + \frac{1}{6} = 1$.
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
print(Rational(1,2) + Rational(1,3) + Rational(1,6))
1
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
Overload the __mul__
function so that you can multiply two rational numbers. Test it on $\frac{1}{3} \times \frac{15}{2} \times \frac{2}{5} = 1$.
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
print(Rational(1,3)*Rational(15,2)*Rational(2,5))
1
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
Overload the __rmul__
function so that you can multiply a rational by an integer. Check that $\frac{1}{2} \times 2 = 1$ and $\frac{1}{2} + (-1) \times \frac{1}{2} = 0$. Also overload the __sub__
function (using previous functions!) so that you can subtract rational numbers and check that $\frac{1}{2} - \frac{1}{2} = 0$.
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __rmul__(self, other):
numerator = self.numerator * other
return Rational(numerator, self.denominator)
def __sub__(a, b):
return a + (-1)*b
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
half = Rational(1,2)
print(2*half)
print(half+(-1)*half)
print(half-half)
1 0 0
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
Overload the __float__
function so that float(q)
returns the floating point approximation to the rational number q
. Test this on $\frac{1}{2}, \frac{1}{3}$, and $\frac{1}{11}$.
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __rmul__(self, other):
numerator = self.numerator * other
return Rational(numerator, self.denominator)
def __sub__(a, b):
return a + (-1)*b
def __float__(a):
return float(a.numerator) / float(a.denominator)
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
print(float(Rational(1,2)))
print(float(Rational(1,3)))
print(float(Rational(1,11)))
0.5 0.3333333333333333 0.09090909090909091
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
Overload the __lt__
function to compare two rational numbers. Create a list of rational numbers where the denominator is $n = 2, \dots, 11$ and the numerator is the floored integer $n/2$, ie n//2
. Use the sorted
function on that list (which relies on the __lt__
function).
class Rational(object):
"""
A rational number.
"""
def __init__(self, numerator, denominator):
n, d = normal_form(numerator, denominator)
self.numerator = n
self.denominator = d
return None
def __add__(a, b):
numerator = a.numerator * b.denominator + b.numerator * a.denominator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __mul__(a, b):
numerator = a.numerator * b.numerator
denominator = a.denominator * b.denominator
return Rational(numerator, denominator)
def __rmul__(self, other):
numerator = self.numerator * other
return Rational(numerator, self.denominator)
def __sub__(a, b):
return a + (-1)*b
def __float__(a):
return float(a.numerator) / float(a.denominator)
def __lt__(a, b):
return a.numerator * b.denominator < a.denominator * b.numerator
def __repr__(self):
max_length = max(len(str(self.numerator)), len(str(self.denominator)))
if self.denominator == 1:
frac = str(self.numerator)
else:
numerator = '\n'+str(self.numerator)+'\n'
bar = max_length*'-'+'\n'
denominator = str(self.denominator)
frac = numerator+bar+denominator
return frac
q_list = [Rational(n//2, n) for n in range(2, 12)]
print(sorted(q_list))
[ 1 - 3, 2 - 5, 3 - 7, 4 - 9, 5 -- 11, 1 - 2, 1 - 2, 1 - 2, 1 - 2, 1 - 2]
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
The Wallis formula for $\pi$ is
\begin{equation} \pi = 2 \prod_{n=1}^{\infty} \frac{ (2 n)^2 }{(2 n - 1) (2 n + 1)}. \end{equation}We can define a partial product $\pi_N$ as
\begin{equation} \pi_N = 2 \prod_{n=1}^{N} \frac{ (2 n)^2 }{(2 n - 1) (2 n + 1)}, \end{equation}each of which are rational numbers.
Construct a list of the first 20 rational number approximations to $\pi$ and print them out. Print the sorted list to show that the approximations are always increasing. Then convert them to floating point numbers, construct a numpy
array, and subtract this array from $\pi$ to see how accurate they are.
def wallis_rational(N):
"""
The partial product approximation to pi using the first N terms of Wallis' formula.
Parameters
----------
N : int
Number of terms in product
Returns
-------
partial : Rational
A rational number approximation to pi
"""
partial = Rational(2,1)
for n in range(1, N+1):
partial = partial * Rational((2*n)**2, (2*n-1)*(2*n+1))
return partial
pi_list = [wallis_rational(n) for n in range(1, 21)]
print(pi_list)
print(sorted(pi_list))
[ 8 - 3, 128 --- 45, 512 --- 175, 32768 ----- 11025, 131072 ------ 43659, 2097152 ------- 693693, 8388608 ------- 2760615, 2147483648 ---------- 703956825, 8589934592 ---------- 2807136475, 137438953472 ------------ 44801898141, 549755813888 ------------ 178837328943, 35184372088832 -------------- 11425718238025, 140737488355328 --------------- 45635265151875, 2251799813685248 ---------------- 729232910488125, 9007199254740992 ---------------- 2913690606794775, 9223372036854775808 ------------------- 2980705490751054825, 36893488147419103232 -------------------- 11912508103174630875, 590295810358705651712 --------------------- 190453061649520333125, 2361183241434822606848 ---------------------- 761284675790187924375, 151115727451828646838272 ------------------------ 48691767863540419643025] [ 8 - 3, 128 --- 45, 512 --- 175, 32768 ----- 11025, 131072 ------ 43659, 2097152 ------- 693693, 8388608 ------- 2760615, 2147483648 ---------- 703956825, 8589934592 ---------- 2807136475, 137438953472 ------------ 44801898141, 549755813888 ------------ 178837328943, 35184372088832 -------------- 11425718238025, 140737488355328 --------------- 45635265151875, 2251799813685248 ---------------- 729232910488125, 9007199254740992 ---------------- 2913690606794775, 9223372036854775808 ------------------- 2980705490751054825, 36893488147419103232 -------------------- 11912508103174630875, 590295810358705651712 --------------------- 190453061649520333125, 2361183241434822606848 ---------------------- 761284675790187924375, 151115727451828646838272 ------------------------ 48691767863540419643025]
/srv/conda/lib/python3.7/site-packages/ipykernel_launcher.py:4: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead. after removing the cwd from sys.path.
import numpy
print(numpy.pi-numpy.array(list(map(float, pi_list))))
[0.47492599 0.29714821 0.21587837 0.16943846 0.1394167 0.11842246 0.10291902 0.09100266 0.08155811 0.07388885 0.06753749 0.06219131 0.05762923 0.05369058 0.05025576 0.04723393 0.04455483 0.0421633 0.04001539 0.03807569]
A candidate for the shortest mathematical paper ever shows the following result:
\begin{equation} 27^5 + 84^5 + 110^5 + 133^5 = 144^5. \end{equation}This is interesting as
This is a counterexample to a conjecture by Euler ... that at least $n$ $n$th powers are required to sum to an $n$th power, $n > 2$.
Using python, check the equation above is true.
lhs = 27**5 + 84**5 + 110**5 + 133**5
rhs = 144**5
print("Does the LHS {} equal the RHS {}? {}".format(lhs, rhs, lhs==rhs))
Does the LHS 61917364224 equal the RHS 61917364224? True
The more interesting statement in the paper is that
\begin{equation} 27^5 + 84^5 + 110^5 + 133^5 = 144^5. \end{equation}[is] the smallest instance in which four fifth powers sum to a fifth power.
Interpreting "the smallest instance" to mean the solution where the right hand side term (the largest integer) is the smallest, we want to use python to check this statement.
You may find the combinations
function from the itertools
package useful.
import numpy
import itertools
The combinations
function returns all the combinations (ignoring order) of r
elements from a given list. For example, take a list of length 6, [1, 2, 3, 4, 5, 6]
and compute all the combinations of length 4:
input_list = numpy.arange(1, 7)
combinations = list(itertools.combinations(input_list, 4))
print(combinations)
[(1, 2, 3, 4), (1, 2, 3, 5), (1, 2, 3, 6), (1, 2, 4, 5), (1, 2, 4, 6), (1, 2, 5, 6), (1, 3, 4, 5), (1, 3, 4, 6), (1, 3, 5, 6), (1, 4, 5, 6), (2, 3, 4, 5), (2, 3, 4, 6), (2, 3, 5, 6), (2, 4, 5, 6), (3, 4, 5, 6)]
We can already see that the number of terms to consider is large.
Note that we have used the list
function to explicitly get a list of the combinations. The combinations
function returns a generator, which can be used in a loop as if it were a list, without storing all elements of the list.
How fast does the number of combinations grow? The standard formula says that for a list of length $n$ there are
\begin{equation} \begin{pmatrix} n \\ k \end{pmatrix} = \frac{n!}{k! (n-k)!} \end{equation}combinations of length $k$. For $k=4$ as needed here we will have $n (n-1) (n-2) (n-3) / 24$ combinations. For $n=144$ we therefore have
n_combinations = 144*143*142*141/24
print("Number of combinations of 4 objects from 144 is {}".format(n_combinations))
Number of combinations of 4 objects from 144 is 17178876.0
Show, by getting python to compute the number of combinations $N = \begin{pmatrix} n \\ 4 \end{pmatrix}$ that $N$ grows roughly as $n^4$. To do this, plot the number of combinations and $n^4$ on a log-log scale. Restrict to $n \le 50$.
from matplotlib import pyplot
%matplotlib inline
n = numpy.arange(5, 51)
N = numpy.zeros_like(n)
for i, n_c in enumerate(n):
combinations = list(itertools.combinations(numpy.arange(1,n_c+1), 4))
N[i] = len(combinations)
pyplot.figure(figsize=(12,6))
pyplot.loglog(n, N, linestyle='None', marker='x', color='k', label='Combinations')
pyplot.loglog(n, n**4, color='b', label=r'$n^4$')
pyplot.xlabel(r'$n$')
pyplot.ylabel(r'$N$')
pyplot.legend(loc='upper left')
pyplot.show()
With 17 million combinations to work with, we'll need to be a little careful how we compute.
One thing we could try is to loop through each possible "smallest instance" (the term on the right hand side) in increasing order. We then check all possible combinations of left hand sides.
This is computationally very expensive as we repeat a lot of calculations. We repeatedly recalculate combinations (a bad idea). We repeatedly recalculate the powers of the same number.
Instead, let us try creating the list of all combinations of powers once.
numpy
array containing all integers in $1, \dots, 144$ to the fifth power.in
keyword).import numpy as np
import itertools
nmax=145
range_to_power = np.arange(1, nmax)**5.0
lhs_combinations = list(itertools.combinations(range_to_power, 4))
Then calculate the sums:
lhs_sums = []
for lhs_terms in lhs_combinations:
lhs_sums.append(np.sum(np.array(lhs_terms)))
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-6-5c877a12f612> in <module> 1 lhs_sums = [] 2 for lhs_terms in lhs_combinations: ----> 3 lhs_sums.append(np.sum(np.array(lhs_terms))) /srv/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py in sum(a, axis, dtype, out, keepdims, initial) 2074 2075 return _wrapreduction(a, np.add, 'sum', axis, dtype, out, keepdims=keepdims, -> 2076 initial=initial) 2077 2078 /srv/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs) 84 return reduction(axis=axis, out=out, **passkwargs) 85 ---> 86 return ufunc.reduce(obj, axis, dtype, out, **passkwargs) 87 88 KeyboardInterrupt:
Finally, loop through the sums and check to see if it matches any possible term on the RHS:
for i, lhs in enumerate(lhs_sums):
if lhs in range_to_power:
rhs_primitive = int(lhs**(0.2))
lhs_primitive = (numpy.array(lhs_combinations[i])**(0.2)).astype(int)
print("The LHS terms are {}.".format(lhs_primitive))
print("The RHS term is {}.".format(rhs_primitive))
The Lorenz system is a set of ordinary differential equations which can be written
\begin{equation} \frac{\text{d} \vec{v}}{\text{d} \vec{t}} = \vec{f}(\vec{v}) \end{equation}where the variables in the state vector $\vec{v}$ are
\begin{equation} \vec{v} = \begin{pmatrix} x(t) \\ y(t) \\ z(t) \end{pmatrix} \end{equation}and the function defining the ODE is
\begin{equation} \vec{f} = \begin{pmatrix} \sigma \left( y(t) - x(t) \right) \\ x(t) \left( \rho - z(t) \right) - y(t) \\ x(t) y(t) - \beta z(t) \end{pmatrix}. \end{equation}The parameters $\sigma, \rho, \beta$ are all real numbers.
Write a function dvdt(v, t, params)
that returns $\vec{f}$ given $\vec{v}, t$ and the parameters $\sigma, \rho, \beta$.
def dvdt(v, t, sigma, rho, beta):
"""
Define the Lorenz system.
Parameters
----------
v : list
State vector
t : float
Time
sigma : float
Parameter
rho : float
Parameter
beta : float
Parameter
Returns
-------
dvdt : list
RHS defining the Lorenz system
"""
x, y, z = v
return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]
Fix $\sigma=10, \beta=8/3$. Set initial data to be $\vec{v}(0) = \vec{1}$. Using scipy
, specifically the odeint
function of scipy.integrate
, solve the Lorenz system up to $t=100$ for $\rho=13, 14, 15$ and $28$.
Plot your results in 3d, plotting $x, y, z$.
import numpy
from scipy.integrate import odeint
v0 = [1.0, 1.0, 1.0]
sigma = 10.0
beta = 8.0/3.0
t_values = numpy.linspace(0.0, 100.0, 5000)
rho_values = [13.0, 14.0, 15.0, 28.0]
v_values = []
for rho in rho_values:
params = (sigma, rho, beta)
v = odeint(dvdt, v0, t_values, args=params)
v_values.append(v)
%matplotlib inline
from matplotlib import pyplot
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = pyplot.figure(figsize=(12,6))
for i, v in enumerate(v_values):
ax = fig.add_subplot(2,2,i+1,projection='3d')
ax.plot(v[:,0], v[:,1], v[:,2])
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
ax.set_title(r"$\rho={}$".format(rho_values[i]))
pyplot.show()
Fix $\rho = 28$. Solve the Lorenz system twice, up to $t=40$, using the two different initial conditions $\vec{v}(0) = \vec{1}$ and $\vec{v}(0) = \vec{1} + \vec{10^{-5}}$.
Show four plots. Each plot should show the two solutions on the same axes, plotting $x, y$ and $z$. Each plot should show $10$ units of time, ie the first shows $t \in [0, 10]$, the second shows $t \in [10, 20]$, and so on.
t_values = numpy.linspace(0.0, 40.0, 4000)
rho = 28.0
params = (sigma, rho, beta)
v_values = []
v0_values = [[1.0,1.0,1.0],
[1.0+1e-5,1.0+1e-5,1.0+1e-5]]
for v0 in v0_values:
v = odeint(dvdt, v0, t_values, args=params)
v_values.append(v)
fig = pyplot.figure(figsize=(12,6))
line_colours = 'by'
for tstart in range(4):
ax = fig.add_subplot(2,2,tstart+1,projection='3d')
for i, v in enumerate(v_values):
ax.plot(v[tstart*1000:(tstart+1)*1000,0],
v[tstart*1000:(tstart+1)*1000,1],
v[tstart*1000:(tstart+1)*1000,2],
color=line_colours[i])
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
ax.set_title(r"$t \in [{},{}]$".format(tstart*10, (tstart+1)*10))
pyplot.show()
This shows the sensitive dependence on initial conditions that is characteristic of chaotic behaviour.
A twin prime is a pair $(p_1, p_2)$ such that both $p_1$ and $p_2$ are prime and $p_2 = p_1 + 2$.
Write a generator that returns twin primes. You can use the generators above, and may want to look at the itertools module together with its recipes, particularly the pairwise
recipe.
Note: we need to first pull in the generators introduced in that notebook
def all_primes(N):
"""
Return all primes less than or equal to N.
Parameters
----------
N : int
Maximum number
Returns
-------
prime : generator
Prime numbers
"""
primes = []
for n in range(2, N+1):
is_n_prime = True
for p in primes:
if n%p == 0:
is_n_prime = False
break
if is_n_prime:
primes.append(n)
yield n
Now we can generate pairs using the pairwise recipe:
from itertools import tee
def pair_primes(N):
"Generate consecutive prime pairs, using the itertools recipe"
a, b = tee(all_primes(N))
next(b, None)
return zip(a, b)
We could examine the results of the two primes directly. But an efficient solution is to use python's filter function. To do this, first define a function checking if the pair are twin primes:
def check_twin(pair):
"""
Take in a pair of integers, check if they differ by 2.
"""
p1, p2 = pair
return p2-p1 == 2
Then use the filter
function to define another generator:
def twin_primes(N):
"""
Return all twin primes
"""
return filter(check_twin, pair_primes(N))
Now check by finding the twin primes with $N<20$:
for tp in twin_primes(20):
print(tp)
(3, 5) (5, 7) (11, 13) (17, 19)
Find how many twin primes there are with $p_2 < 1000$.
Again there are many solutions, but the itertools recipes has the quantify
pattern. Looking ahead to exercise 3 we'll define:
def pi_N(N):
"""
Use the quantify pattern from itertools to count the number of twin primes.
"""
return sum(map(check_twin, pair_primes(N)))
pi_N(1000)
Let $\pi_N$ be the number of twin primes such that $p_2 < N$. Plot how $\pi_N / N$ varies with $N$ for $N=2^k$ and $k = 4, 5, \dots 16$. (You should use a logarithmic scale where appropriate!)
We've now done all the hard work and can use the solutions above.
import numpy
from matplotlib import pyplot
%matplotlib inline
N = numpy.array([2**k for k in range(4, 17)])
twin_prime_fraction = numpy.array(list(map(pi_N, N))) / N
pyplot.semilogx(N, twin_prime_fraction)
pyplot.xlabel(r"$N$")
pyplot.ylabel(r"$\pi_N / N$")
pyplot.show()
For those that have checked Wikipedia, you'll see Brun's theorem which suggests a specific scaling, that $\pi_N$ is bounded by $C N / \log(N)^2$. Checking this numerically on this data:
pyplot.semilogx(N, twin_prime_fraction * numpy.log(N)**2)
pyplot.xlabel(r"$N$")
pyplot.ylabel(r"$\pi_N \times \log(N)^2 / N$")
pyplot.show()