This file is indexed.

/usr/include/eclib/sieve_search.h is in libec-dev 20160101-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
// sieve_search.h: declarations of classes point_processor and qsieve
//////////////////////////////////////////////////////////////////////////
//
// Copyright 1990-2012 John Cremona
// 
// This file is part of the eclib package.
// 
// eclib is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at your
// option) any later version.
// 
// eclib 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 General Public License
// for more details.
// 
// You should have received a copy of the GNU General Public License
// along with eclib; if not, write to the Free Software Foundation,
// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
// 
//////////////////////////////////////////////////////////////////////////
 
// File        : sieve_search.h
// Author      : Sophie Labour
//		 Adaption of M. Stoll's code
// Last change : SL, Apr 21 1999, first definitions
// Last change : JC, Dec 2 1999, adapted for direct use in mwrank et al.
// Last change : JC, Aug 19 2002, adapted for stdc++ library

#ifndef SIEVE_SEARCH_H
#define SIEVE_SEARCH_H

class point_processor { //An abstract class to be used as interface in qsieve
 public:
  virtual int process(const bigint& x,const bigint& y,const bigint& z)=0;
  // where x/z is the point found and y the square root of f(a,b);
  // #!# different if the degree is even or odd, poly monic or not;
  // true homogeneous coordinates are [a*b^n-1:y:b^n] where n=degree/2 if 
  // degree even or degree odd & polyn monic (b is a square so b^n is 
  // (sqrt(b))^degree), n=(degree+1)/2 otherwise 
  // if no_check option, process(x,0,z) is called
  // the int return value is to signal the search to stop (if set to 1)
  virtual ~point_processor() {;}
};

class point_printer : public point_processor {
 public:
  point_printer() {};
  ~point_printer() {};
  int process(const bigint& xx, const bigint& yy, const bigint& zz) 
    {cout<<"x= "<<xx<<" y= "<<yy<<" z= "<<zz<<" is a point."<<endl; return 0;}
};

class point_counter : public point_processor {
  int tally;
 public:
  point_counter() {tally=0;};
  ~point_counter() {};
  int process(const bigint& xx, const bigint& yy, const bigint& zz) 
    {tally++; return 0;}
  int get_tally() {return tally;}
};

// *****************************************************
//          class qsieve
// *****************************************************

#define QS_VERSION "This is an adaptation of ratpoints-1.4 by M. Stoll"

#define QS_MAX_DEGREE 10  // max. degree of f(x) 
#define QS_DEFAULT_SIZE 10     // Default value for the array_size option 
#define QS_DEFAULT_TECH 1      // Default value for the tech option 

#define QS_LONG_LENGTH (8 * sizeof(unsigned long))
#define QS_LONG_SHIFT ((QS_LONG_LENGTH == 16) ? 4 : \
                    (QS_LONG_LENGTH == 32) ? 5 : \
		    (QS_LONG_LENGTH == 64) ? 6 : 0)
#define QS_LONG_MASK (~(-1L<<QS_LONG_SHIFT))
#define QS_HALF_MASK ((bit_array)(~((unsigned long)(-1L) / 3)))

     //#define PERCENT 0.6
     //#define SIEVE_PRIMES 19
#define QS_NUM_PRIMES 53 //42 //53
#define QS_MAX_PRIME 251 //191 //251
     //#define QS_MAX_PRIME_EVEN 256
#define QS_DEFAULT_HEIGHT 10

typedef unsigned long bit_array;
#define bit_zero ((bit_array)0)
#define all_ones ((bit_array)(-1L))
typedef struct {long p; double r; long n;} entry;
typedef struct {double low; double up;} interval;
typedef struct {long p; bit_array *ptr;} sieve_spec;


class qsieve {
  bigint c[QS_MAX_DEGREE+1];
  point_processor* curve;
  int degree;
  int verbose;
  
  //"global variables"
  char init_made;
  static long prime[QS_NUM_PRIMES];
  bit_array* bits;
  char** squares;
 //squares[pn][x] = 1 if x is square mod prime[pn], 0 if not 

