source('../Chapter3/utils.R')
source('utils.R')

## Get the matrix B from Example 4.17:
##
n = 6
n = 20
B= build_tridiagonal_matrix(n, 0, 0.5, 0.5)
if( n<= 10 ){
    print('B=')
    print(B)
}

## Lets get an upper bound on |lambda|:
##
abs_lambda_max = norm(B, '1')
print(sprintf('abs_lambda_max= %10.6f', abs_lambda_max) )

## For a range of possible values how many eigenvalues are less than this value:
##
xs = seq(-abs_lambda_max, abs_lambda_max, length.out=3*n)
nscs = c()
for(x in xs){
    res = evaluate_symmetric_tridiagonal_characteristic_polynomial(x, diag(B), diag(B[-nrow(B), -1]))
    nscs = c(nscs, number_of_sign_changes(res$ps))
}
num_ev_lt_x = data.frame(x=xs, number_of_eigenvalues_lt_x=nscs)
print(sprintf('Number_of_eigenvalues less than x (n= %d):', n))
print(num_ev_lt_x)

root_fn = function(x){
    res = evaluate_symmetric_tridiagonal_characteristic_polynomial(x, diag(B), diag(B[-nrow(B), -1]))
    res$ps[length(res$ps)] ## return the value of p_n(x)
}
root_fn_prime = function(x){
    res = evaluate_symmetric_tridiagonal_characteristic_polynomial(x, diag(B), diag(B[-nrow(B), -1])) ## <- really inefficient ... calling the same working function twice
    res$pprimes[length(res$pprimes)] ## return the value of p_n'(x)
}

## Lets call Newton's algorithm on each interval where we have a new root:
##
newton_roots = c()
extra_root = which(diff(num_ev_lt_x$number_of_eigenvalues_lt_x)==1) + 1
stopifnot(length(extra_root)==n) ## we need to have n roots
for(ii in 1:length(extra_root)){
    x_left = num_ev_lt_x[extra_root[ii]-1, ]$x
    x_right = num_ev_lt_x[extra_root[ii], ]$x
    xn = mean(c(x_left, x_right))
    res = newtons_method(xn, root_fn, root_fn_prime)
    newton_roots = c(newton_roots, res[1])
}
print('Eigenvalues found with Newtons method=')
print(newton_roots)