bisection_method = function(a, b, root_fn, xtol=1.e-6, ftol=1e-3){
    # Algo 3.1: EPage 91
    #
    fa = root_fn(a)
    fb = root_fn(b)
    stopifnot( fa*fb<0 )

    n=0
    done = FALSE
    while( !done ){
        w = b-a # how wide is our interval
        m = 0.5 * ( a + b )
        fm = root_fn(m)
        if( fa*fm <= 0 ){
            b = m
        }else{
            a = m
        }
        if( w<xtol ){
            done = TRUE
        }
        n = n+1
    }
    m = 0.5 * ( a + b )
    fm = root_fn(m)

    return( c(m, fm, n) )
}


regula_falsi = function(a, b, root_fn, xtol=1.e-6, ftol=1.e-3){
    # Algo 3.2: EPage 92
    #
    fa = root_fn(a)
    fb = root_fn(b)
    stopifnot( fa*fb<0 )
    xsize = 0.5*(abs(a) + abs(b))
    fsize = 0.5*(abs(fa) + abs(fb))

    n=0
    done = FALSE
    while( !done ){
        w = ( fb*a - fa*b ) / ( fb - fa )
        fw = root_fn(w)
        if( fa*fw <= 0 ){
            b = w
        }else{
            a = w
        }
        if( abs(b-a)<xtol*xsize || (abs(fw)/fsize)<ftol ){
            done = TRUE
        }
        n = n+1
    }
    w = ( fb*a - fa*b ) / ( fb - fa )
    fw = root_fn(w)

    return( c(w, fw, n) )
}


modified_regula_falsi = function(an, bn, root_fn, xtol=1.e-6, ftol=1.e-3){
    # Algo 3.3: EPage 94
    #
    F = root_fn(an)
    G = root_fn(bn)
    stopifnot( F*G<0 )
    xsize = 0.5*(abs(an)+abs(bn))
    fsize = 0.5*(abs(F)+abs(G))
    wn = an
    fwn = root_fn(wn)
    
    n=0
    done = FALSE
    while( !done ){
        wnp1 = ( G*an - F*bn ) / ( G - F )
        fwnp1 = root_fn(wnp1)
        if( root_fn(an)*fwnp1 <= 0 ){
            anp1 = an
            bnp1 = wnp1
            G = fwnp1
            if( fwn*fwnp1 > 0 ){ F = F/2 }
        }else{
            anp1 = wnp1
            bnp1 = bn
            F = fwnp1
            if( fwn*fwnp1 > 0 ){ G = G/2 }
        }
        an = anp1
        bn = bnp1
        wn = wnp1
        fwn = fwnp1
        if( abs(bn-an)<xtol*xsize || (abs(fwn)/fsize)<ftol ){
            done = TRUE
        }
        n = n+1
    }

    return( c(wn, fwn, n) )
}


secant_method = function(xnm1, xn, root_fn, xtol=1.e-6, ftol=1.e-3){
    # Algo 3.4: EPage 94
    #
    fnm1 = root_fn(xnm1)
    fn = root_fn(xn)
    xsize = 0.5*(abs(xnm1)+abs(xn))
    fsize = 0.5*(abs(fnm1)+abs(fn))
    
    n=0
    done = FALSE
    while( !done ){
        dx = fn*(xnm1-xn)/(fn-fnm1) # See Exercise 3.2-6
        xnp1 = xn + dx
        fnp1 = root_fn(xnp1)
        
        # Update root estimate:
        #
        xnm1 = xn
        fnm1 = fn
        xn = xnp1
        fn = fnp1

        # Check for convergence:
        #
        if( abs(dx)<xtol*xsize || (abs(fnp1)/fsize)<ftol ){
            done = TRUE
        }
        n = n+1
    }

    return( c(xnp1, fnp1, n) )
}


newtons_method = function(xn, root_fn, root_fn_prime, xtol=1.e-6, ftol=1.e-3, ntol=1000){
    # Algo 3.5: EPage 95
    #
    fn = root_fn(xn)
    xsize = abs(xn)
    fsize = abs(fn)
    
    n=0
    done = FALSE
    while( !done ){
        fn_prime = root_fn_prime(xn)
        stopifnot( abs(fn_prime)>1.e-9 ) # add an assert to flag if f'(x_n)=0
        xnp1 = xn - fn / fn_prime 

        # Check for convergence:
        #
        if( abs(xn-xnp1)<xtol*xsize || (abs(fn)/fsize)<ftol || n>=ntol ){
            done = TRUE
        }

        # Update root estimate:
        #
        xn = xnp1
        fn = root_fn(xn)

        n = n+1
    }

    return( c(xn, fn, n) )
}


