/******************************************
**    Combinatorial Matrix    **
**
**    Author: Gerardo Gomez-Ruano
**    Date: July 2002, Mexico City
**    mailto: ggr7@hotmail.com
**
**    This code is provided gratis without any guarantees or
**    warrantees. Proper attributions for the use of this code
**    are welcomed. Proprietary modifications of the code
**    are not permitted. I hope this program is useful to you.
**
**    Function: compute all binary combination vectors, that is,
**    compute all possible vectors of dimension  1X n such
**    that they have k components with a 1 (one) and
**    n-k components with a 0 (zero).
**
**    Inputs: n -- the length of the vectors (ones and zeros)
**            k -- the number of ones in the vectors (ones)
**
**    Output: nCk vectors of dimension 1 x n
**            where each row is a combination
**
**    e.g.1      input: n=4, k=2
**                 output:
**                               1 1 0 0
**                               1 0 1 0
**                               1 0 0 1
**                               0 1 1 0
**                               0 1 0 1
**                               0 0 1 1
**
**     e.g.2     input: n=3, k=1
**                  output:
**                                1 0 0
**                                0 1 0
**                                0 0 1
**
** Comments: The program is in a loose form, so that any
** new operation can be introduced for each vector, like:
** multiply each possible combination vector i with the constant
** vector z and if the product is equal to w then save it in a matrix A.
**         It is worth saying that this program, unlike others on
**  combinations, cannot be found in other mathematical
**  applications, like Mathematica.
**  It is a useful application for identifying individuals with group data
**  like: you have a list of 10 individuals, with their respective
**  incomes Ii (I=1,2,...,10),
** and you know the number of members of family X ( say, 3) and
** Familiy X's aggregate income Ix. Who of all the individuals are
** the members of  Family X?
** For this, you need any combination of three individuals out of ten
** whose aggregate income is equal to Familiy X's. Thus, you
**  multiply the constant vector of individual incomes by each
**  possible combination vector and check if the result is equal to
**  Family X's Income.
**
**  The estimated time (at the start of the program) is a proxy for
**  my computer, it's not
**  necessarily useful to you, though it could be.
******************************************/

new;
cls;
format /rdn 10,1;
print " Combinatorial Matrix Version 3.0, written by Gerardo Gomez-Ruano
This program receives arguments n and k (with n>=k>=0) and constructs the 
n!/[k!(n-k)!] line vectors of dimension 1 X n, such that each one is a 
distinct possible combination of k  'one' entries and n-k 'zero' entries.
Type n";
n=con(1,1);
Print "Type k, with k<=n";
request:
k=con(1,1);
if k>n ;
print "it is imperative that k<=n, please introduce k again" ;
goto request ;
endif;


vectores= n!/(k!*(n-k)!) ;
print " number of final vectors =";
print vectores;

vectoresp=2^n;
print "number of potential vectors =";
print vectoresp;

tiempo=(1.597991212 + 0.002666199874*vectores + 0.0001319554622*vectoresp - 
0.2175601307*k*(n-k))/60;
print "estimated time (minutes) :";
print tiempo;

tiempo2=(-5.813347569 + 0.001903439275*vectores + 0.0001334369273*vectoresp 
+ 1.981826546e-09*(vectores*vectoresp))/60;
print tiempo2;


print "ready to go on? Press any key";
waitc;
t0=hsec;
/* Begin vector construction */

uno={1};
cero={0};

	for i (1, vectoresp, 1 );

	/* computation of w for potential vector i */
	h=i;
	w=0;

		for v (1, n, 1 );
		z=2^(n-v);
			if h<=z;
			w=w+1;
			else;
			h=h-z;
			endif;

		endfor;

	/* End of computation of w */

	/* validity proof for the potential vector i and construction if valid */

		if k==w;
		vector={};
		h=i;

			for v (1, n, 1 );
			z=2^(n-v);
				if h<=z;
				vector=vector~uno;
				else;
				vector=vector~cero;
				h=h-z;
				endif;

			endfor;
			print vector;
		endif;

	endfor;

print "Minutes elapsed  " ftos((hsec-t0)/6000,"%*.*lf",10,5);
print "miliseconds per final vector          " 
ftos((hsec-t0)/(100*vectores),"%*.*lf",10,5);
print "miliseconds per potential vector " 
ftos((hsec-t0)/(100*vectoresp),"%*.*lf",10,5);