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 : Michael Jacobson (MJ)
15 // Changes : See CVS log
16 //
17 //==============================================================================================
18
19
20 #ifndef LIDIA_PARTIAL_RELATION_H_GUARD_
21 #define LIDIA_PARTIAL_RELATION_H_GUARD_
22
23
24 #ifndef LIDIA_BIGINT_H_GUARD_
25 # include "LiDIA/bigint.h"
26 #endif
27 #ifndef LIDIA_MATH_VECTOR_H_GUARD_
28 # include "LiDIA/math_vector.h"
29 #endif
30 #ifndef LIDIA_SORT_VECTOR_H_GUARD_
31 # include "LiDIA/sort_vector.h"
32 #endif
33 #ifndef LIDIA_MATRIX_H_GUARD_
34 # include "LiDIA/matrix.h"
35 #endif
36 #ifndef LIDIA_QUADRATIC_NUMBER_STANDARD_H_GUARD_
37 # include "LiDIA/quadratic_number_standard.h"
38 #endif
39 #include <cstring>
40 #include <cstdlib>
41
42
43
44 #ifdef LIDIA_NAMESPACE
45 namespace LiDIA {
46 # define IN_NAMESPACE_LIDIA
47 #endif
48
49
50
51 //
52 // Class: partial_relation
53 //
54 // This class represents a partial relation found from a sieving algorithm.
55 // It consists of an integer representing the value of the large prime,
56 // an index which distinguishes an instance of this class from other
57 // instances, and a string representing the entire relation.
58 //
59
60 #define REL_STR_LEN 1000
61 #define LINE_LEN 1000
62
63 class partial_relation
64 {
65 protected:
66
67 int lp; // value of the large prime
68 int ex; // exponent of the large prime (1 = negative)
69 bigint H; // index
70 char str[REL_STR_LEN]; // string representing the relation
71 quadratic_number_standard g; // the principal ideal generator
72
73
74 public:
75
partial_relation()76 partial_relation()
77 {
78 lp = ex = 0;
79 H.assign_zero();
80 memset(str, '\0', REL_STR_LEN);
81 g.assign_one();
82 }
83
~partial_relation()84 ~partial_relation()
85 {}
86
get_large()87 int get_large() const
88 {
89 return lp;
90 }
91
get_exp()92 int get_exp() const
93 {
94 return ex;
95 }
96
get_index()97 bigint get_index() const
98 {
99 return H;
100 }
101
get_min()102 quadratic_number_standard get_min() const
103 {
104 return g;
105 }
106
strcopy(char * st)107 void strcopy(char *st)
108 {
109 strcpy(st, str);
110 }
111
assign_zero()112 void assign_zero()
113 {
114 lp = ex = 0;
115 H.assign_zero();
116 memset(str, '\0', REL_STR_LEN);
117 g.assign_one();
118 }
119
assign(bigint & n_H,char * n_str,quadratic_number_standard & q)120 void assign(bigint & n_H, char *n_str, quadratic_number_standard & q)
121 {
122 lp = 1;
123 ex = 0;
124 H = n_H;
125 strcpy(str, n_str);
126 g.assign(q);
127 }
128
assign(int n_lp,int n_ex,bigint & n_H,char * n_str,quadratic_number_standard & q)129 void assign(int n_lp, int n_ex, bigint & n_H, char *n_str, quadratic_number_standard & q)
130 {
131 lp = n_lp;
132 ex = n_ex;
133 H = n_H;
134 strcpy(str, n_str);
135 g.assign(q);
136 }
137
138 partial_relation operator = (const partial_relation & B)
139 {
140 lp = B.lp;
141 ex = B.ex;
142 H.assign(B.H);
143 strcpy(str, B.str);
144 g.assign(B.g);
145 return *this;
146 }
147
148
149
150 friend int compare (const partial_relation & P1, const partial_relation & P2);
151
is_zero()152 bool is_zero()
153 {
154 return (lp == 0);
155 }
156
is_full_relation()157 bool is_full_relation() const
158 {
159 return (lp == 1);
160 }
161
162
163
164 void add_to_matrix(matrix< int > & M,
165 base_vector< quadratic_number_standard > &vec_qn,
166 lidia_size_t j);
167
168
169
170 friend void swap(partial_relation & A, partial_relation & B);
171
172 friend std::istream & operator >> (std::istream & in, partial_relation & A);
173
174 friend std::ostream & operator << (std::ostream & out, const partial_relation & A);
175
176 friend bool make_full(partial_relation & C,
177 partial_relation & A,
178 partial_relation & B,
179 bigint & N,
180 int size_exp);
181
182 };
183
184
185
186 inline void
add_to_matrix(matrix<int> & M,base_vector<quadratic_number_standard> & vec_qn,lidia_size_t j)187 partial_relation::add_to_matrix (matrix< int > & M,
188 base_vector< quadratic_number_standard > &vec_qn,
189 lidia_size_t j)
190 {
191
192 lidia_size_t i;
193 int e;
194 char helpstr[LINE_LEN], *pp, *p;
195
196 // write the relation to the last column of M
197 strcpy(helpstr, str);
198 pp = helpstr;
199 p = strtok(pp, " \n");
200 while (p != NULL) {
201 e = atoi(p);
202 if (!e)
203 break;
204 p = strtok(NULL, " \n");
205 i = static_cast<lidia_size_t>(std::atoi(p));
206 M.sto(i, j, e);
207 p = strtok(NULL, " \n");
208 }
209
210 // write minimum to jth entry of vec
211 vec_qn[j].assign(g);
212 }
213
214
215
216 inline std::istream &
217 operator >> (std::istream & in, partial_relation & A)
218 {
219 char line[LINE_LEN], *p;
220
221
222 in.getline(line, LINE_LEN);
223 if (!in.eof()) {
224 p = strchr(line, '@');
225 if (p) {
226 // reading a parital relation
227 p = strtok(line, "@");
228 A.lp = atoi(p);
229 p = strtok(NULL, "@");
230 A.ex = atoi(p);
231 p = strtok(NULL, ":");
232 }
233 else {
234 // reading a full relation
235 A.lp = 1;
236 A.ex = 0;
237 p = strtok(line, ":");
238 }
239 string_to_bigint(p, A.H);
240 p = strtok(NULL, ":");
241 ++p;
242 strcpy(A.str, p);
243 p = strtok(NULL, ":");
244 ++p;
245 string_to_quadratic_number_standard(p, A.g);
246 }
247 else
248 A.assign_zero();
249
250 return in;
251 }
252
253
254
compare(const partial_relation & P1,const partial_relation & P2)255 inline int compare (const partial_relation & P1, const partial_relation & P2)
256 {
257 if (&P1 == &P2) {
258 return 0;
259 }
260 if (P1.lp == P2.lp) {
261 return 0;
262 }
263 return ((P1.lp < P2.lp) ? -1 : 1);
264 }
265
266
267
268 inline bool operator == (const partial_relation & P1, const partial_relation & P2)
269 {
270 return (compare(P1, P2) == 0);
271 }
272
273
274
275 inline bool operator < (const partial_relation & P1, const partial_relation & P2)
276 {
277 return (compare(P1, P2) < 0);
278 }
279
280
281
282 inline bool operator <= (const partial_relation & P1, const partial_relation & P2)
283 {
284 return (compare(P1, P2) <= 0);
285 }
286
287
288
289 inline bool operator > (const partial_relation & P1, const partial_relation & P2)
290 {
291 return (compare(P1, P2) > 0);
292 }
293
294
295
296 inline bool operator >= (const partial_relation & P1, const partial_relation & P2)
297 {
298 return (compare(P1, P2) >= 0);
299 }
300
301
302
303 inline std::ostream &
304 operator << (std::ostream & out, const partial_relation & A)
305 {
306 if (A.is_full_relation())
307 out << A.H << " :" << A.str << ": " << A.g << "\n";
308 else
309 out << A.lp << " @ " << A.ex << " @ " << A.H << " : " << A.str << ": " << A.g << "\n";
310 return out;
311 }
312
313
314
315 inline void
swap(partial_relation & A,partial_relation & B)316 swap(partial_relation & A, partial_relation & B)
317 {
318 partial_relation C;
319
320
321 C = A;
322 A = B;
323 B = C;
324 }
325
326
327
328 inline bool
make_full(partial_relation & C,partial_relation & A,partial_relation & B,bigint & N,int size_exp)329 make_full(partial_relation & C, partial_relation & A, partial_relation & B, bigint & N, int size_exp)
330 {
331 if (A.lp == B.lp) {
332 char helpstr[REL_STR_LEN], *p, *pp;
333 bigint h_inv, temp;
334 int e;
335 lidia_size_t i;
336 bool neg, zero;
337 short *exp = NULL;
338
339 exp = new short[size_exp];
340 memset(exp, '\0', size_exp*sizeof(short));
341
342 neg = (A.ex == B.ex);
343
344 // do exponent arithmetic
345 strcpy(helpstr, A.str);
346 pp = helpstr;
347 p = strtok(pp, " \n");
348 while (p != NULL) {
349 e = atoi(p);
350 if (!e)
351 break;
352 p = strtok(NULL, " \n");
353 exp[atoi(p)] += e;
354 p = strtok(NULL, " \n");
355 }
356
357 zero = true;
358 strcpy(helpstr, B.str);
359 pp = helpstr;
360 p = strtok(pp, " \n");
361 while (p != NULL) {
362 e = atoi(p);
363 if (!e)
364 break;
365 p = strtok(NULL, " \n");
366 if (neg)
367 exp[atoi(p)] -= e;
368 else
369 exp[atoi(p)] += e;
370 if (exp[atoi(p)])
371 zero = false;
372 p = strtok(NULL, " \n");
373 }
374
375 // exit if relation is zero
376 if (zero) {
377 C.assign_zero();
378 delete [] exp;
379 return true;
380 }
381
382 temp = xgcd_left(h_inv, B.H, N);
383 C.H = (h_inv * A.H) % N;
384
385 C.ex = 0;
386 C.lp = 1;
387 memset(C.str, '\0', REL_STR_LEN);
388 for (i = 0; i < size_exp; i++) {
389 if (exp[i])
390 sprintf(C.str, "%s %d %d", C.str, exp[i], i);
391 }
392 sprintf(C.str, "%s 0", C.str);
393
394 if (N.is_gt_zero()) {
395 if (neg)
396 divide(C.g, A.g, B.g);
397 else
398 multiply(C.g, A.g, B.g);
399 }
400 else
401 C.g.assign_one();
402
403 delete [] exp;
404
405 return true;
406 }
407 else
408 return false;
409 }
410
411
412
413 #ifdef LIDIA_NAMESPACE
414 } // end of namespace LiDIA
415 # undef IN_NAMESPACE_LIDIA
416 #endif
417
418
419
420 #endif // LIDIA_PARTIAL_RELATION_H_GUARD_
421