1 /* linbox/ring/modular.h 2 * Copyright (C) 1999-2001 William J Turner, 3 * 2001 Bradford Hovinen 4 * Copyright (C) 2011 LinBox 5 * 6 * Written by William J Turner <wjturner@math.ncsu.edu>, 7 * Bradford Hovinen <hovinen@cis.udel.edu> 8 * 9 * ------------------------------------ 10 * 2002-04-10 Bradford Hovinen <hovinen@cis.udel.edu> 11 * 12 * LargeModular is now replace by a class Givaro::Modular parameterized on the element 13 * type. So, the old LargeModular is equivalent to Givaro::Modular<integer>. All other 14 * interface details are exactly the same. 15 * 16 * Renamed from large-modular.h to modular.h 17 * ------------------------------------ 18 * 19 * 20 * ========LICENCE======== 21 * This file is part of the library LinBox. 22 * 23 * LinBox is free software: you can redistribute it and/or modify 24 * it under the terms of the GNU Lesser General Public 25 * License as published by the Free Software Foundation; either 26 * version 2.1 of the License, or (at your option) any later version. 27 * 28 * This library is distributed in the hope that it will be useful, 29 * but WITHOUT ANY WARRANTY; without even the implied warranty of 30 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 31 * Lesser General Public License for more details. 32 * 33 * You should have received a copy of the GNU Lesser General Public 34 * License along with this library; if not, write to the Free Software 35 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 36 * ========LICENCE======== 37 *. 38 */ 39 40 #ifndef __LINBOX_field_modular_unsigned_H 41 #define __LINBOX_field_modular_unsigned_H 42 43 //Dan Roche 7-2-04 44 #ifndef __LINBOX_MIN 45 #define __LINBOX_MIN(a,b) ( (a) < (b) ? (a) : (b) ) 46 #endif 47 48 namespace LinBox { /* uint8_t */ 49 50 /*! Specialization of FieldAXPY for uint8_t modular field */ 51 52 template<class Compute_t> 53 class FieldAXPY<Givaro::Modular<uint8_t,Compute_t > > { 54 public: 55 56 typedef uint8_t Element; 57 typedef uint64_t Abnormal; 58 typedef Givaro::Modular<uint8_t, Compute_t> Field; 59 FieldAXPY(const Field & F)60 FieldAXPY (const Field &F) : 61 _k (((uint64_t) -1LL) / ((F.characteristic() - 1) * (F.characteristic() - 1))), 62 _field (&F), 63 _y (0), 64 i (_k) 65 { 66 } 67 FieldAXPY(const FieldAXPY & faxpy)68 FieldAXPY (const FieldAXPY &faxpy) : 69 _k (faxpy._k), 70 _field (faxpy._field), 71 _y (0), 72 i (_k) 73 {} 74 75 FieldAXPY<Field> &operator = (const FieldAXPY &faxpy) 76 { 77 _field = faxpy._field; 78 _y = faxpy._y; 79 _k = faxpy._k; 80 return *this; 81 } 82 mulacc(const Element & a,const Element & x)83 inline uint64_t& mulacc (const Element &a, const Element &x) 84 { 85 uint32_t t = (uint32_t) a * (uint32_t) x; 86 87 if (!i--) { 88 i = int(_k); 89 return _y = _y % (uint32_t) field().characteristic() + t; 90 } 91 else 92 return _y += t; 93 } 94 accumulate(const Element & t)95 inline uint64_t& accumulate (const Element &t) 96 { 97 98 if (!i--) { 99 i = int(_k); 100 return _y = _y % (uint32_t) field().characteristic() + t; 101 } 102 else 103 return _y += t; 104 } 105 get(Element & y)106 inline Element &get (Element &y) const 107 { 108 const_cast<FieldAXPY<Field>*>(this)->_y %= (uint32_t) field().characteristic(); 109 if ((int32_t) _y < 0) const_cast<FieldAXPY<Field>*>(this)->_y += field().characteristic(); 110 y = (uint8_t) _y; 111 const_cast<FieldAXPY<Field>*>(this)->i = int(_k); 112 return y; 113 } 114 assign(const Element y)115 inline FieldAXPY &assign (const Element y) 116 { 117 _y = y; 118 i = int(_k); 119 return *this; 120 } 121 reset()122 inline void reset() 123 { 124 _y = 0; 125 } 126 field()127 inline const Field & field() const { return *_field; } 128 129 public: 130 uint64_t _k; 131 132 private: 133 const Field *_field; 134 uint64_t _y; 135 int64_t i; 136 }; 137 138 //! Specialization of DotProductDomain for unsigned short modular field 139 140 template <class Compute_t> 141 class DotProductDomain<Givaro::Modular<uint8_t, Compute_t> > : public VectorDomainBase<Givaro::Modular<uint8_t, Compute_t> > { 142 public: 143 144 typedef uint8_t Element; 145 typedef Givaro::Modular<uint8_t, Compute_t> Field; 146 DotProductDomain()147 DotProductDomain(){} DotProductDomain(const Field & F)148 DotProductDomain (const Field &F) : 149 VectorDomainBase<Field> (F) 150 {} 151 using VectorDomainBase<Field>::field; 152 using VectorDomainBase<Field>::faxpy; 153 154 protected: 155 template <class Vector1, class Vector2> dotSpecializedDD(Element & res,const Vector1 & v1,const Vector2 & v2)156 inline Element &dotSpecializedDD (Element &res, const Vector1 &v1, const Vector2 &v2) const 157 { 158 typename Vector1::const_iterator i = v1.begin (); 159 typename Vector2::const_iterator j = v2.begin (); 160 161 typename Vector1::const_iterator iterend = v1.begin () + (ptrdiff_t)(v1.size() % faxpy()._k); 162 163 uint64_t y = 0; 164 165 for (; i != iterend; ++i, ++j) 166 y += (uint64_t) *i * (uint64_t) *j; 167 168 y %= (uint64_t) field().characteristic(); 169 170 for (; iterend != v1.end (); j += (ptrdiff_t)faxpy()._k) { 171 typename Vector1::const_iterator iter_i = iterend; 172 typename Vector2::const_iterator iter_j; 173 174 iterend += (ptrdiff_t)faxpy()._k; 175 176 for (iter_j = j; iter_i != iterend; ++iter_i, ++iter_j) 177 y += (uint64_t) *iter_i * (uint64_t) *j; 178 179 y %= (uint64_t) field().characteristic(); 180 } 181 182 return res = (uint8_t) y; 183 } 184 185 template <class Vector1, class Vector2> dotSpecializedDSP(Element & res,const Vector1 & v1,const Vector2 & v2)186 inline Element &dotSpecializedDSP (Element &res, const Vector1 &v1, const Vector2 &v2) const 187 { 188 typename Vector1::first_type::const_iterator i_idx = v1.first.begin (); 189 typename Vector1::second_type::const_iterator i_elt = v1.second.begin (); 190 191 uint64_t y = 0; 192 193 if (v1.first.size () < faxpy()._k) { 194 for (; i_idx != v1.first.end (); ++i_idx, ++i_elt) 195 y += (uint64_t) *i_elt * (uint64_t) v2[*i_idx]; 196 197 return res = uint8_t (y % (uint64_t) field().characteristic()); 198 } 199 else { 200 typename Vector1::first_type::const_iterator iterend = v1.first.begin () +(ptrdiff_t)( v1.first.size() % faxpy()._k); 201 202 for (; i_idx != iterend; ++i_idx, ++i_elt) 203 y += (uint64_t) *i_elt * (uint64_t) v2[*i_idx]; 204 205 y %= (uint64_t) field().characteristic(); 206 207 while (iterend != v1.first.end ()) { 208 typename Vector1::first_type::const_iterator iter_i_idx = iterend; 209 typename Vector1::second_type::const_iterator iter_i_elt = i_elt; 210 211 iterend += (ptrdiff_t)faxpy()._k; 212 i_elt += (ptrdiff_t)faxpy()._k; 213 214 for (; iter_i_idx != iterend; ++iter_i_idx, ++iter_i_elt) 215 y += (uint64_t) *iter_i_elt * (uint64_t) v2[*iter_i_idx]; 216 217 y %= (uint64_t) field().characteristic(); 218 } 219 220 return res = (uint8_t) y; 221 } 222 } 223 }; 224 225 //! Specialization of MVProductDomain for uint8_t modular field 226 227 template<class Compute_t> 228 class MVProductDomain<Givaro::Modular<uint8_t,Compute_t> > { 229 public: 230 231 typedef uint8_t Element; 232 typedef Givaro::Modular<uint8_t,Compute_t> Field; 233 234 protected: 235 template <class Vector1, class Matrix, class Vector2> mulColDense(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v)236 inline Vector1 &mulColDense 237 (const VectorDomain<Field> &VD, Vector1 &w, const Matrix &A, const Vector2 &v) const 238 { 239 return mulColDenseSpecialized (VD, w, A, v, typename VectorTraits<typename Matrix::Column>::VectorCategory ()); 240 } 241 242 private: 243 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::DenseVectorTag)244 Vector1 &mulColDenseSpecialized 245 (const VectorDomain<Field> &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 246 VectorCategories::DenseVectorTag) const 247 { 248 linbox_check (A.coldim () == v.size ()); 249 linbox_check (A.rowdim () == w.size ()); 250 251 typename Matrix::ConstColIterator i = A.colBegin (); 252 typename Vector2::const_iterator j, j_end; 253 typename Matrix::Column::const_iterator k; 254 std::vector<uint32_t>::iterator l, l_end; 255 256 if (_tmp.size () < w.size ()) 257 _tmp.resize (w.size ()); 258 259 std::fill (_tmp.begin (), _tmp.begin () + (ptrdiff_t)w.size (), 0); 260 261 l_end = _tmp.begin () +(ptrdiff_t) w.size (); 262 263 do { 264 j = v.begin (); 265 j_end = j + __LINBOX_MIN (A->coldim (), VD.faxpy()._k); 266 267 for (; j != j_end; ++j, ++i) 268 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) 269 *l += *k * *j; 270 271 j_end += __LINBOX_MIN (A->coldim () - (j_end - v.begin ()), VD.faxpy()._k); 272 273 for (l =_tmp.begin (); l != l_end; ++l) 274 *l %= VD.field ().characteristic(); 275 276 } while (j_end != v.end ()); 277 278 typename Vector1::iterator w_j; 279 280 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 281 *w_j = *l; 282 283 return w; 284 } 285 286 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseSequenceVectorTag)287 Vector1 &mulColDenseSpecialized 288 (const VectorDomain<Field> &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 289 VectorCategories::SparseSequenceVectorTag) const 290 { 291 linbox_check (A.coldim () == v.size ()); 292 linbox_check (A.rowdim () == w.size ()); 293 294 typename Matrix::ConstColIterator i = A.colBegin (); 295 typename Vector2::const_iterator j, j_end; 296 typename Matrix::Column::const_iterator k; 297 std::vector<uint32_t>::iterator l, l_end; 298 299 if (_tmp.size () < w.size ()) 300 _tmp.resize (w.size ()); 301 302 std::fill (_tmp.begin (), _tmp.begin () + (ptrdiff_t)w.size (), 0); 303 304 l_end = _tmp.begin () + (ptrdiff_t)w.size (); 305 306 307 do { 308 j = v.begin (); 309 j_end = j + __LINBOX_MIN (A->coldim (), VD.faxpy()._k); 310 311 for (; j != j_end; ++j, ++i) 312 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) 313 _tmp[k->first] += k->second * *j; 314 315 j_end += __LINBOX_MIN (A->coldim () - (j_end - v.begin ()), VD.faxpy()._k); 316 317 for (l =_tmp.begin (); l != l_end; ++l) 318 *l %= VD.field ().characteristic(); 319 320 } while (j_end != v.end ()); 321 322 typename Vector1::iterator w_j; 323 324 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 325 *w_j = *l; 326 327 return w; 328 } 329 330 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseAssociativeVectorTag)331 Vector1 &mulColDenseSpecialized 332 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 333 VectorCategories::SparseAssociativeVectorTag) const 334 { 335 linbox_check (A.coldim () == v.size ()); 336 linbox_check (A.rowdim () == w.size ()); 337 338 typename Matrix::ConstColIterator i = A.colBegin (); 339 typename Vector2::const_iterator j, j_end; 340 typename Matrix::Column::const_iterator k; 341 std::vector<uint32_t>::iterator l, l_end; 342 343 if (_tmp.size () < w.size ()) 344 _tmp.resize (w.size ()); 345 346 std::fill (_tmp.begin (), _tmp.begin () + (ptrdiff_t)w.size (), 0); 347 348 l_end = _tmp.begin () +(ptrdiff_t) w.size (); 349 350 do { 351 j = v.begin (); 352 j_end = j + __LINBOX_MIN (A->coldim (), VD.faxpy()._k); 353 354 for (; j != j_end; ++j, ++i) 355 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) 356 _tmp[k->first] += k->second * *j; 357 358 j_end += __LINBOX_MIN (A->coldim () - (j_end - v.begin ()), VD.faxpy()._k); 359 360 for (l =_tmp.begin (); l != l_end; ++l) 361 *l %= VD.field ().characteristic(); 362 363 } while (j_end != v.end ()); 364 365 typename Vector1::iterator w_j; 366 367 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 368 *w_j = *l; 369 370 return w; 371 } 372 373 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseParallelVectorTag)374 Vector1 &mulColDenseSpecialized 375 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 376 VectorCategories::SparseParallelVectorTag) const 377 { 378 linbox_check (A.coldim () == v.size ()); 379 linbox_check (A.rowdim () == w.size ()); 380 381 typename Matrix::ConstColIterator i = A.colBegin (); 382 typename Vector2::const_iterator j, j_end; 383 typename Matrix::Column::first_type::const_iterator k_idx; 384 typename Matrix::Column::second_type::const_iterator k_elt; 385 std::vector<uint32_t>::iterator l, l_end; 386 387 if (_tmp.size () < w.size ()) 388 _tmp.resize (w.size ()); 389 390 std::fill (_tmp.begin (), _tmp.begin () + (ptrdiff_t)w.size (), 0); 391 392 l_end = _tmp.begin () + (ptrdiff_t)w.size (); 393 394 do { 395 j = v.begin (); 396 j_end = j + (ptrdiff_t)__LINBOX_MIN (uint64_t (A.coldim ()), VD.faxpy()._k); 397 398 for (; j != j_end; ++j, ++i) 399 for (k_idx = i->first.begin (), k_elt = i->second.begin (), l = _tmp.begin (); 400 k_idx != i->first.end (); 401 ++k_idx, ++k_elt, ++l) 402 _tmp[*k_idx] += *k_elt * *j; 403 404 j_end += (ptrdiff_t) __LINBOX_MIN (uint64_t (A.coldim () - (size_t)(j_end - v.begin ())), VD.faxpy()._k); 405 406 for (l =_tmp.begin (); l != l_end; ++l) 407 *l %= VD.field ().characteristic(); 408 409 } while (j_end != v.end ()); 410 411 typename Vector1::iterator w_j; 412 typedef typename Vector1::value_type val_t ; 413 414 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 415 *w_j = (val_t) *l; 416 417 return w; 418 } 419 420 421 mutable std::vector<uint32_t> _tmp; 422 }; 423 424 } 425 426 namespace LinBox { /* uint16_t */ 427 428 /*! Specialization of FieldAXPY for uint16_t modular field */ 429 template<class Compute_t> 430 class FieldAXPY<Givaro::Modular<uint16_t,Compute_t> > { 431 public: 432 433 typedef uint16_t Element; 434 typedef Givaro::Modular<uint16_t,Compute_t> Field; 435 FieldAXPY(const Field & F)436 FieldAXPY (const Field &F) : 437 _k (((uint64_t) -1LL) / ((F.characteristic() - 1) * (F.characteristic() - 1))), 438 _field (&F), 439 _y (0), 440 i (_k) 441 {} 442 FieldAXPY(const FieldAXPY & faxpy)443 FieldAXPY (const FieldAXPY &faxpy) : 444 _k (faxpy._k), _field (faxpy._field), _y (0), i (_k) 445 {} 446 447 FieldAXPY<Field > &operator = (const FieldAXPY &faxpy) 448 { 449 _field = faxpy._field; 450 _y = faxpy._y; 451 _k = faxpy._k; 452 return *this; 453 } 454 mulacc(const Element & a,const Element & x)455 inline uint64_t& mulacc (const Element &a, const Element &x) 456 { 457 uint64_t t = (uint64_t) ((long long) a * (long long) x); 458 459 if (!i--) { 460 i = (int)_k; 461 return _y = _y % (uint64_t) field().characteristic() + t; 462 } 463 else 464 return _y += t; 465 } 466 accumulate(const Element & t)467 inline uint64_t& accumulate (const Element &t) 468 { 469 if (!i--) { 470 i = (int)_k; 471 return _y = _y % (uint64_t) field().characteristic() + t; 472 } 473 else 474 return _y += t; 475 } 476 get(Element & y)477 inline Element &get (Element &y) const 478 { 479 const_cast<FieldAXPY<Field>*>(this)->_y %= (uint64_t) field().characteristic(); 480 if ((int64_t) _y < 0) const_cast<FieldAXPY<Field>*>(this)->_y += field().characteristic(); 481 y = (uint16_t) _y; 482 const_cast<FieldAXPY<Field>*>(this)->i = int(_k); 483 return y; 484 } 485 assign(const Element y)486 inline FieldAXPY &assign (const Element y) 487 { 488 _y = y; 489 i = (int)_k; 490 return *this; 491 } 492 reset()493 inline void reset() 494 { 495 _y = 0; 496 } 497 field()498 inline const Field & field() const {return *_field;} 499 500 public: 501 uint64_t _k; 502 503 private: 504 const Field *_field; 505 uint64_t _y; 506 int64_t i; 507 }; 508 509 //! Specialization of DotProductDomain for unsigned short modular field 510 511 template<class Compute_t> 512 class DotProductDomain<Givaro::Modular<uint16_t,Compute_t> > : public VectorDomainBase<Givaro::Modular<uint16_t,Compute_t> > { 513 public: 514 515 typedef uint16_t Element; 516 typedef Givaro::Modular<uint16_t,Compute_t> Field; 517 DotProductDomain()518 DotProductDomain () {} DotProductDomain(const Field & F)519 DotProductDomain (const Field &F) : 520 VectorDomainBase<Field > (F) 521 {} 522 using VectorDomainBase<Field>::field; 523 using VectorDomainBase<Field>::faxpy; 524 525 protected: 526 template <class Vector1, class Vector2> dotSpecializedDD(Element & res,const Vector1 & v1,const Vector2 & v2)527 inline Element &dotSpecializedDD (Element &res, const Vector1 &v1, const Vector2 &v2) const 528 { 529 typename Vector1::const_iterator i = v1.begin (); 530 typename Vector2::const_iterator j = v2.begin (); 531 532 typename Vector1::const_iterator iterend = v1.begin () + (ptrdiff_t)(v1.size() % faxpy()._k); 533 534 uint64_t y = 0; 535 536 for (; i != iterend; ++i, ++j) 537 y += (uint64_t) *i * (uint64_t) *j; 538 539 y %= (uint64_t) field().characteristic(); 540 541 for (; iterend != v1.end (); j += faxpy()._k) { 542 typename Vector1::const_iterator iter_i = iterend; 543 typename Vector2::const_iterator iter_j; 544 545 iterend += faxpy()._k; 546 547 for (iter_j = j; iter_i != iterend; ++iter_i, ++iter_j) 548 y += (uint64_t) *iter_i * (uint64_t) *j; 549 550 y %= (uint64_t) field().characteristic(); 551 } 552 553 return res = (uint16_t) y; 554 } 555 556 557 template <class Vector1, class Vector2> dotSpecializedDSP(Element & res,const Vector1 & v1,const Vector2 & v2)558 inline Element &dotSpecializedDSP (Element &res, const Vector1 &v1, const Vector2 &v2) const 559 { 560 typename Vector1::first_type::const_iterator i_idx = v1.first.begin (); 561 typename Vector1::second_type::const_iterator i_elt = v1.second.begin (); 562 563 uint64_t y = 0; 564 565 if (v1.first.size () < faxpy()._k) { 566 for (; i_idx != v1.first.end (); ++i_idx, ++i_elt) 567 y += (uint64_t) *i_elt * (uint64_t) v2[*i_idx]; 568 569 return res = (uint16_t) (y % (uint64_t) field().characteristic()); 570 } 571 else { 572 typename Vector1::first_type::const_iterator iterend = v1.first.begin () +(ptrdiff_t)( v1.first.size() % faxpy()._k ); 573 574 for (; i_idx != iterend; ++i_idx, ++i_elt) 575 y += (uint64_t) *i_elt * (uint64_t) v2[*i_idx]; 576 577 y %= (uint64_t) field().characteristic(); 578 579 while (iterend != v1.first.end ()) { 580 typename Vector1::first_type::const_iterator iter_i_idx = iterend; 581 typename Vector1::second_type::const_iterator iter_i_elt = i_elt; 582 583 iterend += faxpy()._k; 584 i_elt += faxpy()._k; 585 586 for (; iter_i_idx != iterend; ++iter_i_idx, ++iter_i_elt) 587 y += (uint64_t) *iter_i_elt * (uint64_t) v2[*iter_i_idx]; 588 589 y %= (uint64_t) field().characteristic(); 590 } 591 592 return res = (Element) y; 593 } 594 } 595 596 }; 597 598 //! Specialization of MVProductDomain for uint16_t modular field 599 600 template<class Compute_t> 601 class MVProductDomain<Givaro::Modular<uint16_t,Compute_t> > { 602 public: 603 604 typedef uint16_t Element; 605 typedef Givaro::Modular<uint16_t,Compute_t> Field; 606 protected: 607 template <class Vector1, class Matrix, class Vector2> mulColDense(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v)608 inline Vector1 &mulColDense 609 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v) const 610 { 611 return mulColDenseSpecialized (VD, w, A, v, VectorTraits<typename Matrix::Column>::VectorCategory ()); 612 } 613 614 private: 615 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::DenseVectorTag)616 Vector1 &mulColDenseSpecialized 617 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 618 VectorCategories::DenseVectorTag) const 619 { 620 linbox_check (A.coldim () == v.size ()); 621 linbox_check (A.rowdim () == w.size ()); 622 623 typename Matrix::ConstColIterator i = A.colBegin (); 624 typename Vector2::const_iterator j = v.begin (), j_end; 625 typename Matrix::Column::const_iterator k; 626 // Dan Roche, 7-1-04 627 // std::vector<uint32_t>::iterator l, l_end; 628 std::vector<uint64_t>::iterator l, l_end; 629 630 if (_tmp.size () < w.size ()) 631 _tmp.resize (w.size ()); 632 633 std::fill (_tmp.begin (), _tmp.begin () +(ptrdiff_t) w.size (), 0); 634 635 l_end = _tmp.begin () +(ptrdiff_t) w.size (); 636 637 do { 638 j = v.begin (); 639 j_end = j + __LINBOX_MIN (A->coldim (), VD.faxpy()._k); 640 641 for (; j != j_end; ++j, ++i) 642 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) 643 *l += *k * *j; 644 645 j_end += __LINBOX_MIN (A->coldim () - (j_end - v.begin ()), VD.faxpy()._k); 646 647 for (l =_tmp.begin (); l != l_end; ++l) 648 *l %= VD.field ().characteristic(); 649 650 } while (j_end != v.end ()); 651 652 typename Vector1::iterator w_j; 653 654 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 655 *w_j = *l; 656 657 return w; 658 } 659 660 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseSequenceVectorTag)661 Vector1 &mulColDenseSpecialized 662 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 663 VectorCategories::SparseSequenceVectorTag) const 664 { 665 linbox_check (A.coldim () == v.size ()); 666 linbox_check (A.rowdim () == w.size ()); 667 668 typename Matrix::ConstColIterator i = A.colBegin (); 669 typename Vector2::const_iterator j, j_end; 670 typename Matrix::Column::const_iterator k; 671 // Dan Roche, 7-1-04 672 // std::vector<uint32_t>::iterator l, l_end; 673 std::vector<uint64_t>::iterator l, l_end; 674 675 if (_tmp.size () < w.size ()) 676 _tmp.resize (w.size ()); 677 678 std::fill (_tmp.begin (), _tmp.begin () +(ptrdiff_t) w.size (), 0); 679 680 l_end = _tmp.begin () +(ptrdiff_t) w.size (); 681 682 do { 683 j = v.begin (); 684 j_end = j + __LINBOX_MIN (A->coldim (), VD.faxpy()._k); 685 686 for (; j != j_end; ++j, ++i) 687 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) 688 _tmp[k->first] += k->second * *j; 689 690 j_end += __LINBOX_MIN (A->coldim () - (j_end - v.begin ()), VD.faxpy()._k); 691 692 for (l =_tmp.begin (); l != l_end; ++l) 693 *l %= VD.field ().characteristic(); 694 695 } while (j_end != v.end ()); 696 697 typename Vector1::iterator w_j; 698 699 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 700 *w_j = *l; 701 702 return w; 703 } 704 705 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseAssociativeVectorTag)706 Vector1 &mulColDenseSpecialized 707 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 708 VectorCategories::SparseAssociativeVectorTag) const 709 { 710 linbox_check (A.coldim () == v.size ()); 711 linbox_check (A.rowdim () == w.size ()); 712 713 typename Matrix::ConstColIterator i = A.colBegin (); 714 typename Vector2::const_iterator j, j_end; 715 typename Matrix::Column::const_iterator k; 716 // Dan Roche, 7-1-04 717 // std::vector<uint32_t>::iterator l, l_end; 718 std::vector<uint64_t>::iterator l, l_end; 719 720 if (_tmp.size () < w.size ()) 721 _tmp.resize (w.size ()); 722 723 std::fill (_tmp.begin (), _tmp.begin () +(ptrdiff_t) w.size (), 0); 724 725 l_end = _tmp.begin () +(ptrdiff_t) w.size (); 726 727 do { 728 j = v.begin (); 729 j_end = j + __LINBOX_MIN (A->coldim (), VD.faxpy()._k); 730 731 for (; j != j_end; ++j, ++i) 732 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) 733 _tmp[k->first] += k->second * *j; 734 735 j_end += __LINBOX_MIN (A->coldim () - (j_end - v.begin ()), VD.faxpy()._k); 736 737 for (l =_tmp.begin (); l != l_end; ++l) 738 *l %= VD.field ().characteristic(); 739 740 } while (j_end != v.end ()); 741 742 typename Vector1::iterator w_j; 743 744 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 745 *w_j = *l; 746 747 return w; 748 } 749 750 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseParallelVectorTag)751 Vector1 &mulColDenseSpecialized 752 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 753 VectorCategories::SparseParallelVectorTag) const 754 { 755 linbox_check (A.coldim () == v.size ()); 756 linbox_check (A.rowdim () == w.size ()); 757 758 typename Matrix::ConstColIterator i = A.colBegin (); 759 typename Vector2::const_iterator j, j_end; 760 typename Matrix::Column::first_type::const_iterator k_idx; 761 typename Matrix::Column::second_type::const_iterator k_elt; 762 // Dan Roche, 7-1-04 763 // std::vector<uint32_t>::iterator l, l_end; 764 std::vector<uint64_t>::iterator l, l_end; 765 766 if (_tmp.size () < w.size ()) 767 _tmp.resize (w.size ()); 768 769 std::fill (_tmp.begin (), _tmp.begin () +(ptrdiff_t) w.size (), 0); 770 771 l_end = _tmp.begin () +(ptrdiff_t) w.size (); 772 773 do { 774 j = v.begin (); 775 //Dan Roche, 7-2-04 776 //j_end = j + __LINBOX_MIN (A->coldim (), VD.faxpy()._k); 777 j_end = j + __LINBOX_MIN (A.coldim (), VD.faxpy()._k); 778 779 for (; j != j_end; ++j, ++i) 780 for (k_idx = i->first.begin (), k_elt = i->second.begin (), l = _tmp.begin (); 781 k_idx != i->first.end (); 782 ++k_idx, ++k_elt, ++l) 783 _tmp[*k_idx] += *k_elt * *j; 784 785 //j_end += __LINBOX_MIN (A->coldim () - (j_end - v.begin ()), VD.faxpy()._k); 786 j_end += __LINBOX_MIN (A.coldim () - (j_end - v.begin ()), VD.faxpy()._k); 787 788 for (l =_tmp.begin (); l != l_end; ++l) 789 *l %= VD.field ().characteristic(); 790 791 } while (j_end != v.end ()); 792 793 typename Vector1::iterator w_j; 794 795 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 796 *w_j = *l; 797 798 return w; 799 } 800 801 mutable std::vector<uint64_t> _tmp; 802 }; 803 804 } 805 806 #include <givaro/modular-integral.h> 807 808 namespace LinBox { /* uint32_t */ 809 810 template<class Field> 811 class DotProductDomain; 812 template<class Field> 813 class FieldAXPY; 814 template<class Field> 815 class MVProductDomain; 816 817 818 /*! Specialization of FieldAXPY for unsigned short modular field */ 819 820 template<class Compute_t> 821 class FieldAXPY<Givaro::Modular<uint32_t, Compute_t> > { 822 public: 823 824 typedef uint32_t Element; 825 typedef Givaro::Modular<uint32_t, Compute_t> Field; 826 FieldAXPY(const Field & F)827 FieldAXPY (const Field &F) : 828 _field (&F), _y(0) 829 { 830 _two_64 = (uint64_t(1) << 32) % uint64_t(F.characteristic()); 831 _two_64 = (_two_64 * _two_64) % uint64_t(F.characteristic()); 832 } 833 FieldAXPY(const FieldAXPY & faxpy)834 FieldAXPY (const FieldAXPY &faxpy) : 835 _two_64 (faxpy._two_64), _field (faxpy._field), _y (0) 836 {} 837 838 FieldAXPY<Field > &operator = (const FieldAXPY &faxpy) 839 { 840 _field = faxpy._field; 841 _y = faxpy._y; 842 _two_64 = faxpy._two_64; 843 return *this; 844 } 845 mulacc(const Element & a,const Element & x)846 inline uint64_t& mulacc (const Element &a, const Element &x) 847 { 848 uint64_t t = (uint64_t) a * (uint64_t) x; 849 _y += t; 850 851 if (_y < t) 852 return _y += _two_64; 853 else 854 return _y; 855 } 856 accumulate(const Element & t)857 inline uint64_t& accumulate (const Element &t) 858 { 859 _y += t; 860 861 if (_y < t) 862 return _y += _two_64; 863 else 864 return _y; 865 } 866 accumulate_special(const Element & t)867 inline uint64_t& accumulate_special (const Element &t) 868 { 869 return _y += t; 870 } 871 get(Element & y)872 inline Element &get (Element &y) const 873 { 874 const_cast<FieldAXPY<Field>*>(this)->_y %= (uint64_t) field().characteristic(); 875 //if ((int64_t) _y < 0) const_cast<FieldAXPY<Field>*>(this)->_y += field().characteristic(); 876 return y = (uint32_t) _y; 877 } 878 assign(const Element y)879 inline FieldAXPY &assign (const Element y) 880 { 881 _y = y; 882 return *this; 883 } 884 reset()885 inline void reset() { 886 _y = 0; 887 } 888 field()889 inline const Field & field() const { return *_field; } 890 891 public: 892 893 uint64_t _two_64; 894 895 private: 896 897 const Field *_field; 898 uint64_t _y; 899 }; 900 901 //! Specialization of DotProductDomain for uint32_t modular field 902 903 template<class Compute_t> 904 class DotProductDomain<Givaro::Modular<uint32_t,Compute_t> > : public VectorDomainBase<Givaro::Modular<uint32_t,Compute_t> > { 905 public: 906 907 typedef uint32_t Element; 908 typedef Givaro::Modular<uint32_t, Compute_t> Field; 909 DotProductDomain()910 DotProductDomain () {} DotProductDomain(const Field & F)911 DotProductDomain (const Field &F) : 912 VectorDomainBase<Field > (F) 913 {} 914 using VectorDomainBase<Field >::field; 915 using VectorDomainBase<Field >::faxpy; 916 917 protected: 918 template <class Vector1, class Vector2> dotSpecializedDD(Element & res,const Vector1 & v1,const Vector2 & v2)919 inline Element &dotSpecializedDD (Element &res, const Vector1 &v1, const Vector2 &v2) const 920 { 921 typename Vector1::const_iterator i; 922 typename Vector2::const_iterator j; 923 924 uint64_t y = 0; 925 uint64_t t; 926 927 for (i = v1.begin (), j = v2.begin (); i < v1.end (); ++i, ++j) { 928 t = (uint64_t) *i * (uint64_t) *j; 929 y += t; 930 931 if (y < t) 932 y += faxpy()._two_64; 933 } 934 935 y %= (uint64_t) field().characteristic(); 936 937 return res = (uint32_t) y; 938 } 939 940 941 template <class Vector1, class Vector2> dotSpecializedDSP(Element & res,const Vector1 & v1,const Vector2 & v2)942 inline Element &dotSpecializedDSP (Element &res, const Vector1 &v1, const Vector2 &v2) const 943 { 944 typename Vector1::first_type::const_iterator i_idx; 945 typename Vector1::second_type::const_iterator i_elt; 946 947 uint64_t y = 0; 948 uint64_t t = 0; 949 950 for (i_idx = v1.first.begin (), i_elt = v1.second.begin (); i_idx != v1.first.end (); ++i_idx, ++i_elt) { 951 t = (uint64_t) *i_elt * (uint64_t) v2[*i_idx]; 952 y += t; 953 if (y < t) 954 y += faxpy()._two_64; 955 } 956 957 y %= (uint64_t) field().characteristic(); 958 959 return res = (uint32_t)y; 960 } 961 962 963 }; 964 965 //! Specialization of MVProductDomain for uint32_t modular field 966 967 template <class Compute_t> 968 class MVProductDomain<Givaro::Modular<uint32_t,Compute_t> > { 969 public: 970 971 typedef uint32_t Element; 972 typedef Givaro::Modular<uint32_t,Compute_t> Field; 973 974 protected: 975 template <class Vector1, class Matrix, class Vector2> mulColDense(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v)976 inline Vector1 &mulColDense 977 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v) const 978 { 979 return mulColDenseSpecialized (VD, w, A, v, typename VectorTraits<typename Matrix::Column>::VectorCategory ()); 980 } 981 982 private: 983 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::DenseVectorTag)984 Vector1 &mulColDenseSpecialized 985 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 986 VectorCategories::DenseVectorTag) const 987 { 988 linbox_check (A.coldim () == v.size ()); 989 linbox_check (A.rowdim () == w.size ()); 990 991 typename Matrix::ConstColIterator i = A.colBegin (); 992 typename Vector2::const_iterator j; 993 typename Matrix::Column::const_iterator k; 994 std::vector<uint64_t>::iterator l; 995 996 uint64_t t; 997 998 if (_tmp.size () < w.size ()) 999 _tmp.resize (w.size ()); 1000 1001 std::fill (_tmp.begin (), _tmp.begin () +(ptrdiff_t) w.size (), 0); 1002 1003 for (j = v.begin (); j != v.end (); ++j, ++i) { 1004 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) { 1005 t = ((uint64_t) *k) * ((uint64_t) *j); 1006 1007 *l += t; 1008 1009 if (*l < t) 1010 *l += VD.faxpy()._two_64; 1011 } 1012 } 1013 1014 typename Vector1::iterator w_j; 1015 typedef typename Vector1::value_type element; 1016 1017 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 1018 *w_j = (element)(*l % VD.field ().characteristic()); 1019 1020 return w; 1021 } 1022 1023 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseSequenceVectorTag)1024 Vector1 &mulColDenseSpecialized 1025 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 1026 VectorCategories::SparseSequenceVectorTag) const 1027 { 1028 linbox_check (A.coldim () == v.size ()); 1029 linbox_check (A.rowdim () == w.size ()); 1030 1031 typename Matrix::ConstColIterator i = A.colBegin (); 1032 typename Vector2::const_iterator j; 1033 typename Matrix::Column::const_iterator k; 1034 std::vector<uint64_t>::iterator l; 1035 1036 uint64_t t; 1037 1038 if (_tmp.size () < w.size ()) 1039 _tmp.resize (w.size ()); 1040 1041 std::fill (_tmp.begin (), _tmp.begin () + (ptrdiff_t) w.size (), 0); 1042 1043 for (j = v.begin (); j != v.end (); ++j, ++i) { 1044 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) { 1045 t = ((uint64_t) k->second) * ((uint64_t) *j); 1046 1047 _tmp[k->first] += t; 1048 1049 if (_tmp[k->first] < t) 1050 _tmp[k->first] += VD.faxpy()._two_64; 1051 } 1052 } 1053 1054 typename Vector1::iterator w_j; 1055 typedef typename Vector1::value_type val_t; 1056 1057 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 1058 *w_j = val_t(*l % VD.field ().characteristic()); 1059 1060 return w; 1061 } 1062 1063 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseAssociativeVectorTag)1064 Vector1 &mulColDenseSpecialized 1065 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 1066 VectorCategories::SparseAssociativeVectorTag) const 1067 { 1068 linbox_check (A.coldim () == v.size ()); 1069 linbox_check (A.rowdim () == w.size ()); 1070 1071 typename Matrix::ConstColIterator i = A.colBegin (); 1072 typename Vector2::const_iterator j; 1073 typename Matrix::Column::const_iterator k; 1074 std::vector<uint64_t>::iterator l; 1075 1076 uint64_t t; 1077 1078 if (_tmp.size () < w.size ()) 1079 _tmp.resize (w.size ()); 1080 1081 std::fill (_tmp.begin (), _tmp.begin () +(ptrdiff_t) w.size (), 0); 1082 1083 for (j = v.begin (); j != v.end (); ++j, ++i) { 1084 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) { 1085 t = ((uint64_t) k->second) * ((uint64_t) *j); 1086 1087 _tmp[k->first] += t; 1088 1089 if (_tmp[k->first] < t) 1090 _tmp[k->first] += VD.faxpy()._two_64; 1091 } 1092 } 1093 1094 typename Vector1::iterator w_j; 1095 1096 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 1097 *w_j = (uint32_t) (uint32_t)*l % VD.field ().characteristic(); 1098 1099 return w; 1100 } 1101 1102 template <class Vector1, class Matrix, class Vector2> mulColDenseSpecialized(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v,VectorCategories::SparseParallelVectorTag)1103 Vector1 &mulColDenseSpecialized 1104 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 1105 VectorCategories::SparseParallelVectorTag) const 1106 { 1107 linbox_check (A.coldim () == v.size ()); 1108 linbox_check (A.rowdim () == w.size ()); 1109 1110 typename Matrix::ConstColIterator i = A.colBegin (); 1111 typename Vector2::const_iterator j; 1112 typename Matrix::Column::first_type::const_iterator k_idx; 1113 typename Matrix::Column::second_type::const_iterator k_elt; 1114 std::vector<uint64_t>::iterator l; 1115 1116 uint64_t t; 1117 1118 if (_tmp.size () < w.size ()) 1119 _tmp.resize (w.size ()); 1120 1121 std::fill (_tmp.begin (), _tmp.begin () +(ptrdiff_t) w.size (), 0); 1122 1123 for (j = v.begin (); j != v.end (); ++j, ++i) { 1124 for (k_idx = i->first.begin (), k_elt = i->second.begin (), l = _tmp.begin (); 1125 k_idx != i->first.end (); 1126 ++k_idx, ++k_elt, ++l) 1127 { 1128 t = ((uint64_t) *k_elt) * ((uint64_t) *j); 1129 1130 _tmp[*k_idx] += t; 1131 1132 if (_tmp[*k_idx] < t) 1133 _tmp[*k_idx] += VD.faxpy()._two_64; 1134 } 1135 } 1136 1137 typename Vector1::iterator w_j; 1138 typedef typename Vector1::value_type val_t; 1139 1140 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 1141 *w_j = val_t(*l % VD.field ().characteristic()); 1142 1143 return w; 1144 } 1145 1146 1147 mutable std::vector<uint64_t> _tmp; 1148 }; 1149 1150 } 1151 1152 1153 namespace LinBox { /* uint64_t */ 1154 1155 template<class Field> 1156 class DotProductDomain; 1157 template<class Field> 1158 class FieldAXPY; 1159 template<class Field> 1160 class MVProductDomain; 1161 1162 /*! Specialization of FieldAXPY for unsigned short modular field */ 1163 1164 template<typename Compute_t> 1165 class FieldAXPY<Givaro::Modular<uint64_t,Compute_t> > { 1166 public: 1167 1168 typedef uint64_t Element; 1169 typedef Givaro::Modular<uint64_t,Compute_t> Field; 1170 FieldAXPY(const Field & F)1171 FieldAXPY (const Field &F) : 1172 _field (&F), _y(0) 1173 { 1174 _two_64 = (uint64_t(1) << 32) % uint64_t(F.characteristic()); 1175 _two_64 = (_two_64 * _two_64) % uint64_t(F.characteristic()); 1176 } 1177 FieldAXPY(const FieldAXPY & faxpy)1178 FieldAXPY (const FieldAXPY &faxpy) : 1179 _two_64 (faxpy._two_64), _field (faxpy._field), _y (0) 1180 {} 1181 1182 FieldAXPY<Field > &operator = (const FieldAXPY &faxpy) 1183 { 1184 _field = faxpy._field; 1185 _y = faxpy._y; 1186 return *this; 1187 } 1188 mulacc(const Element & a,const Element & x)1189 inline uint64_t& mulacc (const Element &a, const Element &x) 1190 { 1191 uint64_t t = (uint64_t) a * (uint64_t) x; 1192 _y += t; 1193 1194 if (_y < t) 1195 return _y += _two_64; 1196 else 1197 return _y; 1198 } 1199 accumulate(const Element & t)1200 inline uint64_t& accumulate (const Element &t) 1201 { 1202 _y += t; 1203 1204 if (_y < t) 1205 return _y += _two_64; 1206 else 1207 return _y; 1208 } 1209 accumulate_special(const Element & t)1210 inline uint64_t& accumulate_special (const Element &t) 1211 { 1212 return _y += t; 1213 } 1214 get(Element & y)1215 inline Element &get (Element &y) const 1216 { 1217 const_cast<FieldAXPY<Field>*>(this)->_y %= (uint64_t) field().characteristic(); 1218 //if ((int64_t) _y < 0) const_cast<FieldAXPY<Field>*>(this)->_y += field().characteristic(); 1219 return y = (uint64_t) _y; 1220 } 1221 assign(const Element y)1222 inline FieldAXPY &assign (const Element y) 1223 { 1224 _y = y; 1225 return *this; 1226 } 1227 reset()1228 inline void reset() { 1229 _y = 0; 1230 } 1231 field()1232 inline const Field & field() const { return *_field; } 1233 1234 public: 1235 1236 uint64_t _two_64; 1237 1238 private: 1239 1240 const Field *_field; 1241 uint64_t _y; 1242 }; 1243 1244 //! Specialization of DotProductDomain for uint64_t modular field 1245 1246 template <typename Compute_t> 1247 class DotProductDomain<Givaro::Modular<uint64_t,Compute_t>> : public VectorDomainBase<Givaro::Modular<uint64_t,Compute_t> > { 1248 public: 1249 1250 typedef uint64_t Element; 1251 typedef Givaro::Modular<uint64_t,Compute_t> Field; 1252 DotProductDomain()1253 DotProductDomain () {} DotProductDomain(const Field & F)1254 DotProductDomain (const Field &F) : 1255 VectorDomainBase<Field > (F) 1256 {} 1257 using VectorDomainBase<Field >::field; 1258 using VectorDomainBase<Field >::faxpy; 1259 1260 protected: 1261 template <class Vector1, class Vector2> dotSpecializedDD(Element & res,const Vector1 & v1,const Vector2 & v2)1262 inline Element &dotSpecializedDD (Element &res, const Vector1 &v1, const Vector2 &v2) const 1263 { 1264 1265 typename Vector1::const_iterator i; 1266 typename Vector2::const_iterator j; 1267 1268 uint64_t y = 0; 1269 uint64_t t; 1270 1271 for (i = v1.begin (), j = v2.begin (); i < v1.end (); ++i, ++j) 1272 { 1273 t = ( (uint64_t) *i ) * ( (uint64_t) *j ); 1274 y += t; 1275 1276 if (y < t) 1277 y += faxpy()._two_64; 1278 } 1279 1280 y %= (uint64_t) field().characteristic(); 1281 return res = (Element)y; 1282 1283 } 1284 1285 template <class Vector1, class Vector2> dotSpecializedDSP(Element & res,const Vector1 & v1,const Vector2 & v2)1286 inline Element &dotSpecializedDSP (Element &res, const Vector1 &v1, const Vector2 &v2) const 1287 { 1288 typename Vector1::first_type::const_iterator i_idx; 1289 typename Vector1::second_type::const_iterator i_elt; 1290 1291 uint64_t y = 0; 1292 uint64_t t; 1293 1294 for (i_idx = v1.first.begin (), i_elt = v1.second.begin (); i_idx != v1.first.end (); ++i_idx, ++i_elt) 1295 { 1296 t = ( (uint64_t) *i_elt ) * ( (uint64_t) v2[*i_idx] ); 1297 y += t; 1298 1299 if (y < t) 1300 y += faxpy()._two_64; 1301 } 1302 1303 1304 y %= (uint64_t) field().characteristic(); 1305 1306 return res = (Element) y; 1307 } 1308 1309 1310 }; 1311 1312 //! Specialization of MVProductDomain for uint64_t modular field 1313 1314 template <typename Compute_t> 1315 class MVProductDomain<Givaro::Modular<uint64_t,Compute_t> > { 1316 public: 1317 1318 typedef uint64_t Element; 1319 typedef Givaro::Modular<uint64_t,Compute_t> Field; 1320 1321 protected: 1322 template <class Vector1, class Matrix, class Vector2> mulColDense(const VectorDomain<Field> & VD,Vector1 & w,const Matrix & A,const Vector2 & v)1323 inline Vector1 &mulColDense 1324 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v) const 1325 { 1326 return mulColDenseSpecialized (VD, w, A, v, typename VectorTraits<typename Matrix::Column>::VectorCategory ()); 1327 } 1328 1329 private: 1330 template <class Vector1, class Matrix, class Vector2> 1331 Vector1 &mulColDenseSpecialized 1332 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 1333 VectorCategories::DenseVectorTag) const; 1334 template <class Vector1, class Matrix, class Vector2> 1335 Vector1 &mulColDenseSpecialized 1336 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 1337 VectorCategories::SparseSequenceVectorTag) const; 1338 template <class Vector1, class Matrix, class Vector2> 1339 Vector1 &mulColDenseSpecialized 1340 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 1341 VectorCategories::SparseAssociativeVectorTag) const; 1342 template <class Vector1, class Matrix, class Vector2> 1343 Vector1 &mulColDenseSpecialized 1344 (const VectorDomain<Field > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 1345 VectorCategories::SparseParallelVectorTag) const; 1346 1347 mutable std::vector<uint64_t> _tmp; 1348 }; 1349 1350 } 1351 1352 1353 #endif // __LINBOX_field_modular_unsigned_H 1354 1355 // Local Variables: 1356 // mode: C++ 1357 // tab-width: 4 1358 // indent-tabs-mode: nil 1359 // c-basic-offset: 4 1360 // End: 1361 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 1362