1 /* Copyright (C) 2007 LinBox
2  * Written by JG Dumas
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 /*! @file algorithms/cra-builder-early-multip.h
25  * @ingroup algorithms
26  * @brief NO DOC
27  */
28 
29 #ifndef __LINBOX_cra_early_multip_H
30 #define __LINBOX_cra_early_multip_H
31 
32 #include "linbox/util/timer.h"
33 #include <stdlib.h>
34 #include "linbox/integer.h"
35 #include "linbox/solutions/methods.h"
36 #include <vector>
37 #include <utility>
38 
39 #include "linbox/algorithms/cra-builder-single.h"
40 #include "linbox/algorithms/cra-builder-full-multip.h"
41 
42 
43 namespace LinBox
44 {
45 
46 	/*!  @brief NO DOC
47 	 * @ingroup CRA
48 	 *
49 	 */
50 
51 	template<class Domain_Type>
52 	struct CRABuilderEarlyMultip : public CRABuilderEarlySingle<Domain_Type>, public CRABuilderFullMultip<Domain_Type> {
53 		typedef Domain_Type			Domain;
54 		typedef typename Domain::Element DomainElement;
55 		typedef CRABuilderEarlyMultip<Domain> 		Self_t;
56 
57 	protected:
58 		// Random coefficients for a linear combination
59 		// of the elements to be reconstructed
60 		std::vector< size_t >      	randv;
61 
resultCRABuilderEarlyMultip62 		Integer& result(Integer &d) { std::cout << "should not be called" << std::endl; return d ;} ; // DON'T TOUCH
63 	public:
64 
65 		CRABuilderEarlyMultip(const size_t EARLY=LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD) :
66 			CRABuilderEarlySingle<Domain>(EARLY), CRABuilderFullMultip<Domain>()
67 		{}
68 
getModulusCRABuilderEarlyMultip69 		Integer& getModulus(Integer& m)
70 		{
71 			CRABuilderEarlySingle<Domain>::getModulus(m);
72 			return m;
73 		}
getResidueCRABuilderEarlyMultip74 		Integer& getResidue(Integer& m)
75 		{
76 			CRABuilderEarlySingle<Domain>::getResidue(m);
77 			return m;
78 		}
79 
80 		template<template<class T> class Vect>
getResidueCRABuilderEarlyMultip81 		Vect<Integer>& getResidue(Vect<Integer>& m)
82 		{
83 			CRABuilderFullMultip<Domain>::getResidue(m);
84 			return m;
85 		}
86 
87 		//! Init
88 		template<template<class T> class Vect>
initializeCRABuilderEarlyMultip89 		void initialize (const Integer& D, const Vect<Integer>& e)
90 		{
91 			srand48(BaseTimer::seed());
92 			randv. resize ( e.size() );
93 			for ( std::vector<size_t>::iterator int_p = randv. begin(); int_p != randv. end(); ++ int_p)
94 				*int_p = ((size_t)lrand48()) % 20000;
95 			Integer z;
96 			dot(z, D, e, randv);
97 			CRABuilderEarlySingle<Domain>::initialize(D, z);
98 			CRABuilderFullMultip<Domain>::initialize(D, e);
99 		}
100 
101 		template<class Vect>
102                 //		template<template <class> class Alloc, template<class, class> class Vect>
initializeCRABuilderEarlyMultip103 		void initialize (const Domain& D, const Vect& e)
104 		{
105 			// Random coefficients for a linear combination
106 			// of the elements to be reconstructed
107 			srand48(BaseTimer::seed());
108 			randv. resize ( e.size() );
109 			for ( std::vector<size_t>::iterator int_p = randv. begin();
110 			      int_p != randv. end(); ++ int_p)
111 				*int_p = ((size_t)lrand48()) % 20000;
112 			DomainElement z;
113 			// Could be much faster
114 			// - do not compute twice the product of moduli
115 			// - reconstruct one element of e until Early Termination,
116 			//   then only, try a random linear combination.
117 			CRABuilderEarlySingle<Domain>::initialize (D, dot(z, D, e, randv));
118 			CRABuilderFullMultip<Domain>::initialize (D, e);
119 		}
120 
121 		template<class OKDomain>
initializeCRABuilderEarlyMultip122 		void initialize (const Domain& D, const BlasVector<OKDomain>& e)
123 		{
124 			// Random coefficients for a linear combination
125 			// of the elements to be reconstructed
126 			srand48(BaseTimer::seed());
127 			randv. resize ( e.size() );
128 			for ( std::vector<size_t>::iterator int_p = randv. begin();
129 			      int_p != randv. end(); ++ int_p)
130 				*int_p = ((size_t)lrand48()) % 20000;
131 			DomainElement z;
132 			// Could be much faster
133 			// - do not compute twice the product of moduli
134 			// - reconstruct one element of e until Early Termination,
135 			//   then only, try a random linear combination.
136 			CRABuilderEarlySingle<Domain>::initialize(D,dot(z, D, e, randv) );
137 			CRABuilderFullMultip<Domain>::initialize(D, e);
138 		}
139 
140 		//! Progress
141 		template<template<class T> class Vect>
progressCRABuilderEarlyMultip142 		void progress (const Integer& D, const Vect<Integer>& e)
143 		{
144 
145 			Integer z;
146 			CRABuilderEarlySingle<Domain>::progress(D, dot(z, D, e, randv));
147 			CRABuilderFullMultip<Domain>::progress(D, e);
148 		}
149 
150 #if 1
151 		template<class Vect>
152                 //		template<template <class> class Alloc, template<class, class> class Vect>
progressCRABuilderEarlyMultip153 		void progress (const Domain& D, const Vect& e)
154 		{
155 			// DomainElement z;
156 			/*!@todo Could be much faster
157 			  - do not compute twice the product of moduli
158 			  - reconstruct one element of e until Early Termination,
159 			  then only, try a random linear combination.
160 			*/
161                         DomainElement z;
162                         CRABuilderEarlySingle<Domain>::progress(D, dot(z, D, e, randv));
163                         CRABuilderFullMultip<Domain>::progress(D, e);
164 		}
165 #endif
166 
167 		template<class OKDomain>
progressCRABuilderEarlyMultip168 		void progress (const Domain& D, const BlasVector<OKDomain>& e)
169 		{
170 			DomainElement z;
171 			/*!@todo Could be much faster
172 			  - do not compute twice the product of moduli
173 			  - reconstruct one element of e until Early Termination,
174 			  then only, try a random linear combination.
175 			*/
176 			CRABuilderEarlySingle<Domain>::progress(D, dot(z, D, e, randv));
177 			CRABuilderFullMultip<Domain>::progress(D, e);
178 		}
179 
180 		//! Result
181 //		template<template <class> class Alloc, template<class, class> class Vect>
182 		template<class Vect>
resultCRABuilderEarlyMultip183                 Vect& result(Vect& d)
184 		{
185 			return CRABuilderFullMultip<Domain>::result(d);
186 		}
187 
resultCRABuilderEarlyMultip188 		BlasVector<Givaro::ZRing<Integer> >& result(BlasVector<Givaro::ZRing<Integer> >& d)
189 		{
190 			return CRABuilderFullMultip<Domain>::result(d);
191 		}
192 
193 		//! terminate
terminatedCRABuilderEarlyMultip194 		bool terminated()
195 		{
196 			return CRABuilderEarlySingle<Domain>::terminated();
197 		}
198 
noncoprimeCRABuilderEarlyMultip199 		bool noncoprime(const Integer& i) const
200 		{
201 			return CRABuilderEarlySingle<Domain>::noncoprime(i);
202 		}
203 
changeVectorCRABuilderEarlyMultip204 		bool changeVector()
205 		{
206 			for ( std::vector<size_t>::iterator int_p = randv. begin();int_p != randv. end(); ++ int_p)
207 				*int_p = ((size_t)lrand48()) % 20000;
208 
209 			std::vector<Integer> e(randv.size());
210 			/* clear CRAEarlySingle; */
211 			CRABuilderEarlySingle<Domain>::occurency_ = 0;
212 			CRABuilderEarlySingle<Domain>::nextM_ = 1UL;
213 			CRABuilderEarlySingle<Domain>::primeProd_ = 1UL;
214 			CRABuilderEarlySingle<Domain>::residue_ = 0;
215 
216 			/* Computation of residue_ */
217             for (auto it = CRABuilderFullMultip<Domain>::shelves_begin();
218                  it != CRABuilderFullMultip<Domain>::shelves_end();
219                  ++it)
220             {
221                 if (it->occupied) {
222 					Integer D = it->mod();
223 					std::vector<Integer> e_v(randv.size());
224 					e_v = it->residue;
225 					Integer z;
226 					dot(z,D, e_v, randv);
227 					Integer prev_residue_ = CRABuilderEarlySingle<Domain>::residue_;
228 					CRABuilderEarlySingle<Domain>::progress(D,z);
229 					if (prev_residue_ == CRABuilderEarlySingle<Domain>::residue_ )
230 						CRABuilderEarlySingle<Domain>::occurency_ += it->count;
231 					if ( CRABuilderEarlySingle<Domain>::terminated() ) {
232 						return true;
233 					}
234                 }
235 			}
236 			return false;
237 		}
238 
239 	protected:
240 
241 		/*! @bug why a dot product here ?
242 		 */
243 		template <class Vect1, class Vect2>
dotCRABuilderEarlyMultip244 		Integer& dot (Integer& z, const Integer& D, const Vect1& v1, const Vect2& v2)
245 		{
246 			z = 0;
247 			typename Vect1::const_iterator v1_p;
248 			typename Vect2::const_iterator v2_p;
249 			for (v1_p  = v1. begin(), v2_p = v2. begin(); v1_p != v1. end(); ++ v1_p, ++ v2_p) {
250 				z = (z + (*v1_p)*(*v2_p))%D;
251 			}
252 			return z;
253 		}
254 
255 		/*! @bug why a dot product here ?
256 		 */
257 		template <class Vect1, class Vect2>
dotCRABuilderEarlyMultip258 		DomainElement& dot (DomainElement& z, const Domain& D,
259 				    const Vect1& v1,
260 				    const Vect2& v2)
261 		{
262 
263 			D.assign(z,D.zero); DomainElement tmp;
264 			typename Vect1::const_iterator v1_p;
265 			typename Vect2::const_iterator v2_p;
266 			for (v1_p  = v1. begin(), v2_p = v2. begin();
267 			     v1_p != v1. end();
268 			     ++ v1_p, ++ v2_p)
269 				D.axpyin(z, (*v1_p), D.init(tmp, (*v2_p)));
270 
271 			//             commentator().report(Commentator::LEVEL_ALWAYS, INTERNAL_DESCRIPTION) << "v: " << v2 << std::endl;
272 			//             commentator().report(Commentator::LEVEL_ALWAYS, INTERNAL_DESCRIPTION) << "z: " << z << std::endl;
273 			return z;
274 		}
275 
276 		template <class Vect2, class OKDomain>
dotCRABuilderEarlyMultip277 		DomainElement& dot (DomainElement& z, const Domain& D,
278 				    const BlasVector<OKDomain>& v1,
279 				    const Vect2& v2)
280 		{
281 
282 			D.assign(z,D.zero); DomainElement tmp;
283 			typename BlasVector<Domain>::const_iterator v1_p;
284 			typename Vect2::const_iterator v2_p;
285 			for (v1_p  = v1. begin(), v2_p = v2. begin();
286 			     v1_p != v1. end();
287 			     ++ v1_p, ++ v2_p)
288 				D.axpyin(z, (*v1_p), D.init(tmp, (*v2_p)));
289 
290 			//             commentator().report(Commentator::LEVEL_ALWAYS, INTERNAL_DESCRIPTION) << "v: " << v2 << std::endl;
291 			//             commentator().report(Commentator::LEVEL_ALWAYS, INTERNAL_DESCRIPTION) << "z: " << z << std::endl;
292 			return z;
293 		}
294 
295 	};
296 }
297 #endif //__LINBOX_cra_early_multip_H
298 
299 // Local Variables:
300 // mode: C++
301 // tab-width: 4
302 // indent-tabs-mode: nil
303 // c-basic-offset: 4
304 // End:
305 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
306