generate an example matrix with specified frobenius form.
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/ring/ntl.h>
#include <NTL/ZZX.h>
using std::string;
using std::list;
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)
{
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;
if ( *code != 'a')
{
p += ZZX(0, a);
}
else
{
int n = atoi(code+1);
R.init(a, n-1);
p += ZZX(n, R.one);
p += ZZX(1, R.one);
p += ZZX(0, a);
}
ZZX q(0, R.one);
for(int i = 0; i < e; ++i) q *= p;
LinBox::DenseVector<Ring> v(R,deg(q)+1);
for (int i = 0; i < v.size(); ++i) v[i] = coeff(q, i);
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)
{
int N,n = M.rowdim();
N = 2*n;
for (int k = 0; k < N; ++k) {
int i = rand()%M.rowdim();
int j = rand()%M.coldim();
if (i == j) continue;
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);
}
}
}
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)
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 Companion<Ring> BB;
int acp; char* avp[ac];
string opts;
stripOptions(acp, avp, opts, ac, av);
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);
if (opts.size() >= 1)
{ if (opts[0] == 'r')
{
DenseMatrix<Ring> B(Z,A.rowdim(), A.coldim());
LinBox::MatrixHom::map (B, A);
scramble(B);
printMatrix(B);
}
if (opts[0] == 'R') { throw LinBox::NotImplementedYet(); }
}
else {
printMatrix (A);
}
return 0 ;
}
A BlasVector<_Field > represents a vector as an array of _Field::Elements.
Dense matrix representation.
Definition: blas-matrix.h:62
If C = DirectSum(A, B) and y = xA and z = wB, then (y,z) = (x,w)C.
Definition: direct-sum.h:58
the integer ring.
Definition: ring/ntl/ntl-zz.h:114
Ring (in fact, a unique factorization domain) of polynomial with coefficients in class NTL_zz_p (inte...
Definition: ntl-lzz_px.h:78
bool isZero(const Element &x) const
Test if an element equals zero.
Definition: ntl-lzz_px.h:222
Matrix Homomorphism A map function converts a matrix on a field/ring to its natural image in another ...
Companion matrix of a monic polynomial.
Definition: companion.h:41