newton_form_nested_multiplication = function(as, cs, z){
    #
    # This is Algorithm 2.1.
    #
    # EPage 49
    #
    # Inputs:
    #
    # as = [a0, a1, ..., an]
    # cs = [c1, c2, ..., cn]
    #
    # Outputs:
    #
    # asprime =[a0', a1', ..., an']
    # csprime =[c1', ..., cn'] = [x, c1, c2, ..., c_{n-1}]
    #
    n = length(as)-1
    stopifnot(n == length(cs))

    as_prime = rep(NA, n+1)
    as_prime[n+1] = as[n+1]
    for( ii in seq(n, 1, -1) ){
        as_prime[ii] = as[ii] + (z - cs[ii]) * as_prime[ii+1]
    }
    cs_prime = c(z, cs[-length(cs)]) # the new centers

    return(list(as_prime=as_prime, cs_prime=cs_prime) )
}


coefficients_of_newton_form = function(xs, ds){
    ##
    ## Implements Algorithm 2.3.
    ##
    ## Inputs:
    ##
    ## xs = [x_0, x_1, ..., x_n]
    ## ds = [f(x_0), f(x_1), ..., f(x_n)]
    ##
    ## Outputs the newton_form of the interpolating coefficients:
    ##
    ## ds = [f[x_0,x_1,...,x_n], f[x_1,x_2,...,x_n], ..., f[x_{n- 1},x_n], f[x_n]]
    ##
    ## The interpreting polynomial is then given by:
    ##
    ## f[x_0,x_1,...,x_n] (x-x_1) (x-x_2) ... (x-x_n) + f[x_1,...,x n] (x-x_2)...(x-x_n) + ... + f[x_{n-1},x_n] (x-x_n) + f[x_n]
    ##
    n = length(xs)-1

    for( k in seq(1, n) ){
        for( i in seq(0, n-k) ){
            ds[i+1] = (ds[(i+1)+1] - ds[i+1])/(xs[(i+k)+1] - xs[i+1])
        }
    }
    return(ds)
}

newton_form_interpolating_polynomial = function(x, y, z_test){
    ##
    ## This is a shortcut code that:
    ## 1) Builds the Newton interpolating polynomial
    ## 2) Evaluates it at the points z_test
    ##

    ##
    ## Get the newton_form coefficients:
    ##
    as = coefficients_of_newton_form(x, y)

    x = rev(x) ## convert into the orderning needed by newton_form_nested_multiplication
    as = rev(as)

    ## Evaluate our interpolating polynomial at some points:
    ##
    f_hat = c()
    for( z in z_test ){
        res = newton_form_nested_multiplication(as, x[-length(x)], z)
        f_hat = c( f_hat, res$as_prime[1] )
    }
    return(f_hat)
}



table = function(xbar, xs, ys, tol=1.e-3, max_steps=20, verbose=FALSE){
    ##
    ## Implements interpolation in a function table. This is the 'table' subprogram in this chapter.
    ##
    ## xbar = single x value we want to evaluate the interpolating polynomial at
    ## xs = [x_0, x_1, ..., x_n]
    ## ys = [f(x_0), f(x_1), ..., f(x_n)]
    ##
    ntable = length(xs) ## the number of points in the function table

    stopifnot(all(diff(xs)>0)) ## xs is assumed to be increasing
    stopifnot(length(xs) ==length(ys)) ## length(xs)==length(ys)
    stopifnot(xbar > xs[1]) ## we must have the point xbar internal to the array: xs[1] < xbar
    stopifnot(xbar < xs[ntable]) ## xbar < xs[length(xs)]

    ## Locate xbar in the x-array such that xs[NEXT-1] < xbar <= xs[NEXT]:
    ##
    NEXT = 1
    while( xs[NEXT] <= xbar ){
        NEXT = NEXT + 1
    }
    stopifnot((xs[NEXT-1] < xbar) & (xbar <= xs[NEXT]))

    ## Use Algorithm 2.4 with the next x_k a1ways the table entry nearest xbar of those not yet used:
    ##
    XK = rep(NA, ntable)
    XK[1] = xs[NEXT]
    NEXTL = NEXT-1
    NEXTR = NEXT+1

    ## A[] holds the divided difference table: f[x_i, ..., x_{k+1}]
    ##
    A = rep(NA, ntable)
    A[1] = ys[NEXT]


    ## This will hold the desired interpolation value p_{k+1}(xbar)
    ##
    TABLE = A[1]

    PSIK = 1

    kp1 = 1
    error = NA
    if( verbose ){ print(sprintf('kp1= %4d; NEXT= %3d; XK[kp1]= %8.3f; TABLE= %8.3g; error= %8.3g', kp1, NEXT, XK[kp1], TABLE, error)) }

    kp1max = min(max_steps, ntable)
    for( kp1 in 2:kp1max ){
        ##
        ## Find the next point to add
        ##
        if(NEXTL == 0){
            NEXT = NEXTR
            NEXTR = NEXTR+1
        } else if(NEXTR > ntable){
            NEXT = NEXTL
            NEXTL = NEXTL-1
        } else if((xbar - xs[NEXTL]) > (xs[NEXTR] - xbar)){
            NEXT = NEXTR
            NEXTR = NEXTR+1
        } else {
            NEXT = NEXTL
            NEXTL = NEXTL-1
        }

        XK[kp1] = xs[NEXT]
        A[kp1] = ys[NEXT]

        for( jj in seq(kp1-1, 1, -1) ){
            A[jj] = (A[jj+1] - A[jj])/(XK[kp1] - XK[jj])
        }

        ## For I=1, ..., KP1 the array A[I] now containts the divided difference of F[X] of order K-1 at XK[I], ..., XK[KP1]
        ##
        PSIK = PSIK*(xbar - XK[kp1-1])
        error = A[1] * PSIK
        if( verbose ){ print(sprintf('kp1= %4d; NEXT= %3d; XK[kp1]= %8.3f; TABLE= %8.3g; error= %8.3g', kp1, NEXT, XK[kp1], TABLE, error)) }
        TABLE = TABLE + error
        if( abs(error) <= tol ){return(TABLE)}
    }
    return(TABLE)
}


test_function_table = function(){
    ##
    ## Test the function 'table' (the one above this one) on some data.
    ##
    y_fn = function(x){ log10(x) } ## these are the functions we will study
    ##y_fn = function(x){ exp(-x) }
    ##y_fn = function(x){ 1/(1+x) }
    ##y_fn = function(x){ 4*x*42 - x + 10 } ## our interpolation should be exact with at least three true (x_i, y_i) points

    ##xs = c(1.0, 1.5, 2.0, 3.0, 3.5, 4.0)
    xs = c(1.0, 2.0, 4.0)

    ys = y_fn(xs)

    x_grid = seq(min(xs)*1.01, max(xs)*.99, length.out=100)
    y_true = y_fn(x_grid)
    y_hat = rep(NA, length(x_grid))
    for( ii in 1:length(x_grid) ){
        y_hat[ii] = table(x_grid[ii], xs, ys, verbose=FALSE)
    }

    plot(xs, ys, type='p', pch=19, col='black', xlab='x', ylab='y(x)')
    lines(x_grid, y_hat, col='red')
    lines(x_grid, y_true, col='green')
    legend('bottomright', c('y_polynomial','y_true'), col=c('red','green'), lty=1)
    grid()
}


make_divided_difference_table = function(xs, ys, n_differences=3){
    ##
    ## xpdf ../../EBook/Elementary_Numerical_Analysis_An_Algorithmic_Approach_3rd_Edition.pdf 55 & 
    ##
    ## In the returned table T we have
    ##
    ## T[i, k] = (f[x_]
    ##
    stopifnot(length(xs)==length(ys)) ## x and y must be the same length
    stopifnot(n_differences>=1)
    n = length(xs)-1
    stopifnot(n_differences<=n)

    n_elements = 2*(n+1) + n_differences*(n+1)
    T = matrix(rep(NA, n_elements), nrow=(n+1), ncol=n_differences+2)
    T[, 1] = xs ## the first column is the x_i for i=0, 1, 2, ..., n
    T[, 2] = ys ## the second column is the y_i for i=0, 1, 2, ..., n
    
    for(jj in 3:(n_differences+2)){
        k = jj-2 ## this is the order of the difference in the jjth column (k=1 => first order difference; k=2 => second order difference etc.)
        for(ii in 1:(n+1-k)){
            ##print(sprintf('book_k= %d; book_i= %d; denom=(x_%d - x_%d)', k, ii-1, ii+k-1, ii-1))
            T[ii, jj] = (T[ii+1, jj-1] - T[ii, jj-1])/(xs[ii+k] - xs[ii])
        }
    }
    
    return(T) 
}


test_make_divided_difference_table = function(){
    ##
    ## some checks that are code is correct
    ##
    
    ## An example from: https://www.geeksforgeeks.org/newtons-divided-difference-interpolation-formula/
    ##
    xs = c(5, 6, 9, 11)
    ys = c(12, 13, 14, 16)
    T = make_divided_difference_table(xs, ys, n_differences=3)
    print(T)

    ## An example 7 from: xpdf /home/wax/EBooks/_Mathematics/_Numerical/Numerical_Mathematics_N_Computing_6th_Edition_Cheney_N_Kincaid.pdf 159 & 
    ##
    xs = c(1, 3/2, 0, 2)
    ys = c(3, 13/4, 3, 5/3)
    T = make_divided_difference_table(xs, ys, n_differences=3)
    print(T)

}