function ret = pairwise_mean (x, dim)
  siz = size (x);
  if (nargin < 2)
    dim = find (siz > 1,1);
    if isequal (dim,[])
      dim = 1;
    endif
  endif
  n = siz(dim);
  nn = n;
  siz(dim) = 1;
  L = prod (size (x)(1:dim-1));
  inds = ((1:L)'+L*n*(0:prod (siz)/L-1))(:);
  k = 0;
  while n > 1
#    x(inds+L*(0:n/2-1)) += x(inds+L*(ceil (n/2):n-1));
#    x(inds+L*(0:(n-1)/2)) /= 2;
    x(inds+L*(0:n/2-1)) += (x(inds+L*(ceil (n/2):n-1)) - x(inds+L*(0:n/2-1))) / 2;
    x(inds+L*floor(n/2)) /= 2;
    n = ceil (n/2);
    ++k;
  end
  ret = reshape (x(inds),siz) * (2^k / nn);
%!test
%! x = randn (10, 11, 12);
%! assert (pairwise_mean (x), mean (x), 2*eps)
%!test
%! x = randn (10, 11, 12);
%! assert (pairwise_mean (x, 1), mean (x, 1), 2*eps)
%!test
%! x = randn (10, 11, 12);
%! assert (pairwise_mean (x, 2), mean (x, 2), 2*eps)
%!test
%! x = randn (10, 11, 12);
%! assert (pairwise_mean (x, 3), mean (x, 3), 2*eps)
%!test
%! x = [1.7e308, 1.5e308];
%! assert (pairwise_mean (x), 1.6e308, 2*eps)
%!test
%! rand ("state",3); x = rand (10, 1) * 1e-323;
%! assert (pairwise_mean (x) > 0, true)
