This file is indexed.

/usr/include/linbox/algorithms/rational-solver-sn.inl is in liblinbox-dev 1.3.2-1.1build2.

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
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
/* Copyright (C) 2011 LinBox
 * Written Bryan Youse <>
 *
 *
 *
 * ========LICENCE========
 * This file is part of the library LinBox.
 *
  * LinBox is free software: you can redistribute it and/or modify
 * it under the terms of the  GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 * ========LICENCE========
 */

// rational-solver-sn.inl
// WARNING: this file is included from INSIDE the definition of class LinBox::RationalSolverSN

#ifdef DEBUGRC
#define debug(x) std::cerr << x << std::endl
#define debugneol(x) std::cerr << x
#define debug2(x, y) debug(x); y
#else
#define debug(x)
#define debug2(x, y)
#define debugneol(x)
#endif

#define log_2(x) (log(x)/M_LN2)

int zw_shift(NumericSolver & NS_S, size_t n, FVector &r, FVector &x)
{
	//  ZW method for calculating shift
	// compute ax
	FVector ax(n);
	NS_S.apply(ax, x);
	// compute ax = ax -r, the negative of residual
	for(size_t i=0; i<n; i++)
		_field.sub(ax[i], ax[i], r[i]);
	// compute possible shift
	int zw_shift_loc;
	double normr1, normr2, normr3, shift1, shift2;
	normr1 = zw_dmax(n, &*r.begin(), 1);
	normr2 = zw_dmax(n, &*ax.begin(), 1);
	normr3 = zw_dmax(n, &*x.begin(), 1);

	//cerr << normr1 << " " << normr2 << " " << normr3 << endl;
	//try to find a good scalar
	if (normr2 <.0000000001)
		zw_shift_loc = 30;
	else {
		shift1 = floor(log_2 (normr1 / normr2)) - 2;
		zw_shift_loc = (int)(30 < shift1 ? 30 : shift1);
	}

	normr3 = normr3 > 2 ? normr3 : 2;
	shift2 = floor(53. / log_2 (normr3));
	zw_shift_loc = (int)(zw_shift_loc < shift2 ? zw_shift_loc : shift2);

	return zw_shift_loc;
}

