/*
 * examples/samplebb.C
 *
 * Copyright (C) 2005, 2010 D Saunders
 * ========LICENCE========
 * This file is part of the library LinBox.
 *
 * LinBox is free software: you can redistribute it and/or modify
 * it under the terms of the  GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 * ========LICENCE========
 */

/** \file examples/samplebb.C
 * @example  examples/samplebb.C
 * \ingroup examples
 * \brief generate an example matrix with specified frobenius form.
 *
 * samplebb takes options and any number of argument triples denoting companion
 * matrix blocks.
 * For example, the call "samplebb -r 7 2 3 a3 1 1" generates a sparsely
 * randomized matrix (because of the '-r' option) matrix which is similar
 * (because of the two triples '7 2 3' and 'a2 1 1') to a direct sum of 3
 * companion matrices for (x-7)^2 plus one companion matrix for x^3 + x + 2, the
 * polynomial denoted by 'a3'.
 *
 *  In general, in the first position of each triple 'aK' denotes the polynomial
 *  x^k + x + K-1 and a number n denotes the polynomial x-n.  The second number
 *  in the triple specifies a power of the polynomial and the third specifies how
 *  many companion matrix blocks for that power of that polynomial.
 *
 *  Possible options are
 *  -r lightly randomized similarity transform, matrix remains sparse.
 *  -R fully randomized similarity transform, matrix becomes dense.
 *
 *  The matrix is written to standard out in SMS format (triples).
 *
 *  For some other examples:
 *  "samplebb  1 1 2  2 1 2  4 1 2  12 1 1  0 1 1" is a 8 by 8 diagonal matrix in smith form,
 *  diag(1,1,2,2,4,4,12,0)
 *
 */

#include <iostream>
#include <string>
#include <list>
#include <vector>
#include <linbox/blackbox/direct-sum.h>
#include <linbox/blackbox/companion.h>
#include <linbox/algorithms/matrix-hom.h>
#include <linbox/ring/ntl.h>
#include <NTL/ZZX.h>
#include <linbox/vector/blas-vector.h>

using std::string;
using std::list;
using LinBox::Companion;
using LinBox::DirectSum;
using LinBox::DenseMatrix;
using LinBox::NTL_ZZ;
using NTL::ZZX;


void stripOptions(int& acp, char* avp[], string& opts, const int ac, char** av)
{
	acp = 0;
	for (int i = 1; i < ac; ++i)
	{
		//std::cout << av[i] << " ";
		if (av[i][0] == '-') {
			for (const char* j = av[i]+1; *j != 0; ++j)
				opts.push_back(*j);
		}

		else {
			avp[acp] = av[i];
			++acp;
		}
	}
}

template <class List, class Ring>
void augmentBB(List& L, char* code, int e, int k, const Ring& R)
{
	typedef typename Ring::Element Int;
	Int a;
	ZZX p;

	// build poly p

	if ( *code != 'a')  // build linear poly
	{
		R.init(a, -atoi(code));
		p += ZZX(0, a);
		p += ZZX(1, R.one);
	}
	else // build long poly
	{
		int n = atoi(code+1);
		R.init(a, n-1);
		p += ZZX(n, R.one);
		p += ZZX(1, R.one);
		p += ZZX(0, a);
	}

	//std::cout << "(code, e, k) =(" << code << ", " << e << ", " << k << ")" << std::endl;
	//std::cout << "Correspoding poly: " << p << std::endl;
	// compute q =  p^e
	ZZX q(0, R.one);
	for(int i = 0; i < e; ++i) q *= p;
	//std::cout <<"Polynomial: " << q << std::endl;

	LinBox::DenseVector<Ring>  v(R,deg(q)+1);
	for (int i = 0; i < v.size(); ++i) v[i] = coeff(q, i);

	// companion matrix of q
	Companion<Ring>* C = new Companion<Ring>(R, v);
	for(int i = 0; i < k; ++i) L.push_back(C);

}

template < class Ring >
void scramble(DenseMatrix<Ring>& M)
{

	Ring R = M.field();

	int N,n = M.rowdim(); // number of random basic row and col ops.
	N = 2*n;

	for (int k = 0; k < N; ++k) {

		int i = rand()%M.rowdim();

		int j = rand()%M.coldim();

		if (i == j) continue;

		// M*i += alpha M*j and Mj* -= alpha Mi*

		typename Ring::Element alpha, beta, x;
		R.init(alpha, rand()%5 - 2);
		R.neg(beta, alpha);

		for (size_t l = 0; l < M.rowdim(); ++l) if (!R.isZero(alpha)) {
			R.mul(x, alpha, M[l][j]);
			R.addin(M[l][i], x);
		}


		for (size_t l = 0; l < M.rowdim(); ++l) if (!R.isZero(alpha)) {
			R.mul(x, beta, M[i][l]);
			R.addin(M[j][l], x);
		}
	}

	/*
	   std::ofstream out("matrix", std::ios::out);

	//M. write(std::cout);

	out << n << " " << n << "\n";

	for (int i = 0; i < n; ++ i) {

	for ( int j = 0; j < n; ++ j) {

	R. write(out, M[i][j]);

	out << " ";
	}

	out << "\n";

	}
	*/

}

template <class Matrix>
void printMatrix (const Matrix& A)
{
	int m = A. rowdim();
	int n = A. coldim();
	typedef typename Matrix::Field Ring;
	typedef typename Ring::Element Element;
	const Ring &r = A.field();
	LinBox::DenseVector<Ring> x(r,m), y(r,n);

	std::cout << m << " " << n <<  " M" << std::endl;
	typename LinBox::DenseVector<Ring>::iterator y_p;
	for (int i = 0; i < m; ++ i) {
		r. assign (x[i], r.one);
		A. applyTranspose(y, x);
		for(y_p = y. begin(); y_p != y. end(); ++ y_p)
			if (! r.isZero(*y_p))
				std::cout << i+1 << " " << y_p - y.begin() + 1 << " " << *y_p << std::endl;
		r. assign (x[i], r.zero);
	}
	std::cout << "0 0 0" << std::endl;
}


int main(int ac, char* av[])
{
	if (ac < 2)
	{	std::cout << "usage: " << av[0] <<
		" options block-groups." << std::endl;
		std::cout << av[0] << " -r 1 2 3 a4 1 1" << std::endl;
		std::cout <<
		"for lightly randomized matrix similar to direct sum of 3 copies of companion " << std::endl
		<< "matrix of (x-1)^2 and one copy of companion matrix of (x^4 + x + 3)^1." << std::endl;
	}

	typedef NTL_ZZ Ring;
	Ring Z;
	typedef Companion<Ring> BB;
	int acp; char* avp[ac];
	string opts;
	stripOptions(acp, avp, opts, ac, av);
	//std::cout << "number of triples: " << acp << std::endl;
	//for (int i = 0; i < acp; ++ i)
	//	std::cout << avp[i];
	//std::cout << std::endl;
	//std::cout << "Begin to ....\n";
	list<BB*> L;

	for (int i = 0; i < acp; i += 3)
		augmentBB(L, avp[i], atoi(avp[i+1]), atoi(avp[i+2]), Z);

	DirectSum<BB> A(L);
	//std::cout <<"Option: " << opts.c_str() << std::endl;

	if (opts.size() >= 1)
	{	if (opts[0] == 'r')
		{
			// into sparse matrix, then 3n row ops with corresponding col ops
			DenseMatrix<Ring> B(Z,A.rowdim(), A.coldim());
			//MatrixDomain<Ring> MD(Z);
			LinBox::MatrixHom::map (B, A);

			scramble(B);
			printMatrix(B);
			// delete B;
		}

		if (opts[0] == 'R') { throw LinBox::NotImplementedYet(); }
		// into dense matrix, then many row ops
		//...

	}
	else {
		printMatrix (A);
	}

	return 0 ;
}

// Local Variables:
// mode: C++
// tab-width: 4
// indent-tabs-mode: nil
// c-basic-offset: 4
// End:
// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s