Valence of sparse matrix over Z or Zp.
#include <givaro/modular.h>
#include <givaro/givintnumtheo.h>
#include <givaro/gf2.h>
#include <linbox/field/field-traits.h>
#include <linbox/blackbox/transpose.h>
#include <linbox/blackbox/compose.h>
#include <linbox/solutions/rank.h>
#include <linbox/solutions/valence.h>
#include <linbox/algorithms/smith-form-sparseelim-local.h>
#include <linbox/util/matrix-stream.h>
#include <linbox/util/error.h>
template<class Field>
unsigned long& TempLRank(unsigned long& r, char * filename, const Field& F)
{
std::ifstream input(filename);
LinBox::SparseMatrix<Field,LinBox::SparseMatrixFormat::SparseSeq> FA(msf);
input.close();
LinBox::Timer tim; tim.start();
LinBox::rankin(r, FA);
tim.stop();
F.write(std::cerr << "Rank over ") << " is " << r << ' ' << tim << " on T" << omp_get_thread_num() << std::endl;
return r;
}
unsigned long& TempLRank(unsigned long& r, char * filename, const LinBox::GF2& F2)
{
std::ifstream input(filename);
input.close();
LinBox::Timer tim; tim.start();
tim.stop();
F2.write(std::cerr << "Rank over ") << " is " << r << ' ' << tim << " on T" << omp_get_thread_num() << std::endl;
return r;
}
unsigned long& LRank(unsigned long& r, char * filename,Givaro::Integer p)
{
if (p == 2) {
LinBox::GF2 F2;
return TempLRank(r, filename, F2);
}
else if (p <= maxmod16) {
typedef Givaro::Modular<int16_t> Field;
Field F(p);
return TempLRank(r, filename, F);
}
else if (p <= maxmod32) {
typedef Givaro::Modular<int32_t> Field;
Field F(p);
return TempLRank(r, filename, F);
}
else if (p <= maxmod53) {
typedef Givaro::Modular<double> Field;
Field F(p);
return TempLRank(r, filename, F);
}
else if (p <= maxmod64) {
typedef Givaro::Modular<int64_t> Field;
Field F(p);
return TempLRank(r, filename, F);
}
else {
typedef Givaro::Modular<Givaro::Integer> Field;
Field F(p);
return TempLRank(r, filename, F);
}
return r;
}
std::vector<size_t>& PRank(std::vector<size_t>& ranks, size_t& effective_exponent, char * filename,Givaro::Integer p, size_t e, size_t intr)
{
effective_exponent = e;
Givaro::Integer maxmod;
if (p <= maxmod) {
typedef Givaro::Modular<int64_t> Ring;
int64_t lp(p);
Givaro::Integer q = pow(p,uint64_t(e)); int64_t lq(q);
if (q >Givaro::Integer(lq)) {
std::cerr << "Power rank might need extra large composite (" << p << '^' << e << ")." << std::endl;
q = p;
for(effective_exponent=1; q <= Ring::maxCardinality(); ++effective_exponent) {
q *= p;
}
q/=p; --effective_exponent;
lq = (int64_t)q;
std::cerr << "First trying: " << lq << " (=" << p << '^' << effective_exponent << ", without further warning this will be sufficient)." << std::endl;
}
Ring F(lq);
std::ifstream input(filename);
LinBox::SparseMatrix<Ring,LinBox::SparseMatrixFormat::SparseSeq > A (ms);
input.close();
LinBox::Permutation<Ring> Q(F,A.coldim());
LinBox::Timer tim; tim.clear(); tim.start();
PGD.
prime_power_rankin( lq, lp, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
tim.stop();
F.write(std::cerr << "Ranks over ") << " are " ;
for(std::vector<size_t>::const_iterator rit=ranks.begin(); rit != ranks.end(); ++rit)
std::cerr << *rit << ' ';
std::cerr << ' ' << tim << " on T" << omp_get_thread_num() << std::endl;
}
else {
std::cerr << "*** WARNING *** Sorry power rank mod large composite not yet implemented" << std::endl;
std::cerr << "*** WARNING *** Assuming integer rank, extra factors in the Smith form could be missing" << std::endl;
ranks.resize(0); ranks.push_back(intr);
}
return ranks;
}
#include <linbox/field/gf2.h>
#include <linbox/algorithms/smith-form-sparseelim-poweroftwo.h>
std::vector<size_t>& PRankPowerOfTwo(std::vector<size_t>& ranks, size_t& effective_exponent, char * filename, size_t e, size_t intr)
{
effective_exponent = e;
if (e > 63) {
std::cerr << "Power rank power of two might need extra large composite (2^" << e << ")." << std::endl;
std::cerr << "First trying: 63, without further warning this will be sufficient)." << std::endl;
effective_exponent = 63;
}
std::ifstream input(filename);
typedef Givaro::ZRing<int64_t> Ring;
Ring F;
LinBox::SparseMatrix<Ring,LinBox::SparseMatrixFormat::SparseSeq > A (ms);
input.close();
LinBox::GF2 F2;
LinBox::Permutation<LinBox::GF2> Q(F2,A.coldim());
LinBox::Timer tim; tim.clear(); tim.start();
PGD.
prime_power_rankin( effective_exponent, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
tim.stop();
F.write(std::cerr << "Ranks over ") << " modulo 2^" << effective_exponent << " are " ;
for(std::vector<size_t>::const_iterator rit=ranks.begin(); rit != ranks.end(); ++rit)
std::cerr << *rit << ' ';
std::cerr << ' ' << tim << " on T" << omp_get_thread_num() << std::endl;
return ranks;
}
std::vector<size_t>& PRankInteger(std::vector<size_t>& ranks, char * filename,Givaro::Integer p, size_t e, size_t intr)
{
typedef Givaro::Modular<Givaro::Integer> Ring;
Givaro::Integer q = pow(p,uint64_t(e));
Ring F(q);
std::ifstream input(filename);
LinBox::SparseMatrix<Ring,LinBox::SparseMatrixFormat::SparseSeq > A (ms);
input.close();
LinBox::Permutation<Ring> Q(F,A.coldim());
LinBox::Timer tim; tim.clear(); tim.start();
PGD.prime_power_rankin( q, p, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
tim.stop();
F.write(std::cerr << "Ranks over ") << " are " ;
for(std::vector<size_t>::const_iterator rit=ranks.begin(); rit != ranks.end(); ++rit)
std::cerr << *rit << ' ';
std::cerr << ' ' << tim << std::endl;
return ranks;
}
std::vector<size_t>& PRankIntegerPowerOfTwo(std::vector<size_t>& ranks, char * filename, size_t e, size_t intr)
{
typedef Givaro::ZRing<Givaro::Integer> Ring;
Ring ZZ;
std::ifstream input(filename);
LinBox::SparseMatrix<Ring,LinBox::SparseMatrixFormat::SparseSeq > A (ms);
input.close();
LinBox::Permutation<Ring> Q(ZZ, A.coldim());
LinBox::Timer tim; tim.clear(); tim.start();
PGD.
prime_power_rankin( e, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
tim.stop();
ZZ.write(std::cerr << "Ranks over ") << " modulo 2^" << e << " are " ;
for(std::vector<size_t>::const_iterator rit=ranks.begin(); rit != ranks.end(); ++rit)
std::cerr << *rit << ' ';
std::cerr << ' ' << tim << std::endl;
return ranks;
}