template <class IMatrix>
int rat_sol(IVector& numx, Int& denx, FVector& xs_int, FVector& xs_frac, IVector &b, IVector &lastb, FVector& r, FVector& lastr, FVector& x, integer &loopBound, IMatrix &IM) {

	debug("Matrix norm: " << mnorm);
	int thresh_expt = 1;  // consider as low as 1 or 2 (at least when n large)
	double threshold = 1.0/(1 << thresh_expt);

	size_t n = r.size();
	// xs_int: integer part of scaled x.
	// xs_frac: fractional part of scaled x.
	// x * 2^shift = xs_int + xs_frac.
	FVector nextx(n), quo(n);

	integer denx_i;
	typename Field::Element one; _field.init(one, 1);

	//  need to save original r for zw_shift calculation
	//  TODO: I took out the ZWSHIFT, still need last r??

	if(denx == 1){
		// compute first approximate solution x
		_numsolver.solve(x, r);

		//writeVec(x, "x from first solve");
		copy(r.begin(), r.end(), lastr.begin());
		copy(b.begin(), b.end(), lastb.begin());
		//-- compute xs_int, xs_frac, r (residue)
		update_xs(xs_int, xs_frac, x);
		if(exact_apply) update_r_exact(b, r, xs_int, IM);
		else update_r(r, xs_int);
	}

	//std::cerr << "while loopbound " << loopBound << std::endl;
	while (_ring.compare(denx, loopBound) < 0) {
		++iterations;
		//  x = DM^{-1}*r
		_numsolver.solve(nextx, r);

		for(size_t i=0; i<n; i++){
			// TODO - analyze logic here
			//  quo[i] = xs_frac[i] / nextx[i] - 1;
			//_field.div(quo[i], xs_frac[i], nextx[i]);
			//_field.subin(quo[i], one);
			_field.sub(quo[i], xs_frac[i], nextx[i]);
		}

		double q = zw_dmax((int)n, &*quo.begin(), 1);

		//writeVec(nextx, "nextx from loop");  //  DEBUG PRINTOUTS
		//writeVec(xs_frac, "xs_frac from loop");
		//writeVec(quo, "quo");
		/* fails to write!!
		   _VDF.write(std::cerr << "nextx from loop: ", nextx) << std::endl;
		   _VDF.write(std::cerr << "xs_frac from loop: ", xs_frac) << std::endl;
		   _VDF.write(std::cerr << "quo minus 1: ", quo) << std::endl;
		   if ( 0 == iterations%500) std::cerr << iterations << " MAX DIFFERENCE: " << q << std::endl;
		   */

		//if (q == 0.0)   (QUIT HERE??)
		if (q < threshold) {
			HIT++;
			// update numx and denx
			denx <<= shift;
			update_num (numx, xs_int);

			if(_VDF.isZero(r))
				return 2;

			// consider increasing the shift for next iteration
			upshift();

			// make x = nextx, then compute new residual, xs_int, xs_frac
			swap(x, nextx);
			copy(r.begin(), r.end(), lastr.begin());
			copy(b.begin(), b.end(), lastb.begin());
			update_xs(xs_int, xs_frac, x);
			if(exact_apply) update_r_exact(b, r, xs_int, IM);
			else update_r(r, xs_int);
		}
		else {  //  q >= threshold, back-off
			MISS++;

			//  point of no return, quit
			if (shift < 2){
				/*
				   std::cerr << "rat_sol failure, no bits in x overlap in nextx." << std::endl;
				   std::cerr << "Iterations: " << iterations << std::endl;
				   writeVec(b, "b", 0, 5); writeVec(r, "r", 0, 5);
				   */
				return -1;
			}

			downshift();
			//-- compute xs_int, xs_frac, r (residue)
			// but use last r as input
			copy(lastr.begin(), lastr.end(), r.begin());
			copy(lastb.begin(), lastb.end(), b.begin());
			update_xs(xs_int, xs_frac, x);
			if(exact_apply) update_r_exact(b, r, xs_int, IM);
			else update_r(r, xs_int);
		}
	}
	return 0;
}// rat_sol

inline void upshift()
{
	switch(sstatus){
		//  exponential increase
	case SHIFT_GROW:
		shift_prev = shift;
		debugneol("G");
		shift_max = shift<<=1;
		break;
	case SHIFT_SEARCH:
		shift_prev = shift;
		debugneol("S");
		shift = (shift + shift_max)>>1;
		searchPeak = true;
		break;
	case SHIFT_PEAK:
		debugneol("P");
		//  maybe increase if we have been successful for a while
		break;
	case SHIFT_MAX:
		debugneol("M");
		//  machine precision-- can go no higher
		break;
	case SHIFT_SHRINK:
		debugneol("H");
		break;
	}
	if(shift > SHIFT_BOUND){
		shift_max = shift = SHIFT_BOUND;
		sstatus = SHIFT_MAX;
	}
	debug("^shift: "  << shift << "  max: " << shift_max << "  prev: " << shift_prev);
}

inline void downshift()
{
	/*
	   shift -= 2;
	   sstatus = SHIFT_PEAK;
	   cerr << "peaked at " << shift << endl;
	   */
	//  back up
	switch(sstatus){
	case SHIFT_GROW:
		debugneol("G");
	case SHIFT_MAX:
		debugneol("M");
	case SHIFT_SEARCH:
		//  TODO - previous shift could fail
		debugneol("S");
		shift_max = shift;
		shift = (shift_prev + shift)>>1;
		sstatus = SHIFT_SEARCH;
		//  searchPeak true means we were going up but got knocked back
		if(shift == shift_prev || searchPeak)
			sstatus = SHIFT_PEAK;
		break;
	case SHIFT_SHRINK:
		debugneol("H");
		shift >>= 1;
		break;
	case SHIFT_PEAK:
		debugneol("P");
		shift -= 1;
		break;
	default:
		break;
	}

	debug("vshift: "  << shift << "  max: " << shift_max << "  prev: " << shift_prev);
} // downshift