fixed_point_iteration = function(xn, root_fn, iter_fn, xtol=1.e-6, ftol=1.e-3, ntol=1000){
    # Algo 3.6 (fixed point iteration): EPage 109
    #
    fn = root_fn(xn)
    xsize = abs(xn)
    fsize = abs(fn)

    n=0
    done = FALSE
    while( !done ){
        xnp1 = iter_fn(xn)
        fnp1 = root_fn(xnp1)

        # Check for convergence:
        #
        if( abs(xn-xnp1)<xtol*xsize || (abs(fn)/fsize)<ftol || n>=ntol ){
            done = TRUE
        }

        # Update root estimate:
        #
        xn = xnp1
        fn = fnp1
        n = n+1
    }

    return( c(xn, fn, n) )
}


aitkens_delta2 = function(xn){
    # Algo 3.7: EPage 112 & 114 (we compute r_n and x_hat_n from x_n)
    #
    # xn[1] = x_0 the initial iterate value
    # xn[N+1] = x_N the last iterate
    #
    N = length(xn)-1

    dxn = c( NaN, diff(xn) ) # = x_n - x_{n-1}
    
    dxnp1 = c( diff(xn), NaN ) # holds = x_{n+1} - x_n

    rn = dxn / dxnp1 # = Second Equation 3.27 on EPage 112

    xnp1 = c( xn[2:(N+1)], NaN )

    xhatn = xnp1 + (xnp1 - xn)/(rn-1) # = First Equation 3.72 on EPage 112

    return( list(rn=rn, xhatn=xhatn) )
}


steffensen_iteration = function(y0, iter_fn, xtol=1.e-6, ntol=1000){
    # Algo 3.8: EPage 114
    #
    x0 = y0
    xsize = abs(x0)

    n=0
    done = FALSE
    while( !done ){
        x1 = iter_fn(x0)
        x2 = iter_fn(x1)

        # Check for convergence and update our root estimate
        #
        d = x2 - x1
        if( d==0 ){
            done = TRUE
            x0 = x1
        }else{
            r = (x1 - x0)/d
            y1 = x2 + d/(r-1)

            if( abs(y1-y0)<xtol*xsize || n>=ntol ){
                done = TRUE
            }
            x0 = y1
        }

        n = n+1
    }
    return( c(x0, n) )
}


number_of_sign_changes = function(an){
    # Computes the number of sign changes in the polynomial with coefficients an (ignoring zeros).
    #
    # The coefficients of the polynomial are stored as the elements of the vector an such that
    #
    # p(x) = a[1] + a[2] x + a[3] x^2 + ... + a[n] x^{n-1} + a[n+1] x^n
    #
    an_no_zeros = an[an!=0]
    n = length(an_no_zeros) - 1

    n_sign_changes = 0
    current_sign = sign(an_no_zeros[n+1])
    for( ni in seq(n, 1, -1) ){
        if( current_sign * an_no_zeros[ni] < 0 ){
            n_sign_changes = n_sign_changes+1
            current_sign = sign(an_no_zeros[ni])
        }
    }

    return(n_sign_changes)
}


polynomial_root_bounds = function(an){
    #
    # Returns some information on root bounds from the discussions on EPage 126 - 128
    #
    # The coefficients of the polynomial are stored as the elements of the vector an such that
    #
    # p(x) = a[1] + a[2] x + a[3] x^2 + ... + a[n] x^{n-1} + a[n+1] x^n
    #
    n = length(an) - 1 # the degree of the polynomial

    # Descarts' rule of signs:
    #
    nu_p = number_of_sign_changes(an) # maximal number of positive roots
    nu_n = number_of_sign_changes(an * rep(-1, n+1)^seq(0, n)) # maximal number of negative roots
    print(sprintf('Descarts rule of sign: degree= %d, max positive roots: %d, max negative roots: %d', n, nu_p, nu_n))

    #
    rho_1 = n * abs(an[1]) / abs(an[2])
    rho_n = ( abs(an[1]) / abs(an[n+1]) )^(1./n)
    print(sprintf('At least one zero of this polynomial satisfies: |x| <= %.4f', min(c(rho_1, rho_n))))

    # All zeros of the original polynomial i.e inside r <= le <= R.
    # One would need to plot P(x) and Q(x) to find estimates for the numbers R and r:
    #
    P_polynomial = abs(an)
    P_polynomial[1:n] = -P_polynomial[1:n] # P(x) has one real positive zero R

    Q_polynomial = abs(an)
    Q_polynomial[1] = -Q_polynomial[1] # Q(x) has one real positive zero r
    #

    r = 1 + max( abs( an[1:n] / an[n+1] ) )
    print(sprintf('Every zero of this polynomial satisfies: |x| <= %.4f', r))

    return( list(P=P_polynomial, Q=Q_polynomial) )
}