  char** is_f_square;//[QS_NUM_PRIMES][QS_MAX_PRIME];
  // is_f_square[pn][x] = 1 if f(x) is square mod prime[pn], 0 if not 
  long* pnn;//[QS_NUM_PRIMES];
  // holds the number of the primes used for the sieve
  long* kpa;//[QS_NUM_PRIMES];    
  // stores suitable multiples of sieving primes, depending on tech

  // ratios2[d] gives the speed ratio of testing/sieving for degree d.
  // The following values are reasonable for polynomials with coefficients
  // the size as in the examples, on a Pentium. 
  static double ratio1_def;
  static double ratios2[QS_MAX_DEGREE+1];
  static double ratio1;  /* This will be set if the -r option is given */
  static double ratio2;  /* This will be set if the -R option is given */

  bit_array*** sieve;//[QS_NUM_PRIMES][QS_MAX_PRIME_EVEN][QS_MAX_PRIME_EVEN];
  //bit k of sieve[pn][b1][a1] = 0 iff a/b is excluded mod p = prime[pnn[pn]], when  a = a1*QS_LONG_LENGTH + k mod p and b = b1 mod p. 
  bit_array*** sieve2;//[QS_NUM_PRIMES][QS_MAX_PRIME_EVEN][QS_MAX_PRIME_EVEN];
  // bit k of sieve2[pn][b1][a1] = 0 iff (2*a + 1)/b is excluded mod p = prime[pnn[pn]], when  a = a1*QS_LONG_LENGTH + k mod p and b = b1 mod p and b is even.

  sieve_spec sieves[QS_NUM_PRIMES];

  long coeffs_mod_p[QS_NUM_PRIMES][QS_MAX_DEGREE+1];
  bigint bc[QS_MAX_DEGREE+1]; //The coefficients of f, multiplied by powers of b
  bigint fff, tmp, tmp2;   //Some multi-precision integer variables 
  long sieve_primes1;  //The number of primes used for the first sieving stage
  long sieve_primes2; //The number of primes used for the both sieving stages
  bit_array *survivors;
  entry prec[QS_NUM_PRIMES];  //This array is used for sorting in order to determine the `best' sieving primes.

  long height;         //The height bound
  long w_height;       //The height bound divided by the word length in bits 
  int points_at_infty; //A flag saying if we should look for points at infty
  long b_low, b_high;  //Lower and upper bounds for the denominator
  int halt_flag;       // gets set when searching can stop
  int use_opt;         //A flag saying if to use optimised sieving for even denominators
  int no_check;        //A flag that says whether to omit the final check
  long array_size;     //The size of the survivors array (in longs: converted from kbytes to longs at beginning of search)
  long num_inter;      //The number of intervals in the search region 
  interval domain[QS_MAX_DEGREE+1];  //This contains the intervals representing the search region
  static long tech;
  int check_denom; //A flag saying if we should check if the denom is divisible by a `forbidden divisor'
  long forbidden[QS_NUM_PRIMES+2]; //The forbidden divisors, a zero-terminated array.
  int use_squares; //A flag that is set when we can use only squares as denominators (odd degree & polynomial monic)
  int odd_nums;    //A flag that says that we need only consider odd numerators
  int compute_bc;

  long num_surv1;  //Used to count the survivors of the first stage
  long num_surv2;  //Used to count the survivors of the second stage
  

  long w_floor(long a,long b)
    { if (a < 0) return -(1 + (-a-1)/b); else return a/b;}; 
  long w_ceil(long a,long b)
    { if (a <= 0) return -((-a)/b); else return 1+(a-1)/b;};

  void init_data();
  void init_all();//init bits and squares
  void init_f();//init is_f_square and sieve

  long sift(long b);
  long sift0(long b, long w_low, long w_high, int use_odd_nums);
  void check_point(bit_array nums, long b, long i, long* total, int use_odd_nums);
  //  void a_search(const long& amin, const long& amax);
  //  void a_simple_search(const long& amin, const long& amax);
  
  void dealloc_sieves();