inline void update_xs(FVector& xs_int, FVector& xs_frac, FVector& x)
{
	Float scalar, tmp;
	int64_t shifted = ((int64_t)1 << shift);
	_field.init(scalar, (double) shifted);

	//  make xs_int and xs_frac such that x*scalar = xs_int + xs_frac.
	for(size_t i = 0; i < xs_int.size(); ++i){
		// TODO: tmp can overflow a double
		tmp = x[i]*scalar;
		xs_int[i] = floor(tmp + 0.5);  // TODO: ceiling?
		// TODO: TRY THIS xs_int[i] = ceiling(tmp);
		xs_frac[i] = tmp - xs_int[i];
	}
	return;
}

inline void update_r(FVector& r, FVector& xs_int)
{
	Float scalar;
	size_t n = r.size();
	int64_t shifted = ((int64_t)1 << shift);
	_field.init(scalar, (double)shifted);
	FVector y(n);

	//update r = r * 2^shift - Mat*xs_int
	_VDF.mulin(r, scalar);
	_numsolver.apply(y, xs_int);
	_VDF.subin(r, y);

	return;
} // update_r

template <class IMatrix>
inline void update_r_exact(IVector& r_exact, FVector& r, FVector& xs_int, IMatrix &IM){
	size_t n = r.size();

	IVector x_i(n), y_i(n);

	typename Ring::Element scalar = ((int64_t)1 << shift);

	// update r = r * 2^shift - Mat*xs_int
	//  r *= 2 ^shift
	_VDR.mulin(r_exact, scalar);

	// determine if exact apply is needed
	double vnorm = zw_dOOnorm(&*xs_int.begin(), (int)n, 1);

	int64_t th = ((int64_t)1 << 52);  // double mantissa
	Float thresh;
	_field.init(thresh, (double)th);

	debugneol("vnorm " << vnorm);

	//  r -= Mat * xs_int
	if(_field.mulin(vnorm, mnorm) < thresh){
		debugneol("Numeric ");
		FVector y(n);
		_numsolver.apply(y, xs_int);
		for(size_t i = 0; i < n; ++i)
			_ring.init(y_i[i], y[i]);
	}
	else{
		SHIFT_BOUND--; // to less this possibility
		debugneol("Exact ");
		for(size_t i = 0; i < n; ++i)
			_ring.init(x_i[i], xs_int[i]);
		IM.apply(y_i, x_i);
	}

	_VDR.subin(r_exact, y_i);

	//  convert exactly computed residue back to double (for next solve)
	typename FVector::iterator rp = r.begin();
	typename IVector::iterator rep = r_exact.begin();
	for(; rp!= r.end(); ++rp, ++rep)
		_field.init(*rp, *rep);

	return;
} // update_r_exact

// no longer called...
inline int HadamardBound(integer& B, FMatrix& DM)
{
	size_t n = DM.rowdim();
	zw_hbound (B, n, n, &*DM.Begin()); // compute the Hadamard bound
	B = B * B;
	double mnorm_loc = zw_dOOnorm(&*DM.Begin(), n, n);

	// [don't know what this comment is about] should be a check for 2 * mnorm + zw_dmax (n, b, 1);
	// TODO what is "b"? from copied code it is the RHS array of doubles
	// zw_max just seems to get abs(max value of b)
	// next line false, just to compile
	double *b;
	B *= 2 * mnorm_loc + zw_dmax (n, b, 1); // [don't know what this factor is about]
	B <<= 1; // [extra factor of 2 for some reason... ]
	return B;
}

