accuracy = 1E-12;

function [pmin, pmax] = proportionCLI(k, n, alpha);
  if k == 0 then pmin = 0; pmax = 1-(alpha/2)^(1/n);
  elseif k == n then pmin = (alpha/2)^(1/n); pmax = 1;
  elseif n < 50 then
    global accuracy;
    // Clopper-Pearson
    // local p, lb, ub;
    lb = 0; ub = k/n;
    while ub-lb > accuracy;
      p = .5*(lb+ub);
      if 1-cdfbin("PQ", k-1, n, p, 1-p) < alpha/2 then lb = p;
      else ub = p;
      end
    end
    pmin = 0.5*(lb+ub);
    lb = k/n; ub = 1;
    while ub-lb > accuracy;
      p = 0.5*(lb+ub);
      if cdfbin("PQ", k, n, p, 1-p) < alpha/2 then ub = p;
      else lb = p;
      end
    end
    pmax = 0.5*(lb+ub);
//  elseif k < 10 then [pmin, pmax] = poissonCLI(k, alpha)/n;
//  elseif n-k < 10 then [p, q] = poissonCLI(n-k, alpha)/n;
//    pmin = 1-q; pmax = 1-p;
  else
    // Wilson
    // local p, z, t;
    p = k/n;
    z = cdfnor("X", 0, 1, 1-alpha/2, alpha/2);
    z = z*z/n; t = sqrt(z*(z/4+p*(1-p)));
    pmin = (p+z/2-t)/(1+z); pmax = (p+z/2+t)/(1+z);
  end
endfunction

function p = poissoncdf(x, lambda);
  if x < 0 then p = 0;
  elseif lambda > 0 then
    // local k;
    p = 1.0;
    for k = floor(x):-1:1; p = 1+p/k*lambda; end
    p = p*exp(-lambda);
  else p = 1;
  end
endfunction

function [lmin, lmax] = poissonCLI(k, alpha);
  if k == 0 then lmin = 0; lmax = -log(alpha/2);
  else
    global accuracy;
    // local lb, ub;
    // lowerbound
    lb = 0.0; ub = 1.0*k;// mauvaise borne supperieure
    while ub-lb > accuracy;
      if 1-poissoncdf(k-1, (lb+ub)/2) < alpha/2 then lb = (lb+ub)/2;
      else ub = (lb+ub)/2;
      end
    end
    lmin = (lb+ub)/2;
    // upperbound
    lb = 0.0; ub = 1.0*k;// mauvaise borne superieure
    while poissoncdf(k, ub) > alpha/2; ub = 2*ub; end
    while ub-lb > accuracy
      if poissoncdf(k, (lb+ub)/2) > alpha/2 then lb = (lb+ub)/2;
      else ub = (lb+ub)/2;
      end
    end
    lmax = (lb+ub)/2;
  end
endfunction
