#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <iterator>

using namespace std;

vector< double > intersection2_f( vector< double >& x,vector< double >& y )
{

	sort( x.begin(), x.end() );
	sort( y.begin(), y.end() );

	vector< double > res;
	
	vector< double >::iterator xx = x.begin();
	vector< double >::iterator yy = y.begin(); 

	for( ;
		yy != y.end(); 
		yy++ )
	{
		for(; *xx < *yy && xx != x.end(); xx++ );

		if ( *xx == *yy )
			res.push_back(*yy);
	}

	return res;
}

#ifdef OCTAVE
#include <octave/oct.h>
Matrix convert( vector< double > & stl_v )
{
	int rows = stl_v.size();
	
	if (rows == 0)
		return Matrix( 0, 0 );

	Matrix oct_mx( rows, 1 );
	
	for( int r=0; r<rows; r++)
		oct_mx.elem(r,0) = stl_v[r];

	return oct_mx;
}

vector<double> convert( Matrix & oct_v )
{
	vector<double> v;
	int rows = oct_v.length();
	
	if (rows == 0)
		return v;	
	
	for( int r=0; r<rows; r++)
		v.push_back( oct_v.elem(r,0) );

	return v;
}


DEFUN_DLD( intersection2,args, ,"m = intersection2( x, y )\n subtracts set X from Y: Y \\ X ")
{
	// parsing parameters
	if ( args.length() != 2 )
	{
		print_usage ("intersection2");
		return octave_value(0);
	}

	Matrix oct_x = args(0).matrix_value();
	Matrix oct_y = args(1).matrix_value();
	
	vector< double > x = convert( oct_x );
	vector< double > y = convert( oct_y );
	
	vector< double > res = intersection2_f( x, y );
	Matrix oct_mx = convert( res );
	return octave_value( oct_mx );
}

#else

// for testing purposes

int main()
{
	vector< double > x(4,0), y(5,0), res;
	x[0] = 64;
	x[1] = 8;
	x[2] = 12;
	x[3] = 30;
	
	y[0] = 1;
	y[1] = 30;
	y[2] = 5;
	y[3] = 12;
	y[4] = 20;

		
	cout << endl;
	res = intersection2_f( x, y );
	for( vector<double>::iterator ii=res.begin(); ii != res.end(); ii++ )
		cout << *ii << endl;

}
#endif