Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

i_partition_next2 ( n, npart, a, ...
function [ npart_new, a_new, mult_new, more_new ] = i_partition_next2 ( n, npart, a, ...
  mult, more )

%% I_PARTITION_NEXT2 computes the partitions of the integer N one at a time.
%
%  Discussion:
%
%    Unlike compositions, order is not important in a partition.  Thus
%    the sequences 3+2+1 and 1+2+3 represent distinct compositions, but
%    not distinct partitions.  Also 0 is never returned as one of the
%    elements of the partition.
%
%    Initialize the program by calling with MORE = FALSE.  On an initialization
%    call, the input values of A, MULT and NPART are not needed.  Thereafter,
%    they should be set to the output values of A_NEW, MULT_NEW and NPART_NEW
%    from the previous call.
%
%  Examples:
%
%    Sample partitions of 6 include:
%
%      6 = 4+1+1 = 3+2+1 = 2+2+2
%
%  Modified:
%
%    20 July 2004
%
%  Parameters:
%
%    Input, integer N, the integer whose partitions are desired.
%
%    Input, integer NPART, the output value of NPART_NEW on the previous call.
%
%    Input, integer A(N), the output value of A_NEW on the previous call.
%
%    Input, integer MULT(N), the output value of MULT_NEW on the previous call.
%
%    Input, logical MORE, is FALSE on the first call, which causes 
%    initialization.  Thereafter, it should be TRUE.
%
%    Output, integer NPART_NEW, the number of distinct, nonzero parts in the
%    output partition.
%
%    Output, integer A_NEW(N).  A_NEW(1:NPART_NEW) the distinct parts
%    of the partition.
%
%    Output, integer MULT_NEW(1:NPART), the multiplicity of the parts.
%
%    Output, logical MORE_NEW is TRUE if there are more partitions available.
%
  more_new = more;
  
  if ( ~more_new )
    npart_new = 1;
    a_new(npart_new) = n;
    mult_new(npart_new) = 1; 
    more_new = ( mult_new(npart_new) ~= n );
    return
  end
 
  npart_new = npart;
  a_new(1:npart_new) = a(1:npart_new);
  mult_new(1:npart_new) = mult(1:npart_new);
  isum = 1;
 
  if ( a_new(npart) <= 1 )
    isum = mult_new(npart_new) + 1;
    npart_new = npart_new - 1;
  end

  iff = a_new(npart_new) - 1;

  if ( mult_new(npart_new) ~= 1 )
    mult_new(npart_new) = mult_new(npart_new) - 1;
    npart_new = npart_new + 1;
  end

  a_new(npart_new) = iff;
  mult_new(npart_new) = 1 + floor ( isum / iff );
  suma = sum(a_new.*mult_new);
  if(suma > n)
      a_hold = a_new(1:end-1);
      mult_hold = mult_new(1:end-1);
      clear a_new mult_new;
      a_new = a_hold;
      mult_new = mult_hold;
  end
  is = mod ( isum, iff );

  if ( 0 < is )
    npart_new = npart_new + 1;
    a_new(npart_new) = is;
    mult_new(npart_new) = 1;
  end
%
%  There are more partitions, as long as we haven't just computed
%  the last one, which is N copies of 1.
%
  more_new = ( mult_new(npart_new) ~= n );

Contact us