1 // -*- C++ -*- 2 //============================================================================================== 3 // 4 // This file is part of LiDIA --- a library for computational number theory 5 // 6 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved. 7 // 8 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/ 9 // 10 //---------------------------------------------------------------------------------------------- 11 // 12 // $Id$ 13 // 14 // Author : Victor Shoup (VS), Thomas Pfahler (TPf) 15 // Changes : See CVS log 16 // 17 //============================================================================================== 18 19 20 #ifndef LIDIA_FP_POLYNOMIAL_UTIL_H_GUARD_ 21 #define LIDIA_FP_POLYNOMIAL_UTIL_H_GUARD_ 22 23 24 #ifndef HEADBANGER 25 26 27 #ifndef LIDIA_FP_POLYNOMIAL_H_GUARD_ 28 # include "LiDIA/Fp_polynomial.h" 29 #endif 30 #ifndef LIDIA_TIMER_H_GUARD_ 31 # include "LiDIA/timer.h" 32 #endif 33 34 35 36 #ifdef LIDIA_NAMESPACE 37 namespace LiDIA { 38 # define IN_NAMESPACE_LIDIA 39 #endif 40 41 42 43 //*************************************************************** 44 // 45 // tools 46 // 47 //*************************************************************** 48 49 50 51 52 class my_timer 53 { 54 timer t; 55 char *msg; 56 public: 57 my_timer(); 58 ~my_timer(); 59 60 void start(const char* message); 61 void stop(); 62 }; 63 64 65 66 67 68 /**************************************************************** 69 70 class poly_matrix 71 72 ****************************************************************/ 73 74 class poly_matrix 75 { 76 private: 77 Fp_polynomial elts[2][2]; 78 79 poly_matrix(); 80 poly_matrix(const poly_matrix&); // disable 81 82 83 void iter_half_gcd(Fp_polynomial& U, Fp_polynomial& V, lidia_size_t d_red); 84 85 86 void multiply_and_kill(poly_matrix& B, poly_matrix& C); 87 // (*this) = B*C, B and C are destroyed 88 public: 89 90 poly_matrix(const bigint &p); 91 ~poly_matrix()92 ~poly_matrix() { } 93 94 void kill(); 95 96 poly_matrix & operator = (const poly_matrix&); 97 operator()98 Fp_polynomial& operator() (lidia_size_t i, lidia_size_t j) 99 { 100 //NO CHECK 101 return elts[i][j]; 102 } 103 operator()104 const Fp_polynomial& operator() (lidia_size_t i, lidia_size_t j) const 105 { 106 //NO CHECK 107 return elts[i][j]; 108 } 109 110 void half_gcd(const Fp_polynomial& U, const Fp_polynomial& V, lidia_size_t d_red); 111 // deg(U) > deg(V), 1 <= d_red <= deg(U)+1. 112 // 113 // This computes a 2 x 2 polynomial matrix M_out such that 114 // M_out * (U, V)^T = (U', V')^T, 115 // where U', V' are consecutive polynomials in the Euclidean remainder 116 // sequence of U, V, and V' is the polynomial of highest degree 117 // satisfying deg(V') <= deg(U) - d_red. 118 119 void xhalf_gcd(Fp_polynomial& U, Fp_polynomial& V, lidia_size_t d_red); 120 // same as above, except that U is replaced by U', and V by V' 121 122 void multiply(Fp_polynomial& U, Fp_polynomial& V) const; 123 // (U, V)^T = (*this) * (U, V)^T 124 125 }; 126 127 128 129 //**************************************************************** 130 // 131 // class poly_argument 132 // 133 //**************************************************************** 134 135 // The routine poly_argument::build (see below) which is implicitly called 136 // by the various compose and update_map routines builds a table 137 // of polynomials. If poly_arg_bound > 0, then only up to 138 // about poly_arg_bound bigints mod p are allocated. Setting this too 139 // low may lead to slower execution. 140 // If poly_arg_bound <= 0, then it is ignored, and space is allocated 141 // so as to maximize speed. 142 // Initially, poly_arg_bound = 0. 143 144 // If a single h is going to be used with many g's 145 // then you should build a poly_argument for h, 146 // and then use the compose routine below. 147 // poly_argument::build computes and stores h, h^2, ..., h^m mod f. 148 // After this pre-computation, composing a polynomial of degree 149 // roughly n with h takes n/m multiplies mod f, plus n^2 150 // scalar multiplies. 151 // Thus, increasing m increases the space requirement and the pre-computation 152 // time, but reduces the composition time. 153 // If poly_arg_bound > 0, a table of size less than m may be built. 154 155 156 //template < class T > class factorization; 157 158 class poly_argument 159 { 160 base_vector< Fp_polynomial > H; 161 162 //PolyargBound can only be changed if no poly_arguments 163 //have been declared yet 164 static lidia_size_t poly_arg_bound; 165 static bool change_enable; 166 167 static 168 void compose_InnerProd(Fp_polynomial& x, const Fp_polynomial& v, 169 lidia_size_t low, lidia_size_t high, 170 const base_vector< Fp_polynomial > & H, lidia_size_t n, 171 bigint *t); 172 173 public: poly_argument()174 poly_argument() 175 { 176 change_enable = false; 177 debug_handler("poly_argument", "poly_argument (void)"); 178 } 179 ~poly_argument()180 ~poly_argument() 181 { 182 debug_handler("poly_argument", "destructor"); 183 } 184 185 static void set_poly_arg_bound(lidia_size_t n); get_poly_arg_bound()186 static lidia_size_t get_poly_arg_bound() 187 { 188 return poly_arg_bound; 189 } 190 191 void build(const Fp_polynomial& h, const Fp_poly_modulus& F, lidia_size_t m); 192 //m must be > 0, otherwise an error is raised 193 194 195 void compose(Fp_polynomial& x, const Fp_polynomial& g, 196 const Fp_poly_modulus& F) const; 197 //called by the various compose-routines below 198 199 200 friend void project_powers(base_vector< bigint > & x, 201 const base_vector< bigint > & a, lidia_size_t k, 202 const poly_argument& H, const Fp_poly_modulus& F); 203 }; 204 205 206 207 //**************************************************************** 208 // 209 // class int_factor 210 // fac_vec " = " vector< int_factor > 211 // 212 //**************************************************************** 213 214 class fac_vec; 215 216 class int_factor 217 { 218 public: 219 lidia_size_t q; // base 220 int a; // exponent 221 int b; 222 double len; 223 224 private: int_factor()225 int_factor() {}; // allow creation only in class fac_vec ~int_factor()226 ~int_factor() {}; 227 int_factor & operator = (const int_factor c) 228 { 229 q = c.q; 230 a = c.a; 231 len = c.len; 232 b = c.b; 233 return *this; 234 } 235 236 friend class fac_vec; 237 }; 238 239 240 241 class fac_vec 242 { 243 int_factor *fvec; 244 int num_factors; 245 246 fac_vec(); // disable 247 fac_vec(const fac_vec &); 248 249 void compute_factorization(lidia_size_t n); 250 void sort(); 251 252 public: 253 fac_vec(lidia_size_t n); 254 ~fac_vec()255 ~fac_vec() 256 { 257 debug_handler("fac_vec", "destructor"); 258 delete[] fvec; 259 } 260 261 void clear(int lo, int hi); 262 //fvec[i].b = 0 for i = lo, .., hi 263 264 const int_factor& operator[] (int n) const 265 { 266 debug_handler("fac_vec", "operator [] (int) const"); 267 return fvec[n]; 268 } 269 270 int_factor& operator[] (int n) 271 { 272 debug_handler("fac_vec", "operator [] (int)"); 273 return fvec[n]; 274 } 275 number_of_factors()276 int number_of_factors() const 277 { 278 debug_handler("fac_vec", "number_of_factors (void)"); 279 return num_factors; 280 } 281 282 lidia_size_t prod(int lo, int hi) const; 283 //returns product of q^a for i = lo, .., hi 284 }; 285 286 287 288 //**************************************************************** 289 // 290 // fractions 291 // 292 //**************************************************************** 293 294 // ************ addition, test for equality ************ 295 296 void add_frac(Fp_polynomial& x, Fp_polynomial& y, const Fp_polynomial& a, 297 const Fp_polynomial& b, const Fp_polynomial& c, const Fp_polynomial& d, 298 const Fp_poly_modulus& F); 299 // computes (x, y) = (a*d + b*c mod f, b*d mod f) 300 // ALIAS RESTRICTION: inputs may not alias outputs 301 302 void subtract_frac(Fp_polynomial& x, Fp_polynomial& y, const Fp_polynomial& a, 303 const Fp_polynomial& b, const Fp_polynomial& c, const Fp_polynomial& d, 304 const Fp_poly_modulus& F); 305 // computes (x, y) = (a*d - b*c mod f, b*d mod f) 306 // ALIAS RESTRICTION: inputs may not alias outputs 307 308 309 void plain_add_frac(Fp_polynomial& x, Fp_polynomial& y, const Fp_polynomial& a, 310 const Fp_polynomial& b, const Fp_polynomial& c, const Fp_polynomial& d, 311 const Fp_poly_modulus& F); 312 313 void plain_subtract_frac(Fp_polynomial& x, Fp_polynomial& y, 314 const Fp_polynomial& a, const Fp_polynomial& b, const Fp_polynomial& c, 315 const Fp_polynomial& d, const Fp_poly_modulus& F); 316 317 // same as above, but slower 318 319 320 321 bool eq_frac(const Fp_polynomial& a, const Fp_polynomial& b, 322 const Fp_polynomial& c, const Fp_polynomial& d, const Fp_poly_modulus& F); 323 // if (a*d - b*c = 0 mod f) then 1 else 0 324 325 326 bool plain_eq_frac(const Fp_polynomial& a, const Fp_polynomial& b, 327 const Fp_polynomial& c, const Fp_polynomial& d, const Fp_poly_modulus& F); 328 // same as above, but slower 329 330 331 332 333 // ************ class eq_frac_info ************ 334 335 // next come data structures and routines for testing if a fixed 336 // a/b is equal to one of several c/d mod f, using a fast probabilistic test 337 // With a very small probability, two unequal fractions may be 338 // claimed to be equal. 339 340 class eq_frac_info 341 { 342 base_vector< bigint > L1, L2; 343 bigint modulus; 344 345 public: 346 347 void build(const Fp_polynomial& a, const Fp_polynomial& b, 348 const Fp_poly_modulus& F); 349 350 bool eq_frac(const Fp_polynomial& c, const Fp_polynomial& d) const; 351 }; 352 353 354 355 #endif // HEADBANGER 356 357 358 359 #ifdef LIDIA_NAMESPACE 360 } // end of namespace LiDIA 361 # undef IN_NAMESPACE_LIDIA 362 #endif 363 364 365 366 #endif // LIDIA_FP_POLYNOMIAL_UTIL_H_GUARD_ 367