source('utils.R')

# 2.1-1:
#
x = 314.15

newton_form = function(x){
    (x-99*pi) * (x-100*pi) * (x-101*pi)
}
power_form = function(x){
    x^3 - 300*pi*x^2 + 29999*pi^2*x - 999900*pi^3
}
print(sprintf('product form= %f; power form= %f', newton_form(x), power_form(x)))

# 2.1-2:
#
p_fn = function(x){
    p = 3 - (x-1)*(4 - (x+1)*(5 - x*(6 - (x+2))))
    return(p)
}
print(sprintf('p(x=2)= %f', p_fn(2)))

as = c( 3, -4, 5, -6, 1 )
cs = c( 1, -1, 0, -2 )
x = 2.0

nf = newton_form_nested_multiplication(as, cs, x)
print(nf)
asprime = nf$as_prime
csprime = nf$cs_prime # the new centers

nf = newton_form_nested_multiplication(asprime, csprime, x)
print(nf)

# 2.1-4:
#
# Using the above results from the previous problem we have
#
asprime = nf$as_prime
csprime = nf$cs_prime
nf = newton_form_nested_multiplication(asprime, csprime, x)
print(nf)

# 2.1-5:
#
# Get the coefficients of the Taylor seriese with respect to x=3
#
x=3

as = c( 3, -4, 5, -6, 1) # the original set of Newton form coefficients for p(x)
cs = c( 1, -1, 0, -2 )
poly_degree = length(cs)

for( neval in 1:poly_degree ) {
    nf = newton_form_nested_multiplication(as, cs, x)
    as = nf$as_prime
    cs = nf$cs_prime # the new centers with on value x appended
}
print('Final centers:')
print(cs) # the final centers (should be all x=3 values)
print('Final Netwon form coefficients:')
print(as) # the Taylor series coefficients


# 2.1-9:
#

# Part (a):
#
exp_via_nested_multiplication = function(x, N=4){
    stopifnot(N>=1)

    p = x/N + 1 # the first factor in the product
    if(N==1){ return(p) }
    for(ii in seq(N-1, 1, -1)){
        p = p * (x/ii) + 1
    }
    return (p)
}

x = 0.01
print(sprintf('x= %f; exp(x)= %f; exp_via_nested_multiplication(x)= %f', x, exp(x), exp_via_nested_multiplication(x)))
x = 0.1
print(sprintf('x= %f; exp(x)= %f; exp_via_nested_multiplication(x)= %f', x, exp(x), exp_via_nested_multiplication(x)))

# Part (b):
#
log_via_nested_multiplication = function(x, N=4){
    stopifnot(N>=1)
    
    r = (x-1)/(x+1)
    
    p = x/(2*N+1) + 1/(2*N-1) # the first factor in the product 
    if(N==1){ return(2*r*p) }
    for(ii in seq(N-1, 1, -1)){
        p = p * r^2 + 1/(2*ii-1)
    }
    return (2*r*p)
}

x = 1.01
print(sprintf('x= %f; log(x)= %f; log_via_nested_multiplication(x)= %f', x, log(x), log_via_nested_multiplication(x)))
x =1.1
print(sprintf('x= %f; log(x)= %f; log_via_nested_multiplication(x)= %f', x, log(x), log_via_nested_multiplication(x)))

# Part (c):
#
asin_via_nested_multiplication = function(x, N=20){
    stopifnot(N>=1)

    p = ( x^2 / (2*N+1) ) * (2*N-1)/(2*N) + 1/(2*N-1)

    if(N==1){ return(x*p) }
    for(ii in seq(N-1, 1, -1)){
        p = p * x^2 * (2*ii-1)/(2*ii) + 1/(2*ii-1)
    }
    return (x*p)
}

x = 0.01
print(sprintf('x= %f; asin(x)= %f; asin_via_nested_multiplication(x)= %f', x, asin(x), asin_via_nested_multiplication(x)))
x = 0.1
print(sprintf('x= %f; asin(x)= %f; asin_via_nested_multiplication(x)= %f', x, asin(x), asin_via_nested_multiplication(x)))
x = 0.5
print(sprintf('x= %f; asin(x)= %f; asin_via_nested_multiplication(x)= %f', x, asin(x), asin_via_nested_multiplication(x)))
x = 0.75
print(sprintf('x= %f; asin(x)= %f; asin_via_nested_multiplication(x)= %f', x, asin(x), asin_via_nested_multiplication(x)))
x = 0.98 # we will need many terms to evaluate asin accurately at this value of x
print(sprintf('x= %f; asin(x)= %f; asin_via_nested_multiplication(x)= %f', x, asin(x), asin_via_nested_multiplication(x, N=500)))