/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
|