Fangyao - Solution 2.sagews

# Problem set 2

reset

<function reset at 0x7fc4d86bb050>

#Exercise 1 (demo): Goldbach conjecture
def my_goldbach(n): #n must be even
    L=prime_range(n/2)
    S=[]
    for c in L:
        if is_prime(n-c):
            S.append([c,n-c])
    return S
my_goldbach(24)

[[5, 19], [7, 17], [11, 13]]

#Exercise 2(a): Eratosthenes algorithm (template)
reset()
#Simple code snippet
#Input: list of integers L and an integer m
# Output: list L with multiples of m removed
def remove_multiple(m,L):
    for c in L:
        if (c % m)==0:
            L.remove(c)
    return L
#Example
L=[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
remove_multiple(2,L)

[3, 5, 7, 9, 11, 13, 15, 17, 19]

remove_multiple(L[0],L)# Think what L[0] is in the input L

[5, 7, 11, 13, 17, 19]

#Exercise 2: Eratosthenes algorithm
def simple_erato(n):#Fill in the blanks
    S=[] # S is a placeholder for primes <=k that will be removed
    L=[j for j in range(2,n+1)]# list of numbers <=n with L[0]=2; L will be changing
    k=floor(sqrt(n)) # Fact:for any composite number less or equal n there is a prime factor <=k       #floor=less or equal than
    while L[0]<=k:
        S.append(L[0]) #Save the current prime in this list
        remove_multiple(L[0],L)
    return S+L
simple_erato(50)

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

[j for j in range(2,51)] #This is how you make alist of integers <=50

[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]

#Alternative template for Eratosthenes algorithm
def fancy_erato(n):
    S=[] # S is a placeholder for primes <k
    L=[j for j in range(2,n+1)]# list of numbers from <=n with L[0]=2; L will be changing
    k=floor(sqrt(n)) # for any composite number less or equal n there is a prime factor <=k
    while L[0]<=k:
        S.append(L[0])
        L=filter(lambda a: a % L[0]!=0, L) # compact form of our function remove_multiple
        # != means "not equal"; when a%L[0]!=0 is False, the element a is filtered (removed)
    S=S+L
    return S

fancy_erato(50)

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

#Exercise 3 (template): Fermat's factorization algorithm is based on the following fact:
#Every odd number N is the difference of squares, N=a^2-b^2=(a-b)*(a+b).
def fermat_factors(N): # N should be odd
    a = math.ceil(sqrt(N))# calling math is important. It returns a as decimal number
    b = math.sqrt(a^2-N)
    print(a,b)
    while floor(b)-b!=0:#The condition should verify if b is an integer    #integer: floor(b)-b!=0
        a = a + 1 # Alternative Python-like form of this command: a+=1
        b =math.sqrt(a^2-N)
        print a,b
    print "a =", a
    print "b =", b
    print "c =", a+b
    print "d =", a-b
#Test
fermat_factors(51)

(8.0, 3.605551275463989) 9.0 5.47722557505 10.0 7.0 a = 10.0 b = 7.0 c = 17.0 d = 3.0

prime_pi(7);Li(7).n()

4 3.71188798588011

#Exercise 4: Counting primes
#Directions
#Sage uses notation prime_pi for Prime Counting Function. Sage also knows function Li.
#Define the function x/log(x). Make plots of the the three functions and store them under names of your choice. Include legend into each plot. Then plot all three plots in one figure using command show.f(x)=1-x; g(x)=x^2-2*x
f(x)=x/log(x)
p=plot(f(x),3,1500,legend_label='f(x)',color='green',thickness=2)
q=plot(prime_pi(x),3,1500,legend_label='pi(x)',color='red',thickness=2)
r=plot(Li(x),3,1500,legend_label='Li(x)',color='blue',thickness=2)
show(p+q+r,figsize=[4,2])

9/63

1/7

#Exercise 8: Constructing one PPT (template)
#Input: a rational number r, 0<r<1
#Output: legs of the ppt defined by
def one_ppt(r):
    m=numerator(r)
    n=denominator(r)
    if m%2==0 and n%2==0:
        print('both m and n are odd')
    else:
        a=2*m+1;
        b=2*n+1;
        return [a,b,N(sqrt(a^2+b^2))]
one_ppt(2/3)

[5, 7, 8.60232526704263]

-------------------------------------------------------------------------------
#lab 1 solution


#Lab 1: Plotting legs (a,b) of Primitive Pythagorean Triples


#1. The function makes a list of pairs [m,n] of rationals m/n from N diagonals of the table of rationals
#//Input:?
#//Output:?
def mn_pairs(N):
    L=[];
    for n in range(2,N+2):
        temp=[[m,n] for m in range(1,n)]
        L.extend(temp)
    return L

#Example
L=mn_pairs(10);L

[[1, 2], [1, 3], [2, 3], [1, 4], [2, 4], [3, 4], [1, 5], [2, 5], [3, 5], [4, 5], [1, 6], [2, 6], [3, 6], [4, 6], [5, 6], [1, 7], [2, 7], [3, 7], [4, 7], [5, 7], [6, 7], [1, 8], [2, 8], [3, 8], [4, 8], [5, 8], [6, 8], [7, 8], [1, 9], [2, 9], [3, 9], [4, 9], [5, 9], [6, 9], [7, 9], [8, 9], [1, 10], [2, 10], [3, 10], [4, 10], [5, 10], [6, 10], [7, 10], [8, 10], [9, 10], [1, 11], [2, 11], [3, 11], [4, 11], [5, 11], [6, 11], [7, 11], [8, 11], [9, 11], [10, 11]]

# 2.The function makes a list of legs [a,b] of PTs and excludes PTs that are not primitive (i.e., pairs with a and b are not mutuallly prime)
#//Input:
#//Output:number of PTs in L,number of PPT in L,list of PPTs
def ppt_legs(L):
    p_triples=[[c[1]^2-c[0]^2,2*c[0]*c[1]] for c in L ]#list of PTs
    pp_triples=[t for t in p_triples if gcd(t[0],t[1]) == 1]#list of PPTs
    return(len(p_triples),len(pp_triples),pp_triples)

n_pt,n_ppt,S=ppt_legs(L);n_pt #Convenient way to get multiple output

55

n_ppt

27

#Example of plotting points in SageMath
point(S,size=30,color='red',figsize=[4,3])

#3. Main code
#//Input: Number of diagonals in the table of rationals
#//Output:?
def plot_ppt(N):
    L=mn_pairs(N)
    n_pt,n_ppt,S=ppt_legs(L)
    G=point(S,size=30,color='red',figsize=[5,4])
    G.show()
    return("Number of PTs:", n_pt, "Number of PTs:",n_ppt)

plot_ppt(20)

('Number of PTs:', 210, 'Number of PTs:', 92)

#Example of alternative plotting points in S using Python libraries
# Get data for plotting points
x=[s[0] for s in S];y=[s[1] for s in S]

import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
figure(num=None, figsize=(5, 3))
figure=plt.plot(x, y,'ro',markersize=4)
plt.show(figure)

#Advanced nice snippet of code for plotting (just FYI)
import numpy as np
import matplotlib.pyplot as plt

#evenly sampled time at 200ms intervals
t = np.arange(0., 5., 0.2)

# red dashes, blue squares and green triangles
plt.plot(t, t, 'r--', t, t**2, 'bs', t, t**3, 'g^')
plt.show()

[<matplotlib.lines.Line2D object at 0x7f8117accb50>, <matplotlib.lines.Line2D object at 0x7f8117acca90>, <matplotlib.lines.Line2D object at 0x7f8117abb150>]

#For instructor_________________________________
#Data for detecting missed points
xblue=[195,187,171,115,75,27];yblue=[28,84,140,252,308,364]
#figure=plt.plot(xblue,yblue,'bo',markersize=4)# Include into the plotting command code above

#4 Find the missing point in the figure created in Part 3 (see directions in Chapter 2)
gcd(147,196)

49