source('utils.R')

# To check our codes, we dup1icate the results in Examp1e 3.14 on EPage 141
#
root_fn = function(x){ x^3 + x - 3 }

x0s = as.complex( c( 0.5, -0.5, 0.0 ) )
k = 3
mullers_res = mullers_iteration(x0s, root_fn, k)
for( ki in 1:k ){
    print(sprintf('Examp1e 3.14: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %.4g n= %5.1f',
                  ki, Re(mullers_res$x[ki]), Im(mullers_res$x[ki]), abs(mullers_res$delta_x[ki]), mullers_res$m[ki]))
}

# 3.7-1 Part (a):
#
root_fn = function(x){ x^7 - 1 }
polynomial_root_bounds(c(-1, 0, 0, 0, 0, 0, 0, 1))

x0s = as.complex( c( 0.25, 0.5, 0.75 ) )
k = 7
mullers_res = mullers_iteration(x0s, root_fn, k)
for( ki in 1:k ){
    print(sprintf('3.7-1 Part (a): Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e n= %5.1f',
                  ki, Re(mullers_res$x[ki]), Im(mullers_res$x[ki]), abs(mullers_res$delta_x[ki]), mullers_res$m[ki]))
}

# 3.7-1 Part (b):
#
root_fn = function(x){ x^4 - 7 * x^3 + 18 * x^2 - 20 * x + 8 }
an = c(8, -20, 18, -7, 1)
polynomial_root_bounds(an)

x0s = as.complex( c( 0.25, 0.5, 0.75 ) )
k = 4
mullers_res = mullers_iteration(x0s, root_fn, k)
for( ki in 1:k ){
    print(sprintf('3.7-1 Part (b): Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e n= %5.1f',
                  ki, Re(mullers_res$x[ki]), Im(mullers_res$x[ki]), abs(mullers_res$delta_x[ki]), mullers_res$m[ki]))
}

# 3.7-1 Part (c):
#
root_fn = function(x){ x^6 + 2 * x^5 + x^4 + 3 * x^3 + 5 * x^2 + x + 1 }
an = c(1, 1, 5, 3, 1, 2, 1)
polynomial_root_bounds(an)

x0s = as.complex( c( 0.25, 0.5, 0.75 ))
k = 6
mullers_res = mullers_iteration(x0s, root_fn, k)
for( ki in 1:k ){
    print(sprintf('3.7-1 Part (c): Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e n= %5.1f',
                  ki, Re(mullers_res$x[ki]), Im(mullers_res$x[ki]), abs(mullers_res$delta_x[ki]), mullers_res$m[ki]))
}


# 3.7-2:
#
root_fn = function(x){ x - tan(x) }

xs = seq(4, 6, length.out=100 )
fns = root_fn(xs)
plot(xs, fns, type='l')
grid()

a = 3*pi/2
f = 0.01
x0s = as.complex( c( a - 3*f*a, a - 2*f*a, a - f*a ) )
mullers_res = mullers_iteration(x0s, root_fn, 1)
print(sprintf('3.7-2: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f',
              1, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))

a = 5*pi/2
f = 0.01
x0s = as.complex( c( a - 3*f*a, a - 2*f*a, a - f*a ) )
mullers_res = mullers_iteration(x0s, root_fn, 1)
print(sprintf('3.7-2: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f',
              2, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))

a = 7*pi/2
f = 0.01
x0s = as.complex( c( a - 3*f*a, a - 2*f*a, a - f*a ) )
mullers_res = mullers_iteration(x0s, root_fn, 1)
print(sprintf('3.7-2: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f',
              3, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))


# 3.7-3:
#
my_C = function(x, n_max=20){
    #
    # See: https://en.wikipedia.org/wiki/Fresnel_integral
    #
    value = 0
    for( n in 0:n_max ){
        numer = (-1)^n * x^(4*n+1) # * (pi/2)^(2*h)
        denom = factorial(2*n) * ( 4*n + 1 )
        term = numer / denom
        value = value + term
    }
    value
}

xs = seq(0, 4, length.out=500)
ys = my_C(xs)
plot(xs, ys, type='l', xlab='x', ylab='c(x)')
grid ()

root_fn = function(x){ my_C(x) - 0.6 }

x0s = as.complex( c( 0.5, 0.75, 1.0 ) )
mullers_res = mullers_iteration(x0s, root_fn, 1)
print(sprintf('3.7-3: Root #%3d: x: %12.6g %12.6g i; abs(delta_x)= %6.3e; n: %5.1f', 1, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))

x0s = as.complex( c( 1.5, 1.75, 2.0 ) )
mullers_res = mullers_iteration(x0s, root_fn, 1)
print(sprintf('3.7-3: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f', 2, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))

x0s = as.complex( c( 2.5, 2.6, 2.7 ) )
mullers_res = mullers_iteration(x0s, root_fn, 1)
print(sprintf('3.7-3: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f', 3, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))


# 3.7-4:
#
my_J1 = function(z, k_max=40){
    value = 0
    for( k in 0:k_max ){
        numer = (-z^2/4)^k
        denom = factorial(k) * factorial(k + 1)
        term = numer / denom
        value = value + term
    }
    (z/2) * value
}

xs = seq(0, 15, length.out=500)
ys = my_J1(xs)

plot(xs, ys, type='l', xlab='x', ylab='J_1(x)')
grid()

x0s = as.complex( c( 3, 4, 5 ) )
mullers_res = mullers_iteration(x0s, my_J1, 1)
print(sprintf('3.7-4: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f',
              1, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))

x0s = as.complex( c( 6, 7, 8 ) )
mullers_res = mullers_iteration(x0s, my_J1, 1)
print(sprintf('3.7-4: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f',
              2, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))

x0s = as.complex( c( 9, 10, 11 ) )
mullers_res = mullers_iteration(x0s, my_J1, 1)
print(sprintf('3.7-4: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f',
              3, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))

x0s = as.complex( c( 11, 12, 13 ) )
mullers_res = mullers_iteration(x0s, my_J1, 1)
print(sprintf('3.7-4: Root #%3d: x= %12.6g %12.6g i; abs(delta_x)= %6.3e; n= %5.1f',
              4, Re(mullers_res$x[1]), Im(mullers_res$x[1]), abs(mullers_res$delta_x[1]), mullers_res$m[1]))