function [y]=coefbin1(n)
// fournit la liste des coefficients binomiaux C_n^k
// pour k=0 jusqu'a n
y=zeros(0:n);
for k=0:n,
	y(k+1)=prod(1:n)/(prod(1:k)*prod(1:(n-k)));
end
endfunction

function [y]=coefbin2(n)
// version avec factorielles, mais sans boucle for
B=[1 cumprod(1:n)]; // ajout de 0!
C=B((n+1):-1:1);    // retournement k! --> (n-k)!
D=B.*C;             // denominateurs k!(n-k)!
y=prod(1:n)./D
endfunction

function [y]=coefbin3(n)
// 
y(1)=1;
y(2:(n+1))=cumprod(n:-1:1)./cumprod(1:n);
endfunction
// 
// cela marche bien pour n <= 170
// A partir de n=171, on a des pbs pour les grandes valeurs de k
// z=coefbin3(171);z(167:170)
// ans  =   1.0E+08 *  !   11.486197    0.3438981   Inf   Inf !
// On peut ameliorer sensiblement en utilisant la fdr de la loi binomiale
// Bin(n,1/2) qui repose sur la fonction Beta


function [y]=coefbin4(n,k)
// calcul de C_n^k a l'aide de la fdr binomiale
if (k>floor(k))|(n>floor(n)) then 
	warning('les parametres doivent etre entiers');
elseif k==0 then y=1;
else y=2.^n.*(cdfbin("PQ",k,n,0.5,0.5)-cdfbin("PQ",k-1,n,0.5,0.5));
end
endfunction
//
// Attention, la fonction cdfbin donne des valeurs inattendues pour
// k non entier, d'ou le test rajoute sur k (et sur n).
//
// Cette fonction est limitee a n<=1023 a cause des puissances de 2
// D'autre part pour les grandes valeurs de k, le resultat n'est pas fiable
// comparer coefbin4(1000,650) et coefbin4(1000,350) pour s'en convaincre
//
// On peut encore ameliorer en se ramenant a k<=n/2
//
function [y]=coefbin5(n,k)
// calcul de C_n^k a l'aide de la fdr binomiale
if k>ceil(n./2) then k=n-k;
end
y=coefbin4(n,k);
endfunction

// Finalement, la  methode basee sur le triangle de Pascal semble preferable
// elle ne repose que sur des additions. Le code tres court ci dessous
// evite de consommer trop de place memoire en stockant tout le tableau des
// C_n^k. On ne garde que la ligne courante.
//
function [y]=coefbinA(n)
// "A" comme additive
// fournit le vecteur ligne des C_n^k, k=0:n
C=[1]; // initialisation n=0
for i=1:n, C=[C 0]+[0 C];end
y=C
endfunction
//
// Limitations: on peut aller jusqu'a n=1029
// pour le voir, taper
// y=coefbinA(1029);max(y)
// puis
// y=coefbinA(1030);max(y)