//update num, *num <- *num * 2^shift + d
inline IVector& update_num (IVector& num, const FVector& d)
{
	size_t n = d.size();
	IVector d_i(n);
	for (size_t i = 0; i < n; ++i) {
		_ring.init(d_i[i], d[i]);
	}
	Int scalar; _ring.init(scalar, 1UL << shift);
	//  TODO - analyze GMP shifting capability
	_VDR.mulin(num, scalar);
	_VDR.addin(num, d_i);
	return num;
}

//update r = r * shift - M d
inline static int update_r_ll (double* r, int n, const double* M, const double* d, int shift)
{
	long long int tmp;
	double* p1;
	const double* p2;
	const double* pd;
	for (p1 = r, p2 = M; p1 != r + n; ++ p1) {
		tmp = (long long int) *p1;
		tmp <<= shift;
		for (pd = d; pd != d + n; ++ pd, ++ p2) {
			tmp -= (long long int)*pd * (long long int) *p2;
		}
		*p1 = tmp;
	}
	return 0;
}

inline static size_t nextPower2(size_t n)
{
	size_t p = 1;
	while(p < n) p <<= 1;
	return p;
}

inline static double highAbs(FMatrix M)
{
	double max = 0;
	typename FMatrix::Iterator ri = M.Begin();
	for(; ri != M.End(); ++ri){
		double tmp = fabs(*ri);
		if(max < tmp) max = tmp;
	}
	return max;
}

inline static double zw_dOOnorm(const double* M, int m, int n)
{
	double norm = 0;
	double old = 0;
	const double* p;
	for (p = M; p != M + (m * n); ) {
		old = norm;
		norm = cblas_dasum (n, p ,1);
		if (norm < old) norm = old;
		p += n;
	}
	return norm;
}

inline static double zw_dmax (const int N, const double* a, const int inc)
{
	return fabs(a[cblas_idamax (N, a, inc)]);
}

inline static int zw_hbound (integer& b, int m, int n, const double* M)
{
	double norm = 0;
	const  double* p;
	integer tmp;
	b = 1;
	for (p = M; p != M + (m * n); ) {
		norm = cblas_dnrm2 (n, p ,1);
		tmp =  norm;
		integer::mulin (b, tmp);
		p += n;
	}

	return 0;
}

/*  out:  vector to print
 *  tag:  prepend this to output
 *  bound:  upper bound on printed value (0 - ignore bound) */
template<class Vec>
std::ostream& writeVec(Vec& out, const char *tag="", integer bound=/*40000000*/0,
		       size_t numEntries=5, std::ostream &os=std::cerr)
{
	os << tag << ": [";

	size_t n = (out.size() < numEntries ? out.size() : numEntries);
	for( size_t i=0; i<n; ++i){
		if(bound && ((integer)(out[i]) > bound || (integer)(out[i]) < -bound)){
			os << " entry over bound ]"  << std::endl;
			return os;
		}
		else
			os << " " << out[i];
	}
	if(out.size() > numEntries)
		os << " ...";
	os << " ]" << std::endl;
	return os;
}

template<class Vec>
void writeVecFile(Vec& out, const char* file)
{
	std::ofstream os;
	os.open(file, std::ios::out);

	typename Vec::const_iterator vi = out.begin();

	for( ; vi != out.end(); ++vi){
		os << *vi << std::endl;
	}
	os.close();
}

// to write diagnostic info to files for further testing
template <class Matrix>
void dumpData(const Matrix &M, const IVector &b, IVector &numx, integer &denx, integer &denBound)
{
#ifdef SPITOUT
	std::ofstream matout;
	matout.open("debug.mat", std::ios::out);
	M.write(matout);
	matout.close();
	writeVecFile(b, "debug.rhs");
	writeVecFile(numx, "debug.num");
	std::ofstream dout;
	dout.open("debug.den", std::ios::out);
	dout << denx << std::endl;
	dout << denBound << std::endl;
	dout << numx.size() << std::endl;
	dout.close();
#endif
}

// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,:0,t0,+0,=s
// Local Variables:
// mode: C++
// tab-width: 8
// indent-tabs-mode: nil
// c-basic-offset: 8
// End: