// Distributions de Kolmogorov (Scilab) 
// 7 decembre 2011

// Calcul ``exact'' de Marsaglia et al.
// kolmogororovcdf(<ks>, <n>), kolmogorovicdf(<p>, <n>)

// Approximation de Kolmogorov-Dudley de (avec petites am\'eliorations)
// klmcdf(<ks>, <n>), klmicdf(<p>, <n>)

// Distribution de Kolmogorov-Smirnov pour les test unilat\'eraux
// kpmcdf(<ks>, <n>), kpmicdf(<p>, <n>)

// puissance de matrices avec conservation d'exposants

function [V, eV] = mPower(A, eA, n);
	if n == 1 then V = A; eV = eA; return;
	else
		[V, eV] = mPower(A, eA, floor(n/2));
		B = V*V; eB = 2*eV;
		if modulo(n, 2) == 0 then
			V = B; eV = eB; return;
		else
			V = A*B; eV = eA; return;
		end
	end
endfunction

function p = kolmogorovcdf(d, n)
//	local int k, m, eH, eQ; et les param\`etres de boucle i, j, g
//	local double h, s, H, eH, Q, eQ;
// Supprimer les prochaines lignes pour une
// pr\'ecision sup\'erieure \`a 7 d\'ecimales
	s = d*d*n;
	if s > 7.24 | (s > 3.76 & n > 99)
		p = 1-2*exp(-(2.000071+.331/sqrt(n)+1.409/n)*s);
		return;
	end
	k = floor((n*d)+1); m = 2*k-1; h = k-n*d;
	H = ones(m, m); Q = zeros(m, m);
	for i = 1:m;
		for j = 1:m;
			if (i+1 < j) then H(i,j) = 0;
			end
		end
	end
	for i = 1:m;
		// le calcul des puissances $h^n$ pourrait \^etre am\'elior\'e
		// au niveau du temps. La question se pose au niveau
		// de la pr\'ecision...
		H(i,1) = H(i,1)-h^i;
		H(m,i) = H(m, i)-h^(m-i+1);
	end
	if 2*h-1 > 0 then H(m,1) = H(m,1)+(2*h-1)^m; end
	for i = 1:m;
		for j = 1:m;
			if i-j+1 > 0 then
				for g = 1:(i-j+1);
					H(i,j) = H(i,j)/g;
				end
			end
		end
	end
	eH = 0; [Q, eQ] = mPower(H, eH, n);
	p = Q(k,k);
	for i = 1:n;
		p = p*i/n;
		if p < 1e-140 then p = p*1e140; eQ = eQ-140; end
	end
	p = p*(10^eQ);
endfunction

function x = kolmogorovicdf(p, n)
//	global accuracy;
	accuracy = 1e-12;
//	local double a, b;
	if p > 1-accuracy then x = 1; return;
	elseif p < accuracy then x = 0; return;
	else
		a = 0; b = 1;
		while b-a >= accuracy;
			if kolmogorovcdf((a+b)/2, n) < p then
				a = (a+b)/2;
			else b = (a+b)/2;
			end
		end
		x = (a+b)/2; return;
	end
endfunction

// test\'e, impl\'ementation asymptotique (Dudley)

//function p = klmcdf(x, n)
//	if x <= 0 then p = 0;
//	elseif x >= 1 then p = 1;
//	else
//		// local z p term k pm;
//		// global accuracy;
//		accuracy = 1e-12;
//		z = -2*(sqrt(n)+0.12+0.11/sqrt(n))^2*x*x;
//		p = 1.0; k = 1; pm = -1; term = 1;
//		while abs(term) > accuracy & k < 100;
//			term = 2*pm*exp(z*k*k);
//			p = p+term; pm = -pm; k = k+1;
//		end
//		p = max(min(p, 1), 0);
//	end
//endfunction
//
//function x = klmicdf(p, n)
////	global accuracy;
//	accuracy = 1e-12;
////	local double a, b;
//	if p > 1-accuracy then x = 1; return;
//	elseif p < accuracy then x = 0; return;
//	else
//		a = 0; b = 1;
//		while b-a >= accuracy;
//			if klmcdf((a+b)/2, n) < p then
//				a = (a+b)/2;
//			else b = (a+b)/2;
//			end
//		end
//		x = (a+b)/2; return;
//	end
//endfunction

function lg = lnGamma(x);
  lg = log(x+4.5)*(x-.5)-x-4.5...
  +log(2.50662827465*(1+76.18009173/x...
      -86.50532033/(x+1)...
      +24.01409822/(x+2)...
      -1.231739516/(x+3)...
      +0.00120858003/(x+4)...
      -0.00000536382/(x+5)));
endfunction

function p = klmcdf(x, n);
	if x <= 1/(2*n) then p = 0;
	elseif x >= 1 then p = 1;
	elseif x <= 1/n then p = 1; for i = 1:n; p = p*i*(2*x-1/n); end
	elseif x >= 1-1/n then p = 1-2*((1-x)**n);
	else
		// local z p term k pm;
		// global accuracy;
		accuracy = 1e-12;
		z = -2*(sqrt(n)+0.12+0.11/sqrt(n))^2*x*x;
		p = 1.0; k = 1; pm = -1; term = 1;
		while abs(term) > accuracy & k < 100;
			term = 2*pm*exp(z*k*k);
			p = p+term; pm = -pm; k = k+1;
		end
		p = max(min(p, 1), 0);
	end
endfunction

function x = klmicdf(p, n);
	 // global accuracy;
	accuracy = 1e-12;
	if p <= accuracy then x = 1/(2*n);
	elseif p >= 1-accuracy then x = 1;
	elseif log(p) <= lnGamma(n+1)-n*log(n) then
		// x = 0.5*(1/n+(p/factorial(n))**(1/n));
		x = 0.5*(1/n+exp((log(p)-lnGamma(n+1))/n));
	elseif p >= 1-2/(n^n) then
		x = 1-exp(log((1-p)/2)/n);
	else
		a = 0; b = 1;
		while b-a >= accuracy;
			if klmcdf((a+b)/2, n) < p then
				a = (a+b)/2;
			else b = (a+b)/2;
			end
		end
		x = (a+b)/2; return;
	end
endfunction

function p = kpmcdf(x, n);
	if x <= 0 then  p = 0;
	elseif x >= 1 then  p = 1;
	else p = 0;
		for k = n:-1:ceil(n*x+1e-4);
			p = p+exp((lnGamma(n+1)-lnGamma(k+1)-lnGamma(n-k+1))...
			+k*log(k/n-x)+(n-k-1)*log(x+1-k/n));
		end
		p = 1-x*p;
	end
endfunction

function x = kpmicdf(p, n);
	 // global accuracy;
	accuracy = 1e-12;
	if p > 1-accuracy then x = 1;
	elseif p < accuracy then  x = 0;
	else // local lb, ub;
	      lb = 0; ub = 1;
		while ub-lb >= 2*accuracy;
			if kpmcdf((lb+ub)/2, n) < p then lb = (lb+ub)/2;
			else ub = (lb+ub)/2; end
		end
		x = (lb+ub)/2;
	end
endfunction

// petit bonus

function y = choose(n, k);
	 y = 1;
	 for ell = 1:k; y =y*(n-ell+1)/ell; end
	 y = round(y);
endfunction