source('utils.R')


# To check our codes, we duplicate the results in Example 3.8 on EPage 128
#
an =  c(-1, 1, -1, -1, 1)
print('For the polynomial with coefficients:')
print(an)
polynomial_root_bounds(an)


# To check our codes, we duplicate the results in Example on Page 112
#
an = c(-1, -1, 1, 1)
print('For the polynomial with coefficients: ')
print(an)
polynomial_root_bounds(an)


# To check our codes, we duplicate the results in Example 3.10 on EPage 130
#
ep_res = evaluate_polynomial(1.1, c(-3, 1, 0, 1))
print(ep_res)

pn_res = polynomial_newton(1.1, c(-3, 1, 0, 1))
print(pn_res)


# To check our codes, we duplicate the results in Example 3.11 on EPage 131
#
ep_res = evaluate_polynomial(1.5, c(-6.8, 10.8, -10.8, 7.4, -3.7, 1))
print(ep_res)

pn_res = polynomial_newton(1.5, c(-6.8, 10.8, -10.8, 7.4, -3.7, 1))
print(pn_res)



# Exercise 3.6-1:
#
an = c(-1, 2, 0, 1)
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

xs = seq(0, 3, length.out=100)
ps = vector_evaluate_polynomial(xs, prb$P)
qs = vector_evaluate_polynomial(xs, prb$Q)

plot(xs, ps, type='l', col='blue', xlab='x', ylab='', main='P(x) and Q(x)')
lines(xs, qs, type='l', col='red')
grid()
legend('topleft', c('P(x)','Q(x)'), lty=c(1, 1), col=c('blue', 'red'))

