Polynom <- function(x,a){
    n <- length(a)
    y <- a[n]; for (i in 1:(n-1)){y <- x*y+a[n-i]}
    y
}
LagrangeIP <- function(x,y){
    n <- length(x)
    if(!(length(y) == n)){stop("Nesouhlasi delky vektoru")}
    A <- array(1,c(n,n))
    for (i in 1:n){
        for (j in 2:n){A[i,j] <- x[i]*A[i,j-1]}
    }
    solve(A,y)
}

Spline2 <- function(x,y,xout=NULL){
    n <- length(x)
    if(!(length(y) == n)){stop("Nesouhlasí délky vektorů")}
    if(min(x[2:n]-x[1:(n-1)]) < 0){stop("Vektor x není uspořádán")}
    A <- array(0,c(3*n-3,3*n-3))
    b <- array(0,3*(n-1))
    A[1,1] <- 1
    for(k in 2:n){        
        i <- 3*k-4
        j <- 3*k-5
        A[i,j]  <- x[k-1]^2; A[i,j+1] <- x[k-1]; A[i,j+2] <- 1
        b[i] <- y[k-1]
        A[i+1,j] <- x[k]^2; A[i+1,j+1] <- x[k]; A[i+1,j+2] <- 1
        b[i+1] <- y[k]
    }
    for(k in 2:(n-1)){
        i <- 3*k-2
        j <- 3*k-5
        A[i,j] <- 2*x[k]; A[i,j+1] <- 1; A[i,j+3] <- -2*x[k]; A[i,j+4] <- -1
    }
    S <- t(matrix(solve(A,b),c(3,n)))
    if(length(xout) == 0){xout <- seq(min(x),max(x),length.out=100)}
    yout <- array(NA,length(xout))
    for(k in 2:n){
        ind <- which((xout >= x[k-1]) & (xout <= x[k]))
        yout[ind] <- Polynom(xout[ind],S[k-1,3:1])
    }
    list(xout=xout,yout=yout,spline=S)
}


Spline3 <- function(x,y,xout=NULL){
    n <- length(x)
    if(!(length(y) == n)){stop("Nesouhlasí délky vektorů")}
    if(min(x[2:n]-x[1:(n-1)]) < 0){stop("Vektor x není uspořádán")}
    A <- array(0,c(4*(n-1),4*(n-1)))
    b <- array(0,4*(n-1))
    A[1,1] <- 6*x[1]; A[1,2] <- 2
    A[4*(n-1),4*n-7] <- 6*x[n]; A[4*(n-1),4*n-6] <- 2
    for(k in 2:n){        
        i <- 4*k-6
        j <- 4*k-7
        A[i,j]  <- x[k-1]^3; A[i,j+1] <- x[k-1]^2; A[i,j+2] <- x[k-1]
        A[i,j+3] <- 1; b[i] <- y[k-1]
        A[i+1,j] <- x[k]^3; A[i+1,j+1] <- x[k]^2; A[i+1,j+2] <- x[k]
        A[i+1,j+3] <- 1; b[i+1] <- y[k]
    }
    for(k in 2:(n-1)){
        i <- 4*(k-1)
        j <- 4*k-7
        A[i,j] <- 3*x[k]^2; A[i,j+1] <- 2*x[k]; A[i,j+2] <- 1
        A[i,j+4] <- -3*x[k]^2; A[i,j+5] <- -2*x[k]; A[i,j+6] <- -1
        A[i+1,j] <- 6*x[k]; A[i+1,j+1] <- 2
        A[i+1,j+4] <- -6*x[k]; A[i+1,j+5] <- -2
    }
    S <- t(matrix(solve(A,b),c(4,n)))
    if(length(xout) == 0){xout <- seq(min(x),max(x),length.out=100)}
    yout <- array(NA,length(xout))
    for(k in 2:n){
        ind <- which((xout >= x[k-1]) & (xout <= x[k]))
        yout[ind] <- Polynom(xout[ind],S[k-1,4:1])
    }
    list(xout=xout,yout=yout,spline=S)
}

################################################################################

N <- 5; x <- seq(1,N); y <- runif(N,0,1) # náhodné hodnoty


x <- c(-pi/2,-pi/3,-pi/4,-pi/6,0,pi/6,pi/4,pi/3,pi/2)
y <- c(-1,-sqrt(3)/2,-sqrt(2)/2,-0.5,0,0.5,sqrt(2)/2,sqrt(3)/2,1)

y <- c(0,0.5,sqrt(2)/2,sqrt(3)/2,1,sqrt(3)/2,sqrt(2)/2,0.5,0)

################################################################################


rX <- max(x)-min(x); Xlim <- c(min(x)-0.25*rX,max(x)+0.25*rX)
rY <- max(y)-min(y); Ylim <- c(min(y)-0.5*rY,max(y)+0.5*rY)
plot(x,y,type="p",bty="l",pch=19,col="red",
     xlab="x",ylab="y",xlim=Xlim,ylim=Ylim)
 
a <- LagrangeIP(x,y)
xi <- seq(Xlim[1],Xlim[2],length.out=500)
points(xi,Polynom(xi,a),type="l")

################################################################################

Spl2 <- Spline2(x,y)
points(Spl2$xout,Spl2$yout,type="l",col="green")

Spl3 <- Spline3(x,y)
points(Spl3$xout,Spl3$yout,type="l",col="blue")