#Nathan Fox
#Homework 1
#I give permission for this work to be posted online

#Read procedures from class
read(`C1.txt`):

#Help procedure
Help:=proc():
    print(` NR(f, x, x0, eps) , CompF(N, F, iters) `):
    print(` CompM(N, iters:=10) , CompM(K, iters:=10) `):
    print(` NormalizeListPoly(L, p) `):
end:

##PROBLEM 4##

#Newton's method.
#
#INPUT:
#an expression f in the variable x
#a variable x
#an initial guess x0
#a small positive number eps
#
#OUTPUT: the root of f(x)=0 to error less than eps.
#
#NOTE: if the function f has no roots or if you end up
#at a point with a horizontal tangent, the behavior is undefined
NR:=proc(f, x, x0, eps) local fp, current, last, num, denom:
    fp:=diff(f, x):
    last:=x0:
    while true do
        current:=last - subs(x=last, f)/subs(x=last, fp):
        if abs(current - last) < eps then
            return current:
        fi:
        last:=current:
    od:
end:

##PROBLEM 5##

#INPUT:
#a positive integer N
#a function F of two positive integer arguments
#a positive integer iters
#
#OUTPUT:
#a list of length N, whose i-th entry is the average running
#time of F(a,b) of iters randomly-drawn pairs of
#positive integers with d (base 10) digits
CompF:=proc(N, F, iters) local L, i, j, r:
    L:=[]:
    for i from 1 to N do
        r:=rand(10^(i-1)..10^i-1):
        L:=[op(L), add(time(F(r(), r())), j=1..iters)/iters]:
    od:
    return L:
end:

#INPUT:
#a positive integer N
#a positive integer iters (optional, default is 10)
#
#OUTPUT:
#a list of length N, whose i-th entry is the average running
#time of M(a,b) (from C1.txt) of iters randomly-drawn pairs of
#positive integers with d (base 10) digits
CompM:=proc(N, iters:=10):
    return CompF(N, M, iters):
end:

#INPUT:
#a positive integer N
#a positive integer iters (optional, default is 10)
#
#OUTPUT:
#a list of length N, whose i-th entry is the average running
#time of K(a,b) (from C1.txt) of iters randomly-drawn pairs of
#positive integers with d (base 10) digits
CompK:=proc(N, iters:=10):
    return CompF(N, K, iters):
end:

#INPUT:
#a list L of numbers
#a number p
#
#OUTPUT:
#the list, where the ith element is divided by i^p
NormalizeListPoly:=proc(L, p) local i:
    return [seq(L[i]/i^p, i=1..nops(L))]:
end:

#CompM(20) returned
#[0., 0.0009000000000, 0.003900000000, 0.005200000000,
#0.01420000000, 0.01760000000, 0.02170000000, 0.02360000000,
#0.04780000000, 0.05950000000, 0.08470000000, 0.07800000000,
#0.08680000000, 0.1054000000, 0.09650000000, 0.1149000000,
#0.1655000000, 0.2073000000, 0.2327000000, 0.2489000000]
#
#NormalizeListPoly(L, 2) (where L is the list above) returned
#[0., 0.0002250000000, 0.0004333333333, 0.0003250000000,
#0.0005680000000, 0.0004888888889, 0.0004428571429,
#0.0003687500000, 0.0005901234568, 0.0005950000000,
#0.0007000000000, 0.0005416666667, 0.0005136094675,
#0.0005377551020, 0.0004288888889, 0.0004488281250,
#0.0005726643599, 0.0006398148148, 0.0006445983380,
#0.0006222500000]
#These numbers are approximately constant, as expected

#CompK(30) returned
#[0., 0.001700000000, 0.003800000000, 0.007300000000,
#0.01380000000, 0.01560000000, 0.02040000000, 0.02330000000,
#0.04690000000, 0.03610000000, 0.03830000000, 0.04610000000,
#0.05270000000, 0.07610000000, 0.06630000000, 0.07240000000,
#0.09610000000, 0.08210000000, 0.08820000000, 0.1152000000,
#0.1002000000, 0.1047000000, 0.1283000000, 0.1152000000,
#0.1415000000, 0.1285000000, 0.1545000000, 0.1341000000,
#0.1624000000, 0.1497000000]
#
#NormalizeListPoly(L, log[2](3)) (where L is the list above)
#returned
#[0., 0.0005666666664, 0.0006661381391, 0.0008111111110,
#0.001076559010, 0.0009115574534, 0.0009336443304,
#0.0008629629623, 0.001441233169, 0.0009387386540,
#0.0008563088236, 0.0008979230469, 0.0009041743311,
#0.001160953162, 0.0009066778060, 0.0008938271596, 0.001077725076,
#0.0008409754308, 0.0008292626013, 0.0009985474878,
#0.0008038957288, 0.0007802918526, 0.0008911259687,
#0.0007479445768, 0.0008611403662, 0.0007348918502,
#0.0008322821727, 0.0006819264961, 0.0007811600718,
#0.0006824015463]
#These numbers are approximately constant, as expected

##PROBLEM 6##
#I played around with this for about an hour and managed to get
#6 multiplications a number of different ways, but I never could
#get 5.  I decided to go to the literature on the problem, and I
#found this paper:
#http://www.ei.rub.de/media/crypto/veroeffentlichungen/2010/08/08/kaweb.pdf
#These authors, three years ago, obtain 6 multiplications for
#this problem, and their algorithm becomes progressively less
#efficient as the number of terms increases.