# Based on the above we have:
#
r = 0.5
R = 1.6
pn_res = polynomial_newton(0.5*(r+R), an)
print(sprintf('Ex 3.6-1 x= %f delta_x= %g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))

# This gives coefficients of q(x) where p(x) = a_0' + (x-z) q(x):
#
qns = evaluate_polynomial(pn_res[1], an)
print(polyroot( qns[2:4] ))


# Exercise 3.6-2 (Part (a)):
#
an = c(-1, 0, 1, -3, 0, 1)
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

# Plot the polynomial function we want to find the roots to:
#
xs = seq(-2, +2, length.out=1000)
pxs = vector_evaluate_polynomial(xs, an)
plot(xs, pxs, type='l', col='black', xlab='x' , ylab='', main='Ex. 3.6-2 Part (a) Polynomial')
grid()

# Plot the polynomials P(x) and Q(x):
#
xs = seq(0, 2, length.out=100)
ps = vector_evaluate_polynomial(xs, prb$P)
qs = vector_evaluate_polynomial(xs, prb$Q)

plot(xs, ps, type='l', col='blue', xlab='x', ylab='', main='P(x) and Q(x)')
lines(xs, qs, type='l', col='red')
grid()
legend('topleft', c('P(x)','Q(x)'), lty=c(1, 1), col=c('blue', 'red'))

# Based on the above we have:
#
r = 0.0
R = 1.75
pn_res = polynomial_newton(0.5*(r+R), an)
print(sprintf('Ex 3.6-2 (a): x= %f delta_x= %g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))

# This gives coefficients of q(x) where p(x) = a_0' + (x-z) q(x):
#
qns = evaluate_polynomial(pn_res[1], an)

an_d1 = qns[2:length(qns)] # the deflated polynomial with the first root removed 
print('For the first deflated polynomial with coefficients: ')
print(an_d1)
prb = polynomial_root_bounds(an_d1)

# Plot the first deflated polynomial:
#
xs = seq(- 2, +2, length.out=100)
pxs = vector_evaluate_polynomial(xs, an_d1)
plot(xs, pxs, type='l', col='black', xlab='x' , ylab='', main='Ex. 3.6-2 Part (a) Polynomial/(x - 1.618035)' )
grid()


# Exercise 3.6-2 (Part (b)):
#
an = c(-1, 3, 0, 1)
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

# Plot the polynomial function we want to find the roots to:
#
xs = seq(-3, +3, length.out=1000)
pxs = vector_evaluate_polynomial(xs, an)
plot(xs, pxs, type='l', col='black', xlab=' x' ,ylab='', main='Ex. 3.6-2 Part (b) Polynomial' )
grid()

# Plot the polynomials P(x) and Q(x):
#
xs = seq(0, 2, length.out=100)
ps = vector_evaluate_polynomial(xs, prb$P)
qs = vector_evaluate_polynomial(xs, prb$Q)
plot(xs, ps, type='l', col='blue', xlab='x', ylab='', main='P(x) and Q(x)')
lines(xs, qs, type='l', col='red')
grid()
legend('topleft', c('P(x)','Q(x)'), lty=c(1, 1), col=c('blue', 'red'))

# Based on the above we have:
#
r = 0.33
R = 1.75
pn_res = polynomial_newton(0.5*(r+R), an)
print(sprintf('Ex 3.6-2 (b): x= %f delta_x= %g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))

# This gives coefficients of q(x) where p(x) = a_0' + (x-z) q(x):
#
qns = evaluate_polynomial(pn_res[1], an)

an_d1 = qns[2:length(qns)] # the deflated polynomial with the first root removed
print('For the first deflated polynomial with coefficients: ')
print(an_d1)

prb = polynomial_root_bounds(an_d1)


# Exercise 3.6-3:
#
an = c(-4.2, -6.3, -0.38, 2.8, 1)
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

# Plot the polynomial function we want to find the roots to:
#
xs = seq(-3, +2, length.out=1000)
xs = seq(-2, -1, length.out=1000)
pxs = vector_evaluate_polynomial(xs, an)
plot(xs, pxs, type='l', col='black', xlab='x', ylab='', main='Ex. 3.6-3 Polynomial')
grid()

pn_res = polynomial_newton(-1.7, an)
print(sprintf('Ex 3.6-3 root #1: x= %+f delta_x= %+g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))
pn_res = polynomial_newton(-1.5, an)
print(sprintf('Ex 3.6-3 root #2: x= %+f delta_x= %+g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))
pn_res = polynomial_newton(-1.1, an)
print(sprintf('Ex 3.6-3 root #3: x= %+f delta_x= %+g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))
pn_res = polynomial_newton(1.2, an)
print(sprintf('Ex 3.6-3 root #4: x= %+f delta_x= %+g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))


# Exercise 3.6-4:
#
an = c(51200, 0, -39712, 0, 7392, 0, -170, 0, 1)
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

true_roots = c( -10, -8, -2, -sqrt(2), +sqrt(2), 2, 8, 10 )
approx_roots = true_roots * 1.1

for( ri in 1:length(true_roots) ){
    x0 = approx_roots[ri]
    pn_res = polynomial_newton(x0, an)
    print(sprintf('Ex 3.6-4 (original): x0= %f; x= %f; delta_x= %g; n= %.1f', x0, pn_res[1], pn_res[3], pn_res[2]))
}

# Change the coefficient in this polynomial:
#
an[3] = -39710
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

for( ri in 1:length(true_roots) ){
    x0 = approx_roots[ri]
    pn_res = polynomial_newton(x0, an)
    print(sprintf('Ex 3.6-4 (modified): x0= %f; x= %f; delta_x= %g; n= %.1f', x0, pn_res[1], pn_res[3], pn_res[2]))
}



# Exercise 3.6-5:
#
an = c(1, -1, 1, -1, 1)
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

# Plot the polynomial function we want to find the roots to:
#
xs = seq(-1, +1.5, length.out=1000)
pxs = vector_evaluate_polynomial(xs, an)
plot(xs, pxs, ylim=c(0, 5), type='l', col='black', xlab='x', ylab='', main='Ex. 3.6-5 Polynomial')
grid()

# Plot the polynomials P(x) and Q(x):
#
xs = seq(0, 2, length.out=100)
ps = vector_evaluate_polynomial(xs, prb$P)
qs = vector_evaluate_polynomial(xs, prb$Q)

plot(xs, ps, type='l', col='blue', xlab='x', ylab='', main='P(x) and Q(x)')
lines(xs, qs, type='l', col='red')
grid()
legend('topleft', c('P(x)','Q(x)'), lty=c(1, 1), col=c('blue', 'red'))

# Based on the above we have:
#
r = 0.5
R = 1.95
pn_res = polynomial_newton(0.5*(r+R), an)
print(sprintf('Ex 3.6-5: x= %f delta_x= %g n= %.1f', pn_res[1], pn_res[3], pn_res[2]))


# Exercise 3.6-6:
#
an = c(6.3504, 0, -5.85, 0, 1)
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

x0 = c(1.2, 2)
maehly_res = maehly_iteration(an, x0)
R = data.frame(x0=x0, x=maehly_res$x)
print(R)


# Exercise 3.6-8:
#
an = c(-5040, 13068, -13132, 6769, -1960, 322, -28, 1) # the first polynomial
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

x0 = c(0.9, 1.9, 2.9, 3.9, 4.9, 5.9, 6.9)
maehly_res = maehly_iteration(an, x0)
R = data.frame(x0=x0, x=maehly_res$x)
print(R)

an = c(-32, 124, -155, 77.5, -15.5, 1) # the second polynomial
print('For the polynomial with coefficients: ')
print(an)
prb = polynomial_root_bounds(an)

x0 = c(0.45, 0.9, 1.8, 3.6, 7.2)
maehly_res = maehly_iteration(an, x0)

R = data.frame(x0=x0, x=maehly_res$x)
print(R)