Homework 8 Part 2.sagews

# Do your scratchwork at the bottom and put your finished code here.
# To receive full credit, these answers at the top must be easy to read.

%md
Who else is in your group? Xiao Hu <br>
Whose worksheet would you like me to grade? Tony Sanchez

Who else is in your group? Xiao Hu
Whose worksheet would you like me to grade? Tony Sanchez

# Here is a sample of how to do some plotting in CoCalc.
x = var("x")
bound = 5000
g = plot(x/log(x),(x,2,bound),color='red')
g = g + plot(Li(x),(x,2,bound),color='black')
for i in range(2,bound,floor(bound/100)):
    g = g + point((i,prime_pi(i)),color='blue')
show(g)

# In part 1 of this homework, we found wrote a function that
# counted points on elliptic curves over Fp in approximately
# O(p) steps.  A mathematician named Schoof found an algorithm
# which counts points on elliptic curves over Fp in O(log(p)^6) steps.
# Is that really an improvement? Graph the function f(x) = x in blue
# and the function g(x) = log(x)^6 in red for 1 <= x <= 10^8.
x = var("x")
bound = 10^8
f = plot(x,(x,1,bound),color='blue')
f = f + plot(log(x)^6,(x,1,bound),color='red')
show(f)

# Fix some non-zero integer value of a and some non-zero integer value of b.
# These values should never change.  Plot y = 2*sqrt(x) and y = -2*sqrt(x).
# Also do the following procedure 500 times (it should take under 30 seconds).
# Choose a random prime p <= 10,000.  If y^2 = x^3 + ax + b
# is an elliptic curve, then use E.cardinality() to count the number of points on
# the elliptic curve, and put (p,number of points - (p+1)) on your plot.  In other
# words, you are plotting the error of our estimate p+1 for the number of points
# for this particular elliptic curve.
# To get full credit, please check to make sure y^2 = x^3 + ax + b is really an
# elliptic curve (as opposed to a singular cubic).
a = -1
b = 5
p = random_prime(10000, lbound = 9000)
x = var("x")
bound = 10000
f = plot(2*sqrt(x),(x,1,bound),color='blue')
f = f + plot(-2*sqrt(x),(x,1,bound),color='blue')


for i in range (0,500):
    p = random_prime(10000, lbound = 3)
    if (4*a^3+27*b^2) % p != 0:
        E = EllipticCurve(GF(p),[a,b])
        f = f + point((p, E.cardinality() -  (p + 1)), color='black')
show(f)

#testing another method
a = -1
b = 5
p = random_prime(10000, lbound = 9000)
x = var("x")
bound = 10000
f = plot(2*sqrt(x),(x,1,bound),color='blue')
f = f + plot(-2*sqrt(x),(x,1,bound),color='blue')


for i in range (0,500):
    p = random_prime(10000)
    try:
        E = EllipticCurve(GF(p),[a,b])
        f = f + point((p, E.cardinality() -  (p + 1)), color='black')
    except:
        pass
show(f)

def newton_method(f,a,b):
    x = var('x')
    df = diff(f,x)
    NewtonProx(x) = x - (f/df)(x)
    xi = a
    for i in range (0,b):
        xi = N(NewtonProx(xi))
        print xi

f(x) = x**4+x**3-1
newton_method(f,2,9)

1.47727272727273 1.11793386558496 0.908134322032024 0.829688956906674 0.819339664556753 0.819172556393136 0.819172513396167 0.819172513396164 0.819172513396164

f(x) = x*exp(x)-2
newton_method(f,0.5,6)

0.975374212950178 0.863359106097814 0.852693923733206 0.852605508032695 0.852605502013726 0.852605502013726

for i in range (0, 1000):
    if is_prime(i):
        for a in range (0,i):
            for b in range (0,i):
                if (4*a^3+27*b^2) % i != 0:
                    print