
function r=binom(k,n)
    if k > n/2
        k=n-k
    end
    p1=1
    for i=2:k
        p1=p1*i
    end
    p2=1
    for i=n-k+1:n
        p2=p2*i
    end
    r=p2/p1
endfunction

// il faut au plus 2*m multiplications avec m=min(k,n-k)
// par une méthode naive en utilisant directement la définition avec la factorielle il aurait fallu n+2*k multiplications

function r=binomrec(k,n)
    if k==0 | k==n
        r=1
    else
        r=binomrec(k-1,n-1)+binomrec(k,n-1)
    end
endfunction

// On utilise dans bernstein la fonction binom, plus performante que binomrec

function R=bernstein(n,f,T)
    R=[]
    for i=1:length(T)
        s=0
        for k=0:n
            s=s+binom(k,n)*T(i)^k*(1-T(i))^(n-k)*f(k/n)
        end
        R=[R,s]
    end
endfunction


function y=f(x)
    y=cos(10*x)
endfunction


function R=bernsteinOP(n,f,T)
    R=[]
    // On calcule d'abord le tableau des coefficients binomiaux
    CBin=ones(n+1,n+1)
    for i=3:n+1
        for k=2:(i-1)
            CBin(k,i)=CBin(k-1,i-1)+CBin(k,i-1)
        end
    end
    
    for i=1:length(T)
        // On calcule ensuite, pour chaque valeur de la liste T, la liste des T(i)^k*(1-T(i))^n-k
        if T(i)==1
            LPuiss=[zeros(1,n),1]
        else
            LPuiss=[(1-T(i))^n]
            for k=1:n
                LPuiss=[LPuiss,LPuiss($)*T(i)/(1-T(i))]
            end
        end
        
        
        // Finalement on calcule la valeur du polynome de bernstein
        s=0
        for k=1:n+1
            s=s+CBin(k,n+1)*LPuiss(k)*f((k-1)/n)
        end
        R=[R,s]
    end
endfunction
