1 /* Copyright (C) 2010 LinBox 2 * 3 * 4 * 5 * ========LICENCE======== 6 * This file is part of the library LinBox. 7 * 8 * LinBox is free software: you can redistribute it and/or modify 9 * it under the terms of the GNU Lesser General Public 10 * License as published by the Free Software Foundation; either 11 * version 2.1 of the License, or (at your option) any later version. 12 * 13 * This library is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 16 * Lesser General Public License for more details. 17 * 18 * You should have received a copy of the GNU Lesser General Public 19 * License along with this library; if not, write to the Free Software 20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 21 * ========LICENCE======== 22 */ 23 24 25 #ifndef __LINBOX_pir_modular_int32_H 26 #define __LINBOX_pir_modular_int32_H 27 28 #include <givaro/modular-integral.h> 29 //#include <linbox/util/debug.h> 30 #include <linbox/vector/vector-domain.h> 31 32 //#include "linbox/ring/modular.h" 33 #ifndef LINBOX_MAX_INT 34 #define LINBOX_MAX_INT 2147483647 35 #endif 36 37 #ifndef LINBOX_MAX_MODULUS 38 #define LINBOX_MAX_MODULUS 1073741824 39 #endif 40 #include "linbox/field/field-traits.h" 41 42 // Namespace in which all LinBox code resides 43 namespace LinBox 44 { 45 46 template< class PIR> 47 class PIRModular; 48 49 template< class Element > 50 class ModularRandIter; 51 52 template<class Field> 53 class DotProductDomain; 54 55 template<class Field> 56 class FieldAXPY; 57 58 template<class Field> 59 class MVProductDomain; 60 61 template <class Ring> 62 struct ClassifyRing; 63 64 template <class Element> 65 struct ClassifyRing<PIRModular<Element> >; 66 67 template <> 68 struct ClassifyRing<PIRModular<int32_t> > { 69 typedef RingCategories::ModularTag categoryTag; 70 }; 71 72 /// \ingroup ring 73 template <> 74 class PIRModular<int32_t> : public Givaro::Modular<int32_t> { 75 76 public: 77 78 friend class FieldAXPY<PIRModular<int32_t> >; 79 80 friend class DotProductDomain<PIRModular<int32_t> >; 81 82 friend class MVProductDomain<PIRModular<int32_t> >; 83 84 typedef int32_t Element; 85 86 typedef Givaro::Modular<int32_t>::RandIter RandIter; 87 88 uint64_t _two_64; 89 90 91 //default modular field,taking 65521 as default modulus 92 PIRModular () : 93 Givaro::Modular<int32_t>(65521) 94 { 95 _two_64 = (uint64_t(1) << 32) % uint64_t(this->characteristic()); 96 _two_64 = (_two_64 * _two_64) % uint64_t(this->characteristic()); 97 98 } 99 100 PIRModular (int32_t value, int32_t exp = 1) : 101 Givaro::Modular<int32_t>(value) 102 { 103 _two_64 = (uint64_t(1) << 32) % uint64_t(this->characteristic()); 104 _two_64 = (_two_64 * _two_64) % uint64_t(this->characteristic()); 105 106 } 107 108 109 /** PIR functions, gcd, xgcd, dxgcd */ 110 111 Element& gcd (Element& g, const Element& a, const Element& b) const 112 { 113 114 GCD (g, a, b); 115 116 return g; 117 118 } 119 120 Element& xgcd (Element& g, Element& s, Element& t, const Element& a, const Element& b) const 121 { 122 123 XGCD (g, s, t, a, b); 124 125 if (s < 0) 126 s += _p; 127 128 if (t < 0) 129 t += _p; 130 131 132 return g; 133 } 134 135 Element& dxgcd (Element& g, Element& s, Element& t, Element& a1, Element& b1, 136 const Element& a, const Element& b) const 137 { 138 139 140 xgcd (g, s, t, a, b); 141 142 if (g != 0) { 143 144 a1 = a / g; 145 146 b1 = b / g; 147 } 148 149 else { 150 151 a1 = s; 152 153 b1 = t; 154 } 155 156 157 return g; 158 159 } 160 161 bool isDivisor (const Element& a, const Element& b) const 162 { 163 164 Element g; 165 166 if (a == 0) return false; 167 168 else if (b == 0) return true; 169 170 else { 171 172 gcd (g, a, _p); 173 174 return (b % g) == 0; 175 } 176 177 } 178 179 Element& div (Element& d, const Element& a, const Element& b) const 180 { 181 182 Element g, s; 183 184 HXGCD (g, s, b, _p); 185 186 Element r; 187 188 r = a % g; 189 190 if (r != 0) throw PreconditionFailed(LB_FILE_LOC,"Div: not divisible"); 191 192 else { 193 194 d = a / g; 195 196 mulin (d, s); 197 } 198 199 return d; 200 201 } 202 203 Element& normal (Element& a, const Element& b) const 204 { 205 206 if (b == 0) return a = 0; 207 else { 208 GCD (a, b, _p); 209 210 return a; 211 } 212 } 213 214 215 Element& gcdin (Element& a, const Element& b) const 216 { 217 218 GCD (a, a, b); 219 220 221 return normalIn(a); // is this efficient? 222 } 223 224 Element& normalIn (Element& a) const 225 { 226 if (a == 0) return a; 227 else { 228 GCD (a, a, _p); 229 230 return a; 231 } 232 233 } 234 235 236 Element& divin (Element& a, const Element& b) const 237 { 238 239 div (a, a, b); 240 241 return a; 242 } 243 244 245 bool isUnit(const Element& a) const 246 { 247 248 Element g; 249 250 GCD (g, a, _p); 251 252 253 // std::cout << a << " is a unit or not " << g; 254 255 // std::cout << "modulus = " << _p <<"\n"; 256 257 return g == 1; 258 259 } 260 261 private: 262 static void GCD (int32_t& g, int32_t a, int32_t b) { 263 264 int32_t u, v, /* q,*/ r; 265 266 if (a < 0) { 267 if (a < -LINBOX_MAX_INT) throw PreconditionFailed(LB_FILE_LOC,"XGCD: integer overflow"); 268 a = -a; 269 270 } 271 272 if (b < 0) { 273 if (b < -LINBOX_MAX_INT) throw PreconditionFailed(LB_FILE_LOC,"XGCD: integer overflow"); 274 b = -b; 275 } 276 277 u = a; v = b; 278 279 while (v != 0) { 280 // q = u / v; 281 r = u % v; 282 u = v; 283 v = r; 284 } 285 286 g = u; 287 288 } 289 290 static void XGCD(int32_t& d, int32_t& s, int32_t& t, int32_t a, int32_t b) { 291 int32_t u, v, u0, v0, u1, v1, u2, v2, q, r; 292 293 int32_t aneg = 0, bneg = 0; 294 295 if (a < 0) { 296 if (a < -LINBOX_MAX_INT) throw PreconditionFailed(LB_FILE_LOC,"XGCD: integer overflow"); 297 a = -a; 298 aneg = 1; 299 } 300 301 if (b < 0) { 302 if (b < -LINBOX_MAX_INT) throw PreconditionFailed(LB_FILE_LOC,"XGCD: integer overflow"); 303 b = -b; 304 bneg = 1; 305 } 306 307 u1 = 1; v1 = 0; 308 u2 = 0; v2 = 1; 309 u = a; v = b; 310 311 while (v != 0) { 312 q = u / v; 313 r = u % v; 314 u = v; 315 v = r; 316 u0 = u2; 317 v0 = v2; 318 u2 = u1 - q*u2; 319 v2 = v1- q*v2; 320 u1 = u0; 321 v1 = v0; 322 } 323 324 if (aneg) 325 u1 = -u1; 326 327 if (bneg) 328 v1 = -v1; 329 330 d = u; 331 s = u1; 332 t = v1; 333 } 334 335 336 static void HXGCD (int32_t& d, int32_t& s, int32_t a, int32_t b) { 337 338 int32_t u, v, u0, u1, u2, q, r; 339 340 int32_t aneg = 0; 341 342 if (a < 0) { 343 if (a < -LINBOX_MAX_INT) throw PreconditionFailed(LB_FILE_LOC,"XGCD: integer overflow"); 344 a = -a; 345 aneg = 1; 346 } 347 348 if (b < 0) { 349 if (b < -LINBOX_MAX_INT) throw PreconditionFailed(LB_FILE_LOC,"XGCD: integer overflow"); 350 b = -b; 351 } 352 353 u1 = 1; 354 u2 = 0; 355 u = a; v = b; 356 357 while (v != 0) { 358 q = u / v; 359 r = u % v; 360 u = v; 361 v = r; 362 u0 = u2; 363 364 u2 = u1 - q*u2; 365 366 u1 = u0; 367 368 } 369 370 if (aneg) 371 u1 = -u1; 372 373 374 d = u; 375 s = u1; 376 377 } 378 379 }; 380 381 template <> 382 class FieldAXPY<PIRModular<int32_t> > { 383 public: 384 385 typedef int32_t Element; 386 typedef int64_t Abnormal; 387 typedef PIRModular<int32_t> Field; 388 389 FieldAXPY (const Field &F) : 390 _field (F),_y(0) 391 {} 392 393 394 FieldAXPY (const FieldAXPY &faxpy) : 395 _field (faxpy._field), _y (0) 396 {} 397 398 // FieldAXPY<PIRModular<int32_t> > &operator = (const FieldAXPY &faxpy) { 399 // _field = faxpy._field; 400 // _y = faxpy._y; 401 // return *this; 402 // } 403 404 inline uint64_t& mulacc (const Element &a, const Element &x) { 405 uint64_t t = (uint64_t) a * (uint64_t) x; 406 _y += t; 407 if (_y < t) 408 return _y += field()._two_64; 409 else 410 return _y; 411 } 412 413 inline uint64_t& accumulate (const Element &t) { 414 _y += (uint64_t)t; 415 if (_y < (uint64_t)t) 416 return _y += field()._two_64; 417 else 418 return _y; 419 } 420 421 inline Element& get (Element &y) { 422 y = Element(_y % (uint64_t) field().characteristic()); 423 return y; 424 } 425 426 inline FieldAXPY &assign (const Element y) { 427 _y = (uint64_t)y; 428 return *this; 429 } 430 431 inline void reset() { 432 _y = 0; 433 } 434 435 inline const Field & field() const { return _field; } 436 437 protected: 438 const Field &_field; 439 uint64_t _y; 440 }; 441 442 443 template <> 444 class DotProductDomain<PIRModular<int32_t> > : public VectorDomainBase<PIRModular<int32_t> > { 445 446 public: 447 typedef int32_t Element; 448 DotProductDomain (const PIRModular<int32_t> &F) : 449 VectorDomainBase<PIRModular<int32_t> > (F) 450 {} 451 452 453 protected: 454 template <class Vector1, class Vector2> 455 inline Element &dotSpecializedDD (Element &res, const Vector1 &v1, const Vector2 &v2) const 456 { 457 458 typename Vector1::const_iterator i; 459 typename Vector2::const_iterator j; 460 461 uint64_t y = 0; 462 uint64_t t; 463 464 for (i = v1.begin (), j = v2.begin (); i < v1.end (); ++i, ++j) { 465 t = ( (uint64_t) *i ) * ( (uint64_t) *j ); 466 y += t; 467 468 if (y < t) 469 y += field()._two_64; 470 } 471 472 y %= (uint64_t) field().characteristic(); 473 return res = Element(y); 474 475 } 476 477 template <class Vector1, class Vector2> 478 inline Element &dotSpecializedDSP (Element &res, const Vector1 &v1, const Vector2 &v2) const 479 { 480 typename Vector1::first_type::const_iterator i_idx; 481 typename Vector1::second_type::const_iterator i_elt; 482 483 uint64_t y = 0; 484 uint64_t t; 485 486 for (i_idx = v1.first.begin (), i_elt = v1.second.begin (); i_idx != v1.first.end (); ++i_idx, ++i_elt) { 487 t = ( (uint64_t) *i_elt ) * ( (uint64_t) v2[*i_idx] ); 488 y += t; 489 490 if (y < t) 491 y += field()._two_64; 492 } 493 494 495 y %= (uint64_t) field().characteristic(); 496 497 return res = (Element)y; 498 } 499 }; 500 501 502 // Specialization of MVProductDomain for int32_t modular field 503 template <> 504 class MVProductDomain<PIRModular<int32_t> > { 505 public: 506 507 typedef int32_t Element; 508 509 protected: 510 template <class Vector1, class Matrix, class Vector2> 511 inline Vector1 &mulColDense 512 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v) const 513 { 514 return mulColDenseSpecialized 515 (VD, w, A, v, typename VectorTraits<typename Matrix::Column>::VectorCategory ()); 516 } 517 518 private: 519 template <class Vector1, class Matrix, class Vector2> 520 Vector1 &mulColDenseSpecialized 521 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 522 VectorCategories::DenseVectorTag) const; 523 template <class Vector1, class Matrix, class Vector2> 524 Vector1 &mulColDenseSpecialized 525 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 526 VectorCategories::SparseSequenceVectorTag) const; 527 template <class Vector1, class Matrix, class Vector2> 528 Vector1 &mulColDenseSpecialized 529 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 530 VectorCategories::SparseAssociativeVectorTag) const; 531 template <class Vector1, class Matrix, class Vector2> 532 Vector1 &mulColDenseSpecialized 533 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 534 VectorCategories::SparseParallelVectorTag) const; 535 536 mutable std::vector<uint64_t> _tmp; 537 }; 538 539 template <class Vector1, class Matrix, class Vector2> 540 Vector1 &MVProductDomain<PIRModular<int32_t> >::mulColDenseSpecialized 541 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 542 VectorCategories::DenseVectorTag) const 543 { 544 545 linbox_check (A.coldim () == v.size ()); 546 linbox_check (A.rowdim () == w.size ()); 547 548 typename Matrix::ConstColIterator i = A.colBegin (); 549 typename Vector2::const_iterator j; 550 typename Matrix::Column::const_iterator k; 551 std::vector<uint64_t>::iterator l; 552 553 uint64_t t; 554 555 if (_tmp.size () < w.size ()) 556 _tmp.resize (w.size ()); 557 558 std::fill (_tmp.begin (), _tmp.begin () + w.size (), 0); 559 560 for (j = v.begin (); j != v.end (); ++j, ++i) { 561 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) { 562 t = ((uint64_t) *k) * ((uint64_t) *j); 563 564 *l += t; 565 566 if (*l < t) 567 *l += VD.faxpy().field()._two_64; 568 } 569 } 570 571 typename Vector1::iterator w_j; 572 573 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 574 *w_j = *l % VD.field().characteristic(); 575 576 return w; 577 } 578 579 template <class Vector1, class Matrix, class Vector2> 580 Vector1 &MVProductDomain<PIRModular<int32_t> >::mulColDenseSpecialized 581 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 582 VectorCategories::SparseSequenceVectorTag) const 583 { 584 linbox_check (A.coldim () == v.size ()); 585 linbox_check (A.rowdim () == w.size ()); 586 587 typename Matrix::ConstColIterator i = A.colBegin (); 588 typename Vector2::const_iterator j; 589 typename Matrix::Column::const_iterator k; 590 std::vector<uint64_t>::iterator l; 591 592 uint64_t t; 593 594 if (_tmp.size () < w.size ()) 595 _tmp.resize (w.size ()); 596 597 std::fill (_tmp.begin (), _tmp.begin () + w.size (), 0); 598 599 for (j = v.begin (); j != v.end (); ++j, ++i) { 600 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) { 601 t = ((uint64_t) k->second) * ((uint64_t) *j); 602 603 _tmp[k->first] += t; 604 605 if (_tmp[k->first] < t) 606 _tmp[k->first] += VD.faxpy().field()._two_64; 607 } 608 } 609 610 typename Vector1::iterator w_j; 611 612 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 613 *w_j = *l % VD.field().characteristic(); 614 615 return w; 616 } 617 618 template <class Vector1, class Matrix, class Vector2> 619 Vector1 &MVProductDomain<PIRModular<int32_t> >::mulColDenseSpecialized 620 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 621 VectorCategories::SparseAssociativeVectorTag) const 622 { 623 624 linbox_check (A.coldim () == v.size ()); 625 linbox_check (A.rowdim () == w.size ()); 626 627 typename Matrix::ConstColIterator i = A.colBegin (); 628 typename Vector2::const_iterator j; 629 typename Matrix::Column::const_iterator k; 630 std::vector<uint64_t>::iterator l; 631 632 uint64_t t; 633 634 if (_tmp.size () < w.size ()) 635 _tmp.resize (w.size ()); 636 637 std::fill (_tmp.begin (), _tmp.begin () + w.size (), 0); 638 639 for (j = v.begin (); j != v.end (); ++j, ++i) { 640 for (k = i->begin (), l = _tmp.begin (); k != i->end (); ++k, ++l) { 641 t = ((uint64_t) k->second) * ((uint64_t) *j); 642 643 _tmp[k->first] += t; 644 645 if (_tmp[k->first] < t) 646 _tmp[k->first] += VD.faxpy().field()._two_64; 647 } 648 } 649 650 typename Vector1::iterator w_j; 651 652 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 653 *w_j = *l % VD.field().characteristic(); 654 655 return w; 656 } 657 658 template <class Vector1, class Matrix, class Vector2> 659 Vector1 &MVProductDomain<PIRModular<int32_t> >::mulColDenseSpecialized 660 (const VectorDomain<PIRModular<int32_t> > &VD, Vector1 &w, const Matrix &A, const Vector2 &v, 661 VectorCategories::SparseParallelVectorTag) const 662 { 663 664 linbox_check (A.coldim () == v.size ()); 665 linbox_check (A.rowdim () == w.size ()); 666 667 typename Matrix::ConstColIterator i = A.colBegin (); 668 typename Vector2::const_iterator j; 669 typename Matrix::Column::first_type::const_iterator k_idx; 670 typename Matrix::Column::second_type::const_iterator k_elt; 671 std::vector<uint64_t>::iterator l; 672 673 uint64_t t; 674 675 if (_tmp.size () < w.size ()) 676 _tmp.resize (w.size ()); 677 678 std::fill (_tmp.begin (), _tmp.begin () + w.size (), 0); 679 680 for (j = v.begin (); j != v.end (); ++j, ++i) { 681 for (k_idx = i->first.begin (), k_elt = i->second.begin (), l = _tmp.begin (); 682 k_idx != i->first.end (); 683 ++k_idx, ++k_elt, ++l) 684 { 685 t = ((uint64_t) *k_elt) * ((uint64_t) *j); 686 687 _tmp[*k_idx] += t; 688 689 if (_tmp[*k_idx] < t) 690 _tmp[*k_idx] += VD.faxpy().field()._two_64; 691 } 692 } 693 694 typename Vector1::iterator w_j; 695 696 for (w_j = w.begin (), l = _tmp.begin (); w_j != w.end (); ++w_j, ++l) 697 *w_j = *l % VD.field().characteristic(); 698 699 return w; 700 } 701 702 703 } 704 705 706 #endif //__LINBOX_pir_modular_int32_H 707 708 // Local Variables: 709 // mode: C++ 710 // tab-width: 4 711 // indent-tabs-mode: nil 712 // c-basic-offset: 4 713 // End: 714 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 715