evaluate_polynomial = function(z, an){
    # Reduces a polynomial as in Equations 3.47 on EPage 128
    #
    # The coefficients of the polynomial are stored as the elements of the vector an such that
    #
    # p(x) = a[1] + a[2] x + a[3] x^2 + ... + a[n] x^{n-1} + a[n+1] x^n
    #
    n = length(an)
    an_prime = rep(NA, n)
    for( ni in seq(n, 1, -1) ){
        if( ni==n ){
            an_prime[ni] = an[ni]
        }else{
            an_prime[ni] = an[ni] + z * an_prime[ni+1]
        }
    }
    return(an_prime)
}

vector_evaluate_polynomial = function(z, an){
    #
    # Evaluates a polynomial over a vector of inputs.
    #
    res = rep(NA, length(z))
    for( zi in 1:length(z) ){
        ep = evaluate_polynomial(z[zi], an)
        res[zi] = ep[1]
    }
    return(res)
}


polynomial_newton = function(xm, an, xtol=1.e-6, mtol=100){
    # Algo 3.9 on EPage 129
    #
    # The coefficients of the polynomial are stored as the elements of the vector an such that
    #
    # p(x) = a[1] + a[2] x + a[3] x^2 + ... + a[n] x^{n-1} + a[n+1] x^n
    #
    n = length(an) - 1 # the degree of the polynomial
    
    an_prime = evaluate_polynomial(xm, an) # compute a_n'

    an_prime_prime = rep(NA, n)

    m = 0
    while( m < mtol ){
        z = xm
        an_prime_prime[n+1] = an_prime[n+1]
        an_prime[n+1] = an[n+1]
        for( k in seq(n, 2, -1) ){
            an_prime[k] = an[k] + z * an_prime[k+1]
            an_prime_prime[k] = an_prime[k] + z * an_prime_prime[k+1]
        }
        an_prime[1] = an[1] + z * an_prime[2]
        delta_x = an_prime[1] / an_prime_prime[2]
        if( abs(delta_x)<xtol*abs(xm) ){ break }
        xm = xm - delta_x

        m = m+1
    }
    return( c(xm, m, delta_x) )
}


maehly_iteration = function(an, x0, xtol=1.e-6, mtol=1000){
    # Discussed on: EPage 135
    # We seek to find the values of k roots of the polynomial with coefficients an
    #
    # The coefficients of the polynomial are stored as the elements of the vector an such that
    #
    # p(x) = a[1] + a[2] x + a[3] x^2 + ... + a[n] x^{n-1} + a[n+1] x^n
    #
    # The initial guesses for each of the roots is given in the vector x0
    #
    n = length(an) - 1 # the degree of the polynomial
    k = length(x0)

    xs = rep(NA, k) # save outputs from each
    ms = rep(NA, k)
    delta_xs = rep(NA, k)

    #
    # Deflate the polynomial one root at a time:
    #

    # Compute the first root by hand:
    #
    ki = 1
    pn_res = polynomial_newton(x0[ki], an, xtol=xtol, mtol=mtol)

    # Save results:
    #
    xs[ki] = pn_res[1]
    ms[ki] = pn_res[2]
    delta_xs[ki] = pn_res[3]

    # Compute each additional root
    #
    for( ki in 2:k ){

        xm = x0[ki]
        m = 0
        while( m < mtol ){

            # Compute p(x_m) and p'(x_m):
            #
            an_prime = evaluate_polynomial(xm, an) # This gives coefficients of q(x) where p(x) = a_0' + (x-z) q(x):
            pxm = an_prime[1]
            qn_prime = evaluate_polynomial(xm, an_prime[2:(n+1)])
            pxm_prime = qn_prime[1]
            #print(sprintf('ki= %d; m= %d; xm= %f; pxm= %f; pxm_prime= %f', ki, m, xm, pxm, pxm_prime))

            # Compute \sum_{r=1}^k p(x_m) / ( x_m - \xi_r ):
            #
            one_over_diff = 1 / ( xm - xs )
            deflated_root_adjustment = pxm * sum( one_over_diff[1:(ki-1)] )
            #print(sprintf('deflated_root_adjustment= %f', deflated_root_adjustment))

            delta_x = pxm / (pxm_prime - deflated_root_adjustment)

            if( abs(delta_x)<xtol*abs(xm) ){ break }
            xm = xm - delta_x

            m = m+1
        }

        # Save results:
        #
        xs[ki] = xm
        ms[ki] = m
        delta_xs[ki] = delta_x
    }
    return( list(x=xs, m=ms, delta_x=delta_xs) )
}

