# COURSE: Master math by coding in Python
## SECTION: Number theory

#### https://www.udemy.com/course/math-with-python/?couponCode=MXC-DISC4ALL
#### INSTRUCTOR: sincxpress.com

Note about this code: Each video in this section of the course corresponds to a section of code below. Please note that this code roughly matches the code shown in the live recording, but is not exactly the same -- the variable names, order of lines, and parameters may be slightly different. 


```python
# import all necessary modules
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
```

# VIDEO: Counting perfect numbers


```python
# From https://en.wikipedia.org/wiki/Perfect_number:
# "A perfect number is a positive integer that is equal to the sum of its positive divisors, excluding itself."


# create a function that tests whether the input is a perfect number
def test4perfect(number):

    # initialize
    divisors = []

    # find integer divisors
    for d in range(1,number):
        if number%d == 0:
            divisors.append(d)
  
    # test for equality
    return sum(divisors)==number
```


```python
test4perfect(30)
```


```python
perfect = []

for number in range(2,10001):
    if test4perfect(number):
        perfect.append(number)

print(perfect)
```

### Exercise


```python
# implement the function as a list-comprehension

def test4perfect(number):
    return number == sum([d for d in range (1,number) if number%d==0])
```


```python
## list comprehension intro

# as a loop
lst = []
for i in range(10):
    lst.append(i**2)
    
# as list comprehension
lst2 = [i**2 for i in range(10)]

print(lst)
print(lst2)
```


```python

```

# VIDEO: Euclid's Pythagorean triplets

$a=m^{2}-n^{2}\\
b=2mn\\
c=m^{2}+n^{2}\\
m>n$




```python
# specify range of numbers
numberrange = np.arange(2,51)

# initialize empty lists
a = []
b = []
c = []

for m in numberrange:
    for n in numberrange:
        a.append( m**2 - n**2 )
        b.append( 2*m*n )
        c.append( m**2 + n**2 )

```


```python
fig,ax = plt.subplots(1,figsize=(10,7))
plt.scatter(a,b,c=c, marker='o',alpha=.5,linewidths=.5,edgecolors='k')
plt.axis('off')
plt.show()
```


```python
check = np.array(a)**2 + np.array(b)**2 - np.array(c)**2

plt.plot(check,'o')
plt.ylabel('$y \quad = \quad a^2+b^2-c^2$')
plt.show()
```


```python

```

### Exercise


```python
# the code written above does not strictly conform to the equations. Find what's wrong and fix it!
```


```python
# initialize empty lists
a = []
b = []
c = []

for m in numberrange:
    for n in range(numberrange[0],m):
        a.append( m**2 - n**2 )
        b.append( 2*m*n )
        c.append( m**2 + n**2 )

fig,ax = plt.subplots(1,figsize=(7,7))
plt.scatter(a,b,c=c, marker='o',alpha=.5,linewidths=.5,edgecolors='k')
plt.axis('off')
plt.show()
```


```python

```


```python

```

# VIDEO: Fermat's theorem

There are no integer solutions to $x^3+y^3=z^3$


```python
ints = np.arange(1,101)

z = np.zeros((len(ints),len(ints)))

# run the simulation
for x in ints:
    for y in ints:
        z[x-1,y-1] = (x**3 + y**3)**(1/3)

```


```python
# check for integers
zInts = z%1 == 0

```


```python

# now visualize
fig,ax = plt.subplots(1,2,figsize=(12,6))
fig.tight_layout(pad=5)

# show the resulting z-values
h = ax[0].imshow(z)
fig.colorbar(h,ax=ax[0],fraction=.045)
ax[0].set_xlabel('y')
ax[0].set_ylabel('x')
ax[0].set_title('$z=\sqrt[3]{x^3+y^3}$',fontsize=16)

# show the boolean integer map
h = ax[1].imshow(zInts,vmin=0,vmax=1)
fig.colorbar(h,ax=ax[1],fraction=.045)
ax[1].set_xlabel('y')
ax[1].set_ylabel('x')
ax[1].set_title('Whether $z$ is int',fontsize=16)

plt.show()
```

### Exercise: "Fermat's stairway to heaven"


```python
# Modify the code to create a map of integer Pythagorean triplets
```


```python

```

# VIDEO: Plotting number sequences



```python
# initialize
s1 = np.array([])
s2 = np.array([])
s3 = np.array([])
s4 = np.array([])
xx = np.arange(1,81)


# generate the sequences
for n in xx:
    s1 = np.append(s1,n/(n+1))
    s2 = np.append(s2,(n+1)/n)
    s3 = np.append(s3,n/sum(np.arange(1,n+1)))
    s4 = np.append(s4,(-1)**n/np.sqrt(n))


# plotting
fig,ax = plt.subplots(1,figsize=(10,7))
plt.plot(xx,s1,label='$a/(a+1)$')
plt.plot(xx,s2,label='$(a+1)/a$')
plt.plot(xx,s3,label='$a/\Sigma(a)$')
plt.plot(xx,s4,label='$-1^a/\sqrt{a}$')
plt.legend(loc='upper right',fontsize=16)
plt.show()
```


```python

```

### Exercise


```python
# The series above are all convergent. Create a divergent series.
```


```python
# repeat without a for-loop!

# the a's
aSeq = np.arange(1,N)

# the sequences
s1 = aSeq/(aSeq+1)
s2 = (aSeq+1)/aSeq
s3 = aSeq/np.cumsum(aSeq)
s4 = (-1)**aSeq/np.sqrt(aSeq)

# examples of divergent sequences
# s3 = np.cumsum(aSeq)/aSeq
# s4 = np.sqrt(aSeq)/(-1)**aSeq

# plotting
fig,ax = plt.subplots(1,figsize=(10,7))
plt.plot(xx,s1,label='$a/(a+1)$')
plt.plot(xx,s2,label='$(a+1)/a$')
plt.plot(xx,s3,label='$a/\Sigma(a)$')
plt.plot(xx,s4,label='$-1^a/\sqrt{a}$')
plt.legend(loc='upper right',fontsize=16)
plt.show()
```


```python

```

# VIDEO: Heron's method of square roots


```python
# Let's test the algorithm in a simple example
num = 100

# initialze the first guess
x = num/3

# now run through the algorithm
for n in range(5):
    x = ( x + num/x )/2
    print(x)

x
```


```python
# define the numbers to compute
nums2sqrt = np.linspace(2,101,50)

# range of algorithm iterations
niterations = np.arange(3,9)


# initialize matrix of estimates
err = np.zeros((len(niterations),len(nums2sqrt)))


# loop over the numbers to compute
for ni,num in enumerate(nums2sqrt):
  
    # loop over number of iterations
    for ii,iters in enumerate(niterations):
  
        # initial guess
        x = num/3
      
        # implement algorithm
        for n in range(iters):
            x = ( x + num/x )/2
      
        # store error magnitude
        err[ii,ni] = abs(x-np.sqrt(num))

```


```python
plt.imshow(-np.log(err),aspect=10,extent=[nums2sqrt[0],nums2sqrt[-1],niterations[-1],niterations[0]])
plt.xlabel('Number')
plt.ylabel('Iterations')
plt.title('-log(error)')
plt.show()
```

### Exercise


```python
# Exercise: Heron's mosquito space ship #13
# select target number to compute square root of
targnum = 13

# range of starting values
starting = np.linspace(-targnum,targnum,40)

# number of iterations (fixed)
niters = 5

# initialize output matrix
sqrtAlgResults = np.zeros((len(starting),niters+1))


# loop over starting numbers
for idx,startnum in enumerate(starting):
    
    # initialize starting number as a list!
    x = [startnum]
    
    # now run through the algorithm
    for n in range(niters):
        betterguess = ( x[n] + targnum/x[n] )/2
        x.append( betterguess )
    
    sqrtAlgResults[idx,:] = x


# finally, plot the results
fig,ax = plt.subplots(1,figsize=(8,6))
plt.plot(np.arange(niters+1),sqrtAlgResults.T,linewidth=2)
plt.xlabel('Algorithm iteration')
plt.ylabel('Numerical approximation')
plt.title('Estimating $\sqrt{%g}$ = %g' %(targnum,sqrtAlgResults[-1,-1]))
plt.show()
```


```python

```


```python

```

# VIDEO: Smooth numbers





```python
# largest number to search in
maxN = 10000

# find the largest prime factor
largestPrimeFact = np.zeros(maxN+1,dtype=int)/0

for i in range(2,maxN+1): # note: start at 2! not zero!
    largestPrimeFact[i] = np.max(sym.primefactors(i))

# show the smooth numbers on a plot
fig,ax = plt.subplots(1,figsize=(10,7))
plt.plot(largestPrimeFact,'k.',markersize=.5)
plt.xlabel('Number')
plt.ylabel('Largest prime factor')
plt.show()
```


```python
# count the number of occurrences of maximum prime factors

# find all unique max primes (basically all the primes)
uniqueMaxPrimes = np.unique(largestPrimeFact)

# count the number of times each max-prime appears
num = np.zeros(len(uniqueMaxPrimes),dtype=int)
for i,u in enumerate(uniqueMaxPrimes):
    num[i] = np.sum(largestPrimeFact==u)

# and visualize the counts by the smooth numbers
plt.plot(uniqueMaxPrimes,num,'k.')
plt.xlabel('Max prime number')
plt.ylabel('Number of occurrences')
plt.xlim([-1,500])
plt.title('%s is the modal max prime factor' %uniqueMaxPrimes[np.argmax(num)])
plt.show()
```

### Exercise


```python
# list all of the 5-smooth numbers up to maxN
numberlist = np.arange(maxN+1)

print('All 5-smooth numbers up to 10000:')
print( numberlist[largestPrimeFact<=5] )

# check against https://oeis.org/A051037
```


```python
sym.primefactors(4374)
```


```python
# adjust the previous code to plot the count of k-smooth numbers

# count the number of times each smooth number appears
smoothCount = np.zeros(len(uniqueMaxPrimes),dtype=int)
for i,u in enumerate(uniqueMaxPrimes):
    smoothCount[i] = len(numberlist[largestPrimeFact<=u])


# and visualize the counts by the smooth numbers
plt.plot(uniqueMaxPrimes,smoothCount,'k.')
plt.xlabel('Smooth number')
plt.ylabel('Number of occurrences')
plt.xlim([-1,500])
plt.title('K-smooth number frequencies')
plt.show()
```


```python

```

ee9c73ad7e699a4cc102be742a1d26f409d734aa
ipynb
Jupyter Notebook
MXC_pymath_numberTheory.ipynb
stefan-cross/mathematics-with-python
40a130b0d7da98c4bd54406c61e99daeca57680f
MIT
MXC_pymath_numberTheory.ipynb
stefan-cross/mathematics-with-python
40a130b0d7da98c4bd54406c61e99daeca57680f
MIT
MXC_pymath_numberTheory.ipynb
stefan-cross/mathematics-with-python
40a130b0d7da98c4bd54406c61e99daeca57680f
MIT
Qwen/Qwen-72B
1. YES
2. YES