Polynom <- function(x,a){
    n <- length(a)
    y <- a[n]
    if (n>1){for (i in 1:(n-1)){y <- x*y+a[n-i]}}
    else{y <- array(y,length(x))}    
    y
}

TPexp <- function(n){ # McLaurin polynomial - exponential function
    a <- array(1,n+1)
    if (n > 0){for (i in 1:n){a[i+1] <- a[i]/i}}
    a
}

TPsin <- function(n){ # McLaurin polynomial - sine function
    a <- array(0,n+1)
    if (n > 0){
        a[2] <- 1
        if (n > 2){
            for (i in 2:((n+1)/2)){a[2*i] <- -a[2*(i-1)]/((2*i-1)*(2*i-2))}}}
    a
}

TPcos <- function(n){ # McLaurin polynomial - cosine function
    a <- array(0,n+1); a[1] <- 1
    if (n > 1){
        for (i in 1:(n/2)){a[2*i+1] <- -a[2*i-1]/(2*i*(2*i-1))}}
    a
}

log1PX <- function(x){log(1+x)} # function ln(1+x)
TPlog <- function(n){ # McLaurin polynomial - ln(1+x)
    a <- array(0,n+1)
    if (n > 1){
        a[2] <- 1
        for (i in 2:n){a[i+1] <- -a[i]*(i-1)/i}}
    else {if (n == 1){a <- c(0,1)}}
    a
}

XP1naAlpha <- function(x){(1+x)^alpha} # function (1+x)^alpha
TPXP1naAlpha <- function(n){ # McLaurin polynomial - (1+x)^alpha
    a <- array(1,n+1)
    if (n > 0){
        a[2] <- alpha
        if (n > 1){for (i in 2:n){a[i+1] <- a[i]*(alpha-i+1)/i}}}
    a
}

    

DisplayTP <- function( # result plot
    f,TP,              # name of function and its McLaurin polynomial
    n=4,               # degree of polznomial 
    xmin=-2,xmax=2     # range of independent value
    ){
    x <- seq(xmin,xmax,length.out=200)
    y <- f(x)
    TPy <- Polynom(x,TP(n))
    ylm <- c(min(y),max(y))
    err <- max(abs(y-TPy))
    plot(x,f(x),type="l",bty="l",xlab="x",ylab="y",lwd=2,col="red",
         ylim=ylm)
    points(x,Polynom(x,TP(n)),type="l",lwd=3)
    text(x[1],ylm[2],sprintf("      error = %g",err),pos=4)
}