function table = derangement_table ( hi, trials )

%*****************************************************************************80
%
  for n = 1 : hi
    table(n) = derangement_prob ( n, trials );
  end

  return
end
function prob = derangement_prob ( n, trials )

%*****************************************************************************80
%
  derangement_num = 0;

  for t = 1 : trials
    p = permutation_random ( n );
    if ( permutation_is_derangement ( n, p ) )
      derangement_num = derangement_num + 1;
    end
  end

  prob = derangement_num / trials;

  return
end
function p = permutation_random ( n )

%*****************************************************************************80
%
  p = ( 1 : n );

  for i = 1 : n - 1

    k = i4_random ( i, n );

    p1 = p(i);
    p(i) = p(k);
    p(k) = p1;

  end

  return
end
function value = i4_random ( lo, hi )

%*****************************************************************************80
%
  r = ( hi + 1 - lo ) * rand ( );
  value = lo + floor ( r );

  return
end
function value = permutation_is_derangement ( n, p )

%*****************************************************************************80
%
  if ( any ( p(1:n) == (1:n) ) )
    value = 0;
  else
    value = 1;
  end

  return
end