 public:
  qsieve() {;}
  qsieve(point_processor* acurve, int deg, vector<bigint> coeff, int verb=0);
  qsieve(point_processor* acurve, int deg, vector<bigint> coeff, bigfloat h_limx, 
	 int verb=0); 
//qsieve(point_processor* acurve, int deg, vector<bigint> coeff, double h_lim, 
//       double up, double low, int verb=1);
  ~qsieve();

  void set_intervals(vector<double> interv,int nb_bnd,int start_low,int pos_x_only=0);
  //the search will be made on the intervals thus defined:
  //nb_bnd bounds are given in interv, the first one being a lower bound 
  //for the first interval if start_low=1, and the upper bound for the first
  // interval (the lower bound being -height) otherwise.
  //we must have interv[i]<interv[i+1], otherwise, behaviour undefined
  //there can be at most QS_MAX_DEGREE+1 intervals.
  //pos_x_only is 1 if only positive solutions are wanted
  //IMPORTANT: If height is to be different than QS_DEFAULT_HEIGHT and
  //constructor was not called with height argument, set_height() must 
  //be called before set_intervals()

  //All the following set_ functions (and the others, up to set_b_high) 
  //should be called (if called at all) BEFORE set_intervals, and all of 
  //those, before calling a search.

  //Sets the (naive logarithmic) height to which the search is performed 
  //for the x coordinate
  void set_height(double h_lim) 
    {height=(long)(floor(exp(h_lim)));w_height=(height-1)/QS_LONG_LENGTH + 1;};

  void set_sieve_primes1(long sp1);
  //sets the number of primes used for the first stage of sieving to sp1, 
  //unless sieve_primes2 has already been set to a smaller value.  
  //Maximum is QS_NUM_PRIMES.

  void set_sieve_primes2(long sp2);
  //sets the total number of primes used for sieving to sp2, unless 
  //sieve_primes1 has already been set to a greater value.  
  //Maximum is QS_NUM_PRIMES.

  //`ratio1' is the ratio of running time of the second versus the first 
  //stage of sieving (per bit). This is used to find the optimal number 
  //of sieving primes for the first stage automatically.
  void set_ratio1(double r1)
    { if (r1<0.0) 
      cout<<"qsieve::set_ratio1():ratio must be >0"<<endl; 
    else ratio1=r1; };
  //`ratio2' is the ratio of running time needed for checking if x-coordinates
  //give rise to points versus one step of the second sieving stage. 
  //This is used to find the optimal number of sieving primes for the second 
  //stage automatically
  void set_ratio2(double r2)
    { if (r2<0.0) 
      cout<<"qsieve::set_ratiow():ratio must be >0"<<endl; 
    else ratio2=r2; };
  //`size' is the size (in kilobytes) used for the array of bits in which 
  //the sieving is done. This affects the amount of memory used by the 
  //program. A smaller value might give better performance if the array can 
  //be held in the cache in its entirety. So it is to be expected that the 
  //effect of this parameter depends heavily on the cache size and speed. 
  //Default is DEFAULT_SIZE = 10. 
  void set_array_size(long as)
    { if (as<0) 
      cout<<"qsieve::set_array_size(): size must be >0"<<endl; 
    else array_size=as;};//size of survivors array in kbytes
  void be_quiet() {verbose=0;};
  //to switch off the check if a survivng candidate is really a point. 
  //This is useful when you use another program and this program does the 
  //test for itself
  void set_no_check() {no_check=1;};
  //If we do not want it
  void no_points_at_infinity() {points_at_infty=0;};
  //lower bound for the denominator (default:1)
  void set_b_low(long b) {if (b>0) b_low=b;};
  //upper bound for the denominator (default:height)
  void set_b_high(long b) {if (b>0) b_high=b;};


  //all perform sieve-assisted fast search according to the info provided:
  //h_lim(logarithmic, default:5), b_low,b_high(bounds z, optional)
  //intervals for x/z, optional, must be set AFTER height
  //the sieve uses height, so it is set every time search is called
  long search();
  long search(double h_lim);
  
};

/* This is a comparison function needed for sorting in order to determine
   the `best' primes for sieving. */
int compare_entries(const void *a, const void *b);

#endif