/usr/include/linbox/blackbox/ntl-sylvester.inl is in liblinbox-dev 1.1.6~rc0-4.1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 | /* -*- mode: C++; tab-width: 6; indent-tabs-mode: t; c-basic-offset: 4 -*- */
/*-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
* ntl-sylvester.h
* Copyright (C) 2003 Austin Lobo, B. David Saunders
* Member functions for the sylvester matrix in one variable
* for polynomials in one variable.
* Linbox version 2003
*-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+*/
#include <iostream>
namespace LinBox
{
/*----------------------------------------------------------------------
* Destructor
*---------------------------------------------------------------------*/
template <class Field>
inline Sylvester<Field>::~Sylvester()
{
#ifdef DBGMSGS
std::cout << "Sylvester::~Sylvester():\tDestroyed a "
<< rowDim << "x"<< colDim << " Sylvester matrix "<< std::endl;
#endif
}// End, ~Sylvester()
/*----------------------------------------------------------------------
* Default constructor.
*---------------------------------------------------------------------*/
template <class Field>
Sylvester<Field>::Sylvester()
{
sysDim = // Default dimension is 0
rowDim = // Default row dim is 0
colDim = 0; // Default col dim is 0
#ifdef DBGMSGS
std::cout << "Sylvester::Sylvester():\tCreated a " << rowDim << "x"<< colDim<<
" Sylvester matrix "<< std::endl;
#endif
}
/*----------------------------------------------------------------------
* Constructor. Builds two rectangular Toeplitz blocks, the top from
* px with deg(qx) rows, and the bottom out qx, with deg(px)
* rows.
*---------------------------------------------------------------------*/
template <class Field>
Sylvester<Field>::Sylvester(
const Field F,
const std::vector<typename Field::Element> &vp,
const std::vector<typename Field::Element> &vq
) : pdata(vp), qdata(vq)
{
// Set the row, col, and system dimensions
rowDim =
colDim =
sysDim = vp.size() + vq.size()-2; // square system with dimension deg(px)+deg(qx)
// Copy the vector to a polynomial representation, for px
pxdata.SetMaxLength( vp.size() );
for (size_t ip=0; ip< vp.size(); ip++ )
SetCoeff( pxdata, ip, vp[ip]);
// Copy the vector to a polynomial representation, for qx
qxdata.SetMaxLength( vq.size() );
for ( size_t iq=0; iq< vq.size(); iq++ )
SetCoeff( qxdata, iq, vq[iq] );
#ifdef DBGMSGS
std::cout << "Sylvester::Sylvester(F,v,w):\tCreated a " << rowDim <<
"x"<< colDim<< " Sylvester matrix "<< std::endl;
#endif
/* Inserted for a timing calculation with the gcd
* double start_time = GetTime();
* ZZ_pX gx;
* GCD( gx, pxdata, qxdata);
*
* double end_time = GetTime();
* cout << gx << endl;
* std::cout << "Time for gcd = " << end_time-start_time << " sec\n";
*/
}// Sylvester()
/*----------------------------------------------------------------------
* Prints to an ouput stream or to stdout by default -- useful for
* debugging
*---------------------------------------------------------------------*/
template <class Field>
void Sylvester<Field>::print(std::ostream& os) const
{
if ( sysDim < 20 )
{
os << rowDim << " " << colDim << std::endl;
for (size_t irow=0; irow < qxdeg(); irow++ )
{
os << "[ ";
for (size_t icoeff=0; icoeff < irow; icoeff++ )
os << "0 ";
for (size_t icoeff = pxdeg(); icoeff > 0; icoeff-- )
os << pdata[ icoeff ] << " ";
os << pdata[ 0 ] << " ";
for (size_t icoeff = sysDim-1-pxdeg()-irow; icoeff > 0; icoeff-- )
os << "0 ";
os << " ]\n";
}
for (size_t irow=0; irow < pxdeg(); irow++ )
{
os << "[ ";
for (size_t icoeff=0; icoeff < irow; icoeff++ )
os << "0 ";
for (size_t icoeff = qxdeg(); icoeff > 0; icoeff-- )
os << qdata[ icoeff ] << " ";
os << qdata[ 0 ] << " ";
for (size_t icoeff = sysDim-1-qxdeg()-irow; icoeff > 0; icoeff-- )
os << "0 ";
os << " ]\n";
}
}
else
{
os << rowDim << " " << colDim << std::endl;
os << pxdata;
os << "\n";
os << qxdata;
os << "\n";
}
}// Print()
/*----------------------------------------------------------------------
* Print to a file. By default print to stdout
*---------------------------------------------------------------------*/
template <class Field>
void Sylvester<Field>::print(char *outFileName) const
{
if ( outFileName == NULL )
print();
else {
std::ofstream o_fp(outFileName, std::ios::out);
o_fp << rowDim << " " << colDim << std::endl;
o_fp << "[";
for (size_t icoeff = 0; icoeff <= pxdeg() ; icoeff++ )
o_fp << pdata[ icoeff ] << " ";
o_fp << "]\n";
o_fp << "[";
for (size_t icoeff = 0; icoeff <= qxdeg(); icoeff++ )
o_fp << qdata[ icoeff ] << " ";
o_fp << "]\n";
o_fp.close();
}
}// Print
/*----------------------------------------------------------------------
* Print to a file in the format specified by the French Linbox group
* for LUDivine. Format is <nrow> <ncol> M
* <row> <col> <entry>
* .....
* 0 0 0
*--------------------------------------------------------------------- */
template <class Field>
void Sylvester<Field>::printcp(char *outFileName) const
{
if ( outFileName == NULL )
print();
else {
std::ofstream os(outFileName, std::ios::out);
os << rowDim << " " << colDim << " M" << std::endl;
size_t irow;
size_t icoeff;
for (irow=0; irow < qxdeg(); irow++ )
{
for (icoeff = pxdeg(); icoeff > 0; icoeff-- )
os << (irow+1) << " " << (irow + pxdeg()- icoeff +1) << " " <<
pdata[ icoeff ] << "\n";
os << (irow+1) << " " << (irow + pxdeg()+1) << " " <<
pdata[ icoeff ] << "\n";
}
for (irow=0; irow < pxdeg(); irow++ )
{
for (icoeff = qxdeg(); icoeff > 0; icoeff-- )
os << (irow+qxdeg()+1) << " " << (irow + pxdeg()- icoeff +1) << " " <<
qdata[ icoeff ] << "\n";
os << (irow+qxdeg()+1) << " " << (irow + pxdeg()+1) << " " <<
qdata[ icoeff ] << "\n";
}
os << "0 0 0\n";
os.close();
}
}
/* *----------------------------------------------------------------------
* Applytranspose. Does 2 FFT's each of degree (deg(qx)+deg(px))
* Speed can be improved to 1 FFT in reverse direction
* by saving pre-computed FFT values
*---------------------------------------------------------------------*/
template <class Field>
template <class OutVector, class InVector>
OutVector& Sylvester<Field>::apply( OutVector &v_out,
const InVector& v_in) const
{
/* uncomment the following lines here and in apply() and swap method names
* to swap the apply() for the faster applyTranspose()
*
* static int flipper = 1;
* if (flipper == 1)
* { cout << "\t sylvester applyTranspose is flipped apply\n";
* flipper++;
* }
*/
/*--------------- Check the size of the output vector ----------*/
if ( v_out.size() != sysdim() )
std::cout << "\tSylvester::apply()\t output vector not correct size, at "
<< v_out.size() << ". System rowdim is" << sysdim() << std::endl;
NTL::ZZ_pX txOut, txIn;
/*--------------- Convert input vector to a polynomial ---------*/
txIn.SetMaxLength( v_in.size() -1 );
for ( size_t i=0; i < v_in.size(); i++ )
SetCoeff( txIn, i, v_in[i] );
/*-------------- Poly multiply the upper Sylvester poly by input ------*/
mul( txOut, txIn, pxdata);
int Nq = qxdeg();
int m = pxdeg();
/*-------------- vout[0..deg(q)-1] <--- txout[deg(qx)...2deg(qx)-1] --- */
for ( int i=0; i < Nq; ++i )
GetCoeff(v_out[i], txOut, m+i);
/*-------------- Poly multiply the lower Sylvester poly by input -----*/
mul( txOut, txIn, qxdata);
int Np = pxdeg();
int n = qxdeg();
/*-------------- vout[deg(qx)..deg(qx)+deg(px)-1] <---
* txout[deg(qx)...deg(qx)+deg(px)-1] --- */
for ( int i=0; i < Np; ++i )
GetCoeff(v_out[Nq+i], txOut, n+i );
return v_out;
}// Apply -- Tested and works 3/17/03
/*----------------------------------------------------------------------
* Apply Transpose: Does 2 polymults of degree deg(qx), and deg(px)
* respectively. It takes approx 55% of the time of
* apply, in apparent violation of Tellegen's theorem
* Actually, both apply and applyT can be done with
* just one FFT per call, by saving previous FFTs
* we use deg(qx) = m; and deg(px) = n; N=m+n below
*---------------------------------------------------------------------*/
template <class Field>
template <class OutVector, class InVector>
OutVector& Sylvester<Field>::applyTranspose( OutVector &v_out,
const InVector& v_in) const
{
/* uncomment the following lines here and in apply() and swap method names
* to swap the apply() for the faster applyTranspose()
*
* static int flipper=1;
* if (flipper == 1)
* { cout << "\tSylvester apply is flipped applyTranspose\n"; flipper++;}
*/
/*--------- Check for size-compatibility of output vectors ------------*/
if ( v_out.size() != sysdim() )
std::cout << "\tSylvester::apply()\t output vector not correct size, at "
<< v_out.size() << ". System rowdim is" << sysdim() << std::endl;
NTL::ZZ_pX txOut, txIn;
NTL::ZZ_p tval;
/*--------- We need to reverse the input vector -----------*/
/*--------- txIn[0...m-1] <--- vin[m-1...0] ----------*/
txIn.SetMaxLength( qxdeg()-1 );
for ( size_t i=0; i < qxdeg(); i++ )
SetCoeff( txIn, qxdeg()-1-i, v_in[i] );
mul( txOut, txIn, pxdata); // Do the poly multiply
/*--------- We need to reverse the output vector -----------*/
/*------- v_out[0..N-1] <--- txOut[N-1...0] --------*/
for (size_t i=0; i < v_out.size(); i++)
GetCoeff(v_out[i], txOut, sysdim()-1-i); // Extract the coeffs
/*--------- We need to reverse the input vector -----------*/
/*--------- txIn[m ... m+n-1] <--- vin[m+n-1 ... m] -----------*/
txIn.SetMaxLength( pxdeg()-1 );
for ( size_t i=0; i < pxdeg(); i++ )
SetCoeff( txIn,pxdeg()-1-i, v_in[qxdeg()+i] );
mul( txOut, txIn, qxdata); // do the poly multiply
/*--------- We need to reverse the input vector -----------*/
/*------- v_out[0..N-1] <--- txOut[N-1...0] --------*/
for (size_t i=0; i < v_out.size(); i++) {
GetCoeff(tval, txOut, sysdim()-1-i);
add( v_out[i], v_out[i], tval ); // add to the accumulated output
}
return v_out;
}
}//End, namespace LinBox
|