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

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

#Help procedure
Help:=proc():
    print(` Hd1(f,d,x,x0) , Hd(f,d,x,x0,N) `):
    print(` TestHouseholder(f,d,x,x0,eps,N) `):
    print(` Basins(c1,c2,c3,eps,compress:=false) `):
end:

##PROBLEM 2##

#Householder's method.
#
#INPUT:
#an expression f in the variable x
#a positive integer d
#a variable x
#an initial guess x0
#
#OUTPUT: one iteration of the Householder method of order d on f
#with initial guess x0
Hd1:=proc(f, d, x, x0)
    return x0 + normal(d*subs(x=x0, diff(1/f, x$(d-1)))/
        subs(x=x0, diff(1/f, x$d))):
end:

#Householder's method.
#
#INPUT:
#an expression f in the variable x
#a positive integer d
#a variable x
#an initial guess x0
#a positive integer N
#
#OUTPUT: N iterations of the Householder method of order d on f
#with initial guess x0
Hd:=proc(f, d, x, x0, N) local x1, i:
    x1:=x0:
    for i from 1 to N do
        x1:=Hd1(f, d, x, x1):
    od:
    return x1:
end:

##PROBLEM 3##

#Test Householder Methods.
#
#INPUT:
#an expression f in the variable x
#a positive integer d
#a variable x
#an initial guess x0
#a small eps
#a positive integer N
#
#OUTPUT: the successive errors
#abs(xn-alpha) for n from 1 to N where alpha is the
#approximation to the root of f(x)=0 found by the bisection method
TestHouseholder:=proc(f, d, x, x0, eps, N) local i, L, alpha, x1:
    alpha:=B(f, x, x0-2, x0+2, eps):
    if alpha = FAIL then
        return FAIL:
    fi:

    L:=[]:
    x1:=x0:
    for i from 1 to N do
        x1:=Hd1(f, d, x, x1):
        L:=[op(L), abs(x1-alpha[1])]:
    od:
    return L:
end:
#Some experimentation shows that the Householder method of order d
#does in fact converge at a rate of d+1

##PROBLEM 5##

#Find basins of attraction for Newton's method
#
#INPUT:
#integers c1, c2, and c3
#resolution eps
#boolean value compress (default is false)
#
#OUTPUT: the list [S1, S2, S3] such that if you apply Newton
#5 times (with floating point) to f=(x-c1)*(x-c2)*(x-c3)
#Si is the set of numbers that are multiples of eps between c1-1
#and c3+1 that converge to ci
#
#If compress is true, the sets will actually be lists of intervals
#
#NOTE: I modified NR1 to return FAIL if it would divide by 0
#NOTE: This may return the wrong point of convergence for a minute
#fraction of values (where 5 iterations isn't enough to tell)
Basins:=proc(c1, c2, c3, eps, compress:=false)
local S, j, startv, endv, i, test, f, N, d1, d2, d3, dm, last:
    startv:=ceil((c1-1)/eps)*eps:
    endv:=floor((c3+1)/eps)*eps:
    f:=(x-c1)*(x-c2)*(x-c3):
    N:=5:
    if compress then
        S:=[[], [], []]:
    else
        S:=[{}, {}, {}]:
    fi:
    for i from startv by eps to endv do
        test:=NR(f, x, i, N):
        if test <> FAIL then
            d1:=abs(test-c1):
            d2:=abs(test-c2):
            d3:=abs(test-c3):
            dm:=min(d1, d2, d3):
            if dm = d1 then
                j:=1:
            elif dm = d2 then
                j:=2:
            else
                j:=3:
            fi:
            if compress then
                if S[j] = [] then
                    S[j]:=[[i, i]]:
                else
                    last:=S[j][-1][-1]:
                    if last = i - eps then
                        S[j][-1][-1]:=i:
                    else
                        S[j]:=[op(S[j]), [i,i]]:
                    fi:
                fi:
            else
                S[j]:=S[j] union {i}:
            fi:
        fi:
    od:
    return S:
end:
#Basins(-3, 1, 4, 1e-4, true) returned
#[
#[[-4.0000, -1.3610], [-0.8488, -0.7953], [-0.7846, -0.7833],
#[2.3529, 2.3530], [2.3544, 2.3605], [2.4100, 2.6942]],
#[[-0.7829, 2.3528]],
#[[-1.3609, -0.8489], [-0.7952, -0.7847], [-0.7832, -0.7830],
#[2.3531, 2.3543], [2.3606, 2.4099], [2.6943, 5.0000]]
#]
#which illustrates the Wikipedian behavior