mullers_iteration = function(x0s, root_fn, k, xtol=1.e-6, ftol=1.e-3, mtol=1000, verbose=FALSE){
    #
    # Discussed on: EPage 138
    #
    # We seek to find the values of k roots of the function f(x)
    #
    # The initial guesses for all roots is given in the complex vector x0
    #
    # Note that the roots found can be complex.
    #
    stopifnot(length(x0s)==3)
    if(!is.complex(x0s)){ x0s = as.complex(x0s) }

    xs = rep(NA, k) # save outputs from each
    ms = rep(NA, k)
    delta_xs = rep(NA, k)

    deflated_fn = function(x, xs, ki){
        if( ki==1 ){ ## no deflations initially
            root_fn(x) 
        }else{ ## deflate after each root is found
            denom = prod((x-xs[1:(ki-1)]))
            if( abs(denom)<1.e-6 ){ denom = 1.e-6 } ## prevent "divide by zero"
            root_fn(x)/denom
        }
    }
    
    for( ki in 1:k ){
        #
        # Search for the next root with these 3 approximations:
        #
        xm2 = x0s[1] # x_{i-2} = x_0 when i=2
        fxm2 = deflated_fn(xm2, xs, ki)

        xm1 = x0s[2] # x_{i-1} = x_1 when i=2
        fxm1 = deflated_fn(xm1, xs, ki)

        x = x0s[3] # x_i = x_2 when i=2
        fx = deflated_fn(x, xs, ki)

        h = x - xm1 # h_i = h_2 when i=2
        hm1 = xm1 - xm2 # h_{i-1} = h_1 when i=2

        f_x_xm1 = (fx - fxm1)/h
        f_xm1_xm2 = (fxm1 - fxm2)/hm1
        #print( c(ki, xm2, fxm2, xm1, fxm1, x, fx, h, hm1) )

        m = 0
        while( m < mtol ){

            f_x_xm1_xm2 = ( f_x_xm1 - f_xm1_xm2 ) / (h + hm1)
            ci = f_x_xm1 + h * f_x_xm1_xm2

            rad = ci^2 - 4 * fx * f_x_xm1_xm2
            #print(sprintf('ki= %d; m= %d, f_x_xm1_xm2= %f %+f; #f %+f; rad: %f %+f', ki, m, Re(f_x_xm1_xm2), Im(f_x_xm1), Re(ci), Im(ci), Re(rad), Im(rad)))
            denom_minus = ci - sqrt(rad)
            denom_plus = ci + sqrt(rad)
            if( abs(denom_minus)>abs(denom_plus) ){
                denom = denom_minus
            }else{
                denom = denom_plus
            }
            hp1 = -2 * fx / denom

            if( verbose ){
                print(sprintf('ki= %5d; m= %5d, x= %12.6g %12.6g; f_x_xm1_xm2= %12.6g %12.6g; ci= %12.6g %12.6g; rad= %12.6g %12.6g; denom= %12.6g %12.6g; abs(hp1)= %12.6g; xtol*abs(x)= %12.6g; abs(fx)= %12.6g; ftol= %12.6g', ki, m,
                              Re(x), Im(x),
                              Re(f_x_xm1_xm2), Im(f_x_xm1_xm2),
                              Re(ci), Im(ci),
                              Re(rad), Im(rad),
                              Re(denom), Im(denom), abs(hp1), xtol*abs(x), abs(fx), ftol))
            }

            x_tol_met = abs(hp1)<xtol*abs(x)
            f_tol_met = abs(fx)<ftol
            if( x_tol_met & f_tol_met ){ break }

            xp1 = x + hp1
            fxp1 = deflated_fn(xp1, xs, ki)
            f_xp1_x = (fxp1 - fx)/hp1

            # Update:
            #
            xm2 = xm1
            fxm2 = fxm1

            xm1 = x
            fxm1 = fx

            x = xp1
            fx = fxp1

            hm1 = h
            h = hp1

            f_xm1_xm2 = f_x_xm1
            f_x_xm1 = f_xp1_x

            m = m+1
        }

        # Save results:
        #
        xs[ki] = x
        ms[ki] = m
        delta_xs[ki] = hp1
        if( verbose ){ print('') }
    }
    return( list(x=xs, m=ms, delta_x=delta_xs) )
}