1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Patrick Theobald (PT)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifndef LIDIA_BIGINT_MATRIX_ALGORITHMS_CC_GUARD_
20 #define LIDIA_BIGINT_MATRIX_ALGORITHMS_CC_GUARD_
21
22
23 #ifndef LIDIA_INFO_H_GUARD_
24 # include "LiDIA/info.h"
25 #endif
26 #ifndef LIDIA_DENSE_FP_MATRIX_KERNEL_H_GUARD_
27 # include "LiDIA/matrix/dense_fp_matrix_kernel.h"
28 #endif
29 #ifndef LIDIA_MODULAR_ARITHMETIC_H_GUARD_
30 # include "LiDIA/matrix/modular_arithmetic.h"
31 #endif
32 #include "LiDIA/precondition_error.h"
33
34
35
36 #ifdef LIDIA_NAMESPACE
37 # ifndef IN_NAMESPACE_LIDIA
38 namespace LiDIA {
39 # endif
40 #endif
41
42
43
44 //
45 // debug defines / error defines
46 //
47
48 extern const char *PRT;
49 extern const char *matrix_error_msg[];
50
51 #define DVALUE LDBL_MATRIX // Debug value
52 #define DMESSAGE "bigint_matrix_algorithms" // Debug message
53 #define EMESSAGE matrix_error_msg // Error message
54
55 //
56 // divide
57 //
58
59 template< class REP1, class REP2, class REP3 >
60 inline void bigint_matrix_algorithms< REP1, REP2, REP3 >::
divide(matrix<bigint> & RES,const matrix<bigint> & A,const bigint & k) const61 divide(matrix< bigint > &RES, const matrix< bigint > &A,
62 const bigint &k) const
63 {
64 register lidia_size_t j, i;
65 bigint TMP;
66
67 for (j = 0; j < A.rows; j++)
68 for (i = 0; i < A.columns; i++) {
69 LiDIA::divide(TMP, rep_modul2.member(A, j, i), k);
70 rep_modul1.sto(RES, j, i, TMP);
71 }
72 }
73
74
75
76 template< class REP1, class REP2, class REP3 >
77 inline void bigint_matrix_algorithms< REP1, REP2, REP3 >::
compwise_divide(matrix<bigint> & RES,const matrix<bigint> & A,const matrix<bigint> & B) const78 compwise_divide(matrix< bigint > &RES, const matrix< bigint > &A,
79 const matrix< bigint > &B) const
80 {
81 register lidia_size_t j, i;
82 bigint TMP;
83
84 for (j = 0; j < RES.rows; j++)
85 for (i = 0; i < RES.columns; i++) {
86 LiDIA::divide(TMP, rep_modul2.member(A, j, i), rep_modul3.member(B, j, i));
87 rep_modul1.sto(RES, j, i, TMP);
88 }
89 }
90
91
92
93 //
94 // remainder
95 //
96
97 template< class REP1, class REP2, class REP3 >
98 void bigint_matrix_algorithms< REP1, REP2, REP3 >::
remainder(matrix<bigint> & RES,const matrix<bigint> & M,const bigint & mod) const99 remainder(matrix< bigint > &RES, const matrix< bigint > & M,
100 const bigint & mod) const
101 {
102 RES.set_representation(M.bitfield.get_representation());
103
104 if (RES.bitfield.get_representation() == matrix_flags::dense_representation) {
105 lidia_size_t i, j;
106 bigint *REStmp, *Mtmp;
107
108 for (i = 0; i < M.rows; i++) {
109 REStmp = RES.value[i];
110 Mtmp = M.value[i];
111 for (j = 0; j < M.columns; j++)
112 LiDIA::best_remainder(REStmp[j], Mtmp[j], mod);
113 }
114 }
115 else {
116 lidia_size_t i, j;
117 bigint *REStmp;
118 lidia_size_t *REStmp2;
119 bigint *Mtmp;
120
121 for (i = 0; i < M.rows; i++) {
122 Mtmp = M.value[i];
123
124 if (RES.allocated[i] < M.value_counter[i]) {
125 delete[] RES.index[i];
126 delete[] RES.value[i];
127
128 REStmp = new bigint[M.value_counter[i]];
129 REStmp2 = new lidia_size_t[M.value_counter[i]];
130 RES.value_counter[i] = M.value_counter[i];
131 RES.allocated[i] = M.value_counter[i];
132 }
133 else {
134 REStmp = RES.value[i];
135 REStmp2 = RES.index[i];
136 RES.value_counter[i] = M.value_counter[i];
137 }
138
139 for (j = 0; j < M.value_counter[i]; j++) {
140 LiDIA::best_remainder(REStmp[j], Mtmp[j], mod);
141 REStmp2[j] = M.index[i][j];
142 }
143 RES.value[i] = REStmp;
144 RES.index[i] = REStmp2;
145 }
146 }
147 }
148
149
150
151 template< class REP1, class REP2, class REP3 >
152 inline void bigint_matrix_algorithms< REP1, REP2, REP3 >::
trans_remainder(matrix<bigint> & RES,const matrix<bigint> & M,const bigint & mod) const153 trans_remainder(matrix< bigint > &RES, const matrix< bigint > & M,
154 const bigint & mod) const
155 {
156
157 lidia_size_t i, j;
158 if (RES.bitfield.get_representation() == matrix_flags::dense_representation) {
159 bigint *REStmp;
160
161 for (i = 0; i < M.columns; i++) {
162 REStmp = RES.value[i];
163 for (j = 0; j < M.rows; j++)
164 LiDIA::best_remainder(REStmp[j], M.value[j][i], mod);
165 }
166 }
167 else {
168 bigint TMP;
169 for (i = 0; i < M.columns; i++) {
170 for (j = 0; j < M.rows; j++) {
171 LiDIA::best_remainder(TMP, M(j, i), mod);
172 RES.sto(i, j, TMP);
173 }
174 }
175 }
176 }
177
178
179
180 template< class REP1, class REP2, class REP3 >
181 void bigint_matrix_algorithms< REP1, REP2, REP3 >::
remainder(matrix<long> & RES,const matrix<bigint> & M,long mod) const182 remainder(matrix< long > &RES, const matrix< bigint > & M, long mod) const
183 {
184 RES.set_representation(M.bitfield.get_representation());
185
186 if (RES.bitfield.get_representation() == matrix_flags::dense_representation) {
187 lidia_size_t i, j;
188 long *REStmp;
189 bigint *Mtmp;
190
191 for (i = 0; i < M.rows; i++) {
192 REStmp = RES.value[i];
193 Mtmp = M.value[i];
194 for (j = 0; j < M.columns; j++)
195 LiDIA::best_remainder(REStmp[j], Mtmp[j], mod);
196 }
197 }
198 else {
199 lidia_size_t i, j;
200 long *REStmp;
201 lidia_size_t *REStmp2;
202 bigint *Mtmp;
203
204 for (i = 0; i < M.rows; i++) {
205 Mtmp = M.value[i];
206
207 if (RES.allocated[i] < M.value_counter[i]) {
208 delete[] RES.index[i];
209 delete[] RES.value[i];
210
211 REStmp = new long[M.value_counter[i]];
212 REStmp2 = new lidia_size_t[M.value_counter[i]];
213 RES.value_counter[i] = M.value_counter[i];
214 RES.allocated[i] = M.value_counter[i];
215 }
216 else {
217 REStmp = RES.value[i];
218 REStmp2 = RES.index[i];
219 RES.value_counter[i] = M.value_counter[i];
220 }
221
222 for (j = 0; j < M.value_counter[i]; j++) {
223 LiDIA::best_remainder(REStmp[j], Mtmp[j], mod);
224 REStmp2[j] = M.index[i][j];
225 }
226 RES.value[i] = REStmp;
227 RES.index[i] = REStmp2;
228 }
229 }
230 }
231
232
233
234 template< class REP1, class REP2, class REP3 >
235 inline void bigint_matrix_algorithms< REP1, REP2, REP3 >::
trans_remainder(matrix<long> & RES,const matrix<bigint> & M,long mod) const236 trans_remainder(matrix< long > &RES, const matrix< bigint > & M,
237 long mod) const
238 {
239 lidia_size_t i, j;
240 if (RES.bitfield.get_representation() == matrix_flags::dense_representation) {
241 long *REStmp;
242
243 for (i = 0; i < M.columns; i++) {
244 REStmp = RES.value[i];
245 for (j = 0; j < M.rows; j++)
246 LiDIA::best_remainder(REStmp[j], M.value[j][i], mod);
247 }
248 }
249 else {
250 long TMP;
251 for (i = 0; i < M.columns; i++) {
252 for (j = 0; j < M.rows; j++) {
253 LiDIA::best_remainder(TMP, M(j, i), mod);
254 RES.sto(i, j, TMP);
255 }
256 }
257 }
258 }
259
260
261
262 //
263 // Kernel
264 //
265
266 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
267 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
kernel1(matrix<bigint> & RES,const matrix<bigint> & A,const bigint & H) const268 kernel1(matrix< bigint > &RES, const matrix< bigint > & A, const bigint &H) const
269 {
270 debug_handler_l(DMESSAGE, "in member - function "
271 "kernel1(const matrix< bigint > &)", DVALUE + 8);
272
273 const modular_arithmetic< DRMK < bigint >,
274 dense_fp_matrix_kernel< long, MR < long > >,
275 dense_fp_matrix_kernel< bigint, MR < bigint > > > Dm_bigint_modul;
276
277 bigint *ZBAtmp, *Atmp;
278 register long i, j;
279 lidia_size_t c = A.columns;
280
281 // Step 1
282 lidia_size_t *linuz = Dm_bigint_modul.lininr1(A, H);
283 lidia_size_t r = linuz[0];
284 if (r == c) {
285 matrix< bigint > RES1(c, 1);
286 RES.assign(RES1);
287 delete[] linuz;
288 return;
289 }
290
291 // Step 2
292 matrix< bigint > ZBA(r, c);
293 for (i = 1; i <= r; i++) {
294 ZBAtmp = ZBA.value[i - 1];
295 Atmp = A.value[linuz[r - i + 1]];
296 for (j = 0; j < c; j++)
297 ZBAtmp[j].assign(Atmp[j]);
298 }
299 delete[] linuz;
300
301 // Step 3
302 matrix< bigint > TRANS(c, c);
303 ZBA.hnf(TRANS);
304
305 // Step 4
306 matrix< bigint > PART2(c, r);
307 if (RES.rows != c)
308 RES.set_no_of_rows(c);
309 if (RES.columns != c - r)
310 RES.set_no_of_columns(c - r);
311
312 //dmodul2.split_h(TRANS, RES, PART2);
313 dmodul.insert_at(RES, 0, 0, TRANS, 0, 0, RES.rows, RES.columns);
314 dmodul.insert_at(PART2, 0, 0, TRANS, 0, TRANS.columns - PART2.columns, PART2.rows, PART2.columns);
315 }
316
317
318
319 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
320 inline void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
kernel2(matrix<bigint> & RES,const matrix<bigint> & A) const321 kernel2(matrix< bigint > &RES, const matrix< bigint > & A) const
322 {
323 debug_handler_l(DMESSAGE, "in member - function "
324 "kernel2(const matrix< bigint > &)", DVALUE + 8);
325
326 register lidia_size_t i;
327 matrix< bigint > B = A;
328 B.hnf_havas(RES);
329
330 for (i = 0; i < A.columns && B.is_column_zero(i); i++);
331
332 if (i == 0) {
333 matrix< bigint > C(A.columns, 1);
334 RES.assign(C);
335 }
336 else
337 RES.set_no_of_columns(i);
338 }
339
340
341
342 //
343 // regular InvImage
344 //
345
346 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
347 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
reginvimage1(matrix<bigint> & RES,const matrix<bigint> & A,const matrix<bigint> & B,const bigint & H) const348 reginvimage1(matrix< bigint > &RES, const matrix< bigint > & A,
349 const matrix< bigint > & B, const bigint &H) const
350 {
351 //
352 // Task: C.reginvimage1(A, B);
353 // = > A * C.column(j) = g(j)*B.column(j), j = 0, ..., B.columns
354 // = > g(j) minimal
355 // Version: 2.0
356 //
357
358 debug_handler_l(DMESSAGE, "reginvimage1(const matrix< bigint > &, const matrix< bigint > &",
359 DVALUE + 8);
360
361 register long i, j;
362 bigint TMP, TMP1;
363
364 // Step 1
365 bigint DET;
366 A.det(DET, H);
367 if (DET == 0) {
368 precondition_error_handler(DET, "det(A)", "det(A) != 0",
369 "void matrix< bigint >::"
370 "reginvimage1(const matrix< bigint > & A, const matrix< bigint > & B)",
371 DMESSAGE, EMESSAGE[11]);
372 return;
373 }
374 matrix< bigint > ADJ;
375
376 const modular_arithmetic< DRMK < bigint >,
377 dense_fp_matrix_kernel< long, MR < long > >,
378 dense_fp_matrix_kernel< bigint, MR < bigint > > > Dm_bigint_modul;
379
380 Dm_bigint_modul.adj1(ADJ, A, H, DET);
381
382 // Step 2
383 matrix< bigint > PROD = ADJ * B;
384 bigint *u = new bigint[B.rows];
385 memory_handler(u, DMESSAGE, "reginvimage1 :: "
386 "Error in memory allocation (u)");
387 bigint *g, phi;
388
389 // Step 3
390 if (RES.rows != B.rows + 1)
391 RES.set_no_of_rows(B.rows + 1);
392 if (RES.columns != B.columns)
393 RES.set_no_of_columns(B.columns);
394 for (i = 0; i < B.columns; i++) {
395 for (j = 0; j < PROD.rows; j++)
396 u[j].assign(PROD.value[j][i]);
397 g = LiDIA::mgcd2(u, PROD.rows);
398 div_rem(phi, TMP, DET, gcd(g[0], DET));
399 if (phi.is_lt_zero())
400 phi.negate();
401
402 // Step 4
403 for (j = 0; j < PROD.rows; j++) {
404 LiDIA::multiply(TMP, phi, u[j]);
405 div_rem(RES.value[j][i], TMP1, TMP, DET);
406 }
407 RES.value[PROD.rows][i].assign(phi);
408 delete[] g;
409 }
410 delete[] u;
411 }
412
413
414
415 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
416 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
reginvimage2(matrix<bigint> & RES,const matrix<bigint> & A,const matrix<bigint> & B,const bigint & H) const417 reginvimage2(matrix< bigint > &RES, const matrix< bigint > & A,
418 const matrix< bigint > & B, const bigint &H) const
419 {
420 debug_handler_l(DMESSAGE, "in member - function "
421 "reginvimage2(const matrix< bigint > &, const matrix< bigint > &", DVALUE + 8);
422
423 register lidia_size_t i, j, len, oldlen;
424 bigint TMP, TMP1;
425
426 // Step 1
427 bigint DET;
428 A.det(DET, H);
429 if (DET == 0) {
430 precondition_error_handler(DET, "det(A)", "det(A) != 0",
431 "void matrix< bigint >::"
432 "reginvimage2(const matrix< bigint > & A, const matrix< bigint > & B)",
433 DMESSAGE, EMESSAGE[11]);
434 return;
435 }
436
437 // Step 2
438 oldlen = B.rows;
439 len = B.rows + 1;
440 LiDIA::multiply(RES, LiDIA::adj(A), B);
441
442 bigint *u = new bigint[len];
443 memory_handler(u, DMESSAGE, "reginvimage :: "
444 "Error in memory allocation (u)");
445 bigint phi;
446
447 // Step 3
448 RES.set_no_of_rows(len);
449 for (i = 0; i < B.columns; i++) {
450 for (j = 0; j < oldlen; j++)
451 u[j].assign(RES.value[j][i]);
452 u[oldlen].assign(DET);
453 LiDIA::mgcd2(TMP1, u, len);
454 div_rem(phi, TMP, DET, TMP1);
455 if (phi.is_lt_zero())
456 phi.negate();
457
458 // Step 4
459 for (j = 0; j < oldlen; j++) {
460 LiDIA::multiply(TMP, phi, u[j]);
461 div_rem(RES.value[j][i], TMP1, TMP, DET);
462 }
463 RES.value[oldlen][i].assign(phi);
464 }
465 delete[] u;
466 }
467
468
469
470 //
471 // Image
472 //
473
474 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
475 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
image1(matrix<bigint> & RES,const matrix<bigint> & A,const bigint & H) const476 image1(matrix< bigint > &RES, const matrix< bigint > & A, const bigint &H) const
477 {
478 bigint *ZBAtmp, *Atmp;
479 register long i, j;
480
481 // Step 1
482 const modular_arithmetic< DRMK < bigint >,
483 dense_fp_matrix_kernel< long, MR < long > >,
484 dense_fp_matrix_kernel< bigint, MR < bigint > > > Dm_bigint_modul;
485
486 lidia_size_t *v = Dm_bigint_modul.lininr1(A, H);
487 lidia_size_t RANG = v[0];
488
489 // Step 2, 3
490 matrix< bigint > ZBA(RANG, A.columns);
491 for (i = 1; i <= RANG; i++) {
492 ZBAtmp = ZBA.value[i - 1];
493 Atmp = A.value[v[RANG - i + 1]];
494 for (j = 0; j < A.columns; j++)
495 ZBAtmp[j].assign(Atmp[j]);
496 }
497 delete[] v;
498
499 // Step 4
500 matrix< bigint > TRANS(A.columns, A.columns);
501 ZBA.hnf(TRANS);
502
503 // Step 5
504 if (RES.rows != A.rows)
505 RES.set_no_of_rows(A.rows);
506 if (RES.columns != RANG)
507 RES.set_no_of_columns(RANG);
508
509 if (A.columns == RANG)
510 LiDIA::multiply(RES, A, TRANS);
511 else {
512 matrix< bigint > M(A.rows, A.columns);
513 LiDIA::multiply(M, A, TRANS);
514 matrix< bigint > PART1(RANG, A.columns - RANG);
515 //dmodul2.split_h(M, PART1, RES);
516 dmodul.insert_at(PART1, 0, 0, M, 0, 0, PART1.rows, PART1.columns);
517 dmodul.insert_at(RES, 0, 0, M, 0, M.columns - RES.columns, RES.rows, RES.columns);
518 }
519 }
520
521
522
523 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
524 inline void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
image2(matrix<bigint> & RES,const matrix<bigint> & A) const525 image2(matrix< bigint > &RES, const matrix< bigint > & A) const
526 {
527 register lidia_size_t i;
528 RES.assign(A);
529 RES.hnf_havas();
530
531 for (i = 0; i < RES.columns && RES.is_column_zero(i); i++);
532
533 if (i != 0)
534 if (i == RES.columns)
535 RES.set_no_of_columns(1);
536 else {
537 matrix< bigint > M(RES);
538 matrix< bigint > PART1(RES.rows, i);
539 RES.set_no_of_columns(RES.columns-i);
540 //dmodul2.split_h(M, PART1, RES);
541 dmodul.insert_at(PART1, 0, 0, M, 0, 0, PART1.rows, PART1.columns);
542 dmodul.insert_at(RES, 0, 0, M, 0, M.columns - RES.columns, RES.rows, RES.columns);
543 }
544 }
545
546
547
548 //
549 // InvImage
550 //
551
552 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
553 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
invimage(matrix<bigint> & RES,const matrix<bigint> & B,const bigint * b) const554 invimage(matrix< bigint > &RES, const matrix< bigint > & B, const bigint * b) const
555 {
556 if (b == NULL)
557 precondition_error_handler(PRT, "b", "b != NULL",
558 "void matrix< bigint >::"
559 "invimage(const matrix< bigint > & B, const bigint * b)",
560 DMESSAGE, EMESSAGE[1]);
561
562 register long i;
563 bigint *tmp;
564
565 // Step 1
566 matrix< bigint > A = B;
567 A.set_no_of_columns(B.columns + 1);
568 for (i = 0; i < B.rows; i++)
569 A.value[i][B.columns].assign(-b[i]);
570
571 kernel1(RES, A, LiDIA::hadamard(A));
572
573 // Step 2
574 if (RES.is_column_zero(0) || RES.is_row_zero(B.columns)) {
575 matrix< bigint > C(B.rows, 1);
576 RES.set_no_of_columns(1);
577 RES.assign(C);
578 return;
579 }
580
581 tmp = new bigint[RES.columns];
582 RES.get_row(tmp, B.columns);
583 bigint *g = LiDIA::mgcd2(tmp, RES.columns);
584 delete[] tmp;
585 if (g[0] > 1) {
586 matrix< bigint > C(B.rows, 1);
587 RES.set_no_of_columns(1);
588 RES.assign(C);
589 return;
590 }
591
592 // Step 3, 4
593 bigint *x = (RES) * &(g[1]);
594 delete[] g;
595
596 // Step 5
597 kernel1(RES, B, LiDIA::hadamard(B));
598 RES.set_no_of_columns(RES.columns + 1);
599 for (i = 0; i < RES.rows; i++)
600 RES.value[i][RES.columns-1].assign(x[i]);
601 delete[] x;
602 }
603
604
605
606 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
607 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
invimage(matrix<bigint> & RES,const matrix<bigint> & B,const math_vector<bigint> & b) const608 invimage(matrix< bigint > &RES, const matrix< bigint > & B, const math_vector< bigint > &b) const
609 {
610 debug_handler_l(DMESSAGE, "in member - function "
611 "invimage(const matrix< bigint > &, const math_vector< bigint > &)", DVALUE + 8);
612
613 register long i;
614 // Step 1
615 matrix< bigint > A = B;
616 A.set_no_of_columns(B.columns + 1);
617
618 if (b.size() != B.rows)
619 precondition_error_handler(b.size(), "b.size", "b.size == B.rows",
620 B.rows, "B.rows", "b.size == B.rows",
621 "void matrix< bigint >::"
622 "invimage(const matrix< bigint > & B, const math_vector< bigint > &b)",
623 DMESSAGE, EMESSAGE[1]);
624
625 bigint *tmp = b.get_data_address();
626 for (i = 0; i < B.rows; i++)
627 A.value[i][B.columns].assign(-tmp[i]);
628 kernel1(RES, A, LiDIA::hadamard(A));
629
630 // Step 2
631 if (RES.is_column_zero(0) || RES.is_row_zero(B.columns)) {
632 matrix< bigint > C(B.rows, 1);
633 RES.set_no_of_columns(1);
634 RES.assign(C);
635 return;
636 }
637
638 tmp = new bigint[RES.columns];
639 RES.get_row(tmp, B.columns);
640 bigint *g = LiDIA::mgcd2(tmp, RES.columns);
641
642 delete[] tmp;
643 if (g[0] > 1) {
644 matrix< bigint > C(B.rows, 1);
645 RES.set_no_of_columns(1);
646 RES.assign(C);
647 return;
648 }
649
650 // Step 3, 4
651 bigint *x = (RES) * &(g[1]);
652 delete[] g;
653
654 // Step 5
655 kernel1(RES, B, LiDIA::hadamard(B));
656 RES.set_no_of_columns(RES.columns + 1);
657 for (i = 0; i < RES.rows; i++)
658 RES.value[i][RES.columns-1].assign(x[i]);
659 delete[] x;
660 }
661
662
663
664 //
665 // Smith normal form
666 //
667
668 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
669 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_hartley(matrix<bigint> & RES) const670 snf_hartley(matrix< bigint > &RES) const
671 {
672 debug_handler_l(DMESSAGE, "in member - function "
673 "snf_hartley()", DVALUE + 8);
674
675 register lidia_size_t startr, startc, TEILBARKEIT;
676 bigint TMP1, TMP2;
677 bigint *tmp;
678 register lidia_size_t xpivot, ypivot, i, j, z;
679
680 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
681
682 // pivot: first non zero
683
684 xpivot = -1;
685 ypivot = -1;
686 for (i = startr; i < RES.rows; i++) {
687 tmp = RES.value[i];
688 for (j = startc; j < RES.columns; j++)
689 if (!tmp[j].is_zero()) {
690 xpivot = i;
691 ypivot = j;
692 i = RES.rows;
693 j = RES.columns;
694 }
695 }
696
697 if (xpivot != -1) {
698 // swap to diagonalposition
699 RES.swap_rows(startr, xpivot);
700 RES.swap_columns(startc, ypivot);
701
702 TEILBARKEIT = 0;
703
704 while (TEILBARKEIT == 0) {
705 TEILBARKEIT = 1;
706
707 // mgcd computation for row
708 for (i = startc + 1; i < RES.columns; i++)
709 if (RES.value[startr][i] % RES.value[startr][startc] != 0) {
710 div_rem(TMP1, TMP2, RES.value[startr][i], RES.value[startr][startc]);
711 for (j = startr; j < RES.rows; j++) {
712 tmp = RES.value[j];
713 LiDIA::multiply(TMP2, tmp[startc], TMP1);
714 LiDIA::subtract(tmp[i], tmp[i], TMP2);
715 }
716 RES.swap_columns(startc, i);
717 i = startc;
718 }
719
720 // mgcd computation for column
721 for (i = startr + 1; i < RES.rows; i++)
722 if (RES.value[i][startc] % RES.value[startr][startc] != 0) {
723 div_rem(TMP1, TMP2, RES.value[i][startc], RES.value[startr][startc]);
724 for (j = startc; j < RES.columns; j++) {
725 LiDIA::multiply(TMP2, RES.value[startr][j], TMP1);
726 LiDIA::subtract(RES.value[i][j], RES.value[i][j], TMP2);
727 }
728 TEILBARKEIT = 0; //perhaps
729 RES.swap_rows(i, startr);
730 i = startr;
731 }
732 }
733
734 // row elimination
735 for (i = startc + 1; i < RES.columns; i++) {
736 div_rem(TMP1, TMP2, RES.value[startr][i], RES.value[startr][startc]);
737 for (j = startr; j < RES.rows; j++) {
738 tmp = RES.value[j];
739 LiDIA::multiply(TMP2, tmp[startc], TMP1);
740 LiDIA::subtract(tmp[i], tmp[i], TMP2);
741 }
742 }
743
744 // column elimination
745 for (i = startr + 1; i < RES.rows; i++) {
746 div_rem(TMP1, TMP2, RES.value[i][startc], RES.value[startr][startc]);
747 for (j = startc; j < RES.columns; j++) {
748 LiDIA::multiply(TMP2, RES.value[startr][j], TMP1);
749 LiDIA::subtract(RES.value[i][j], RES.value[i][j], TMP2);
750 }
751 }
752
753 // modulo test
754 for (i = startr + 1; i < RES.rows; i++)
755 for (j = startc + 1; j < RES.columns; j++)
756 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
757 for (z = 0; z < RES.columns; z++)
758 LiDIA::add(RES.value[startr][z], RES.value[startr][z], RES.value[i][z]);
759 i = RES.rows;
760 j = RES.columns;
761 startc = 0;
762 startr = 0;
763 }
764 }
765 }
766
767 // diagonal >= 0
768 for (i = 0; i < RES.columns && i < RES.rows; i++)
769 if (RES.value[i][i].is_lt_zero())
770 RES.value[i][i].negate();
771 }
772
773
774
775 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
776 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_hartley(matrix<bigint> & RES,matrix<bigint> & T1,matrix<bigint> & T2) const777 snf_hartley(matrix< bigint > &RES, matrix< bigint > & T1, matrix< bigint > & T2) const
778 {
779 debug_handler_l(DMESSAGE, "in member - function "
780 "snf_hartley(matrix< bigint > &, matrix< bigint > &)", DVALUE + 8);
781
782 register lidia_size_t startr, startc, TEILBARKEIT;
783 bigint TMP1, TMP2;
784 lidia_size_t xpivot, ypivot;
785 register lidia_size_t i, j, z;
786
787 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
788
789 // pivot: first non zero
790 xpivot = -1;
791 ypivot = -1;
792 for (i = startr; i < RES.rows; i++)
793 for (j = startc; j < RES.columns; j++)
794 if (RES.value[i][j] != 0) {
795 xpivot = i;
796 ypivot = j;
797 i = RES.rows;
798 j = RES.columns;
799 }
800 if (xpivot != -1) {
801 // swap to diagonalposition
802 RES.swap_rows(startr, xpivot);
803 T1.swap_rows(startr, xpivot);
804 RES.swap_columns(startc, ypivot);
805 T2.swap_columns(startc, ypivot);
806
807 TEILBARKEIT = 0;
808
809 while (TEILBARKEIT == 0) {
810 TEILBARKEIT = 1;
811
812 // mgcd computation for row
813 for (i = startc + 1; i < RES.columns; i++)
814 if (RES.value[startr][i] % RES.value[startr][startc] != 0) {
815 div_rem(TMP1, TMP2, RES.value[startr][i], RES.value[startr][startc]);
816 for (j = startr; j < RES.rows; j++) {
817 LiDIA::multiply(TMP2, RES.value[j][startc], TMP1);
818 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
819 }
820 for (j = 0; j < RES.columns; j++) {
821 LiDIA::multiply(TMP2, T2.value[j][startc], TMP1);
822 LiDIA::subtract(T2.value[j][i], T2.value[j][i], TMP2);
823 }
824 RES.swap_columns(startc, i);
825 T2.swap_columns(startc, i);
826 i = startc;
827 }
828
829 // mgcd computation for column
830 for (i = startr + 1; i < RES.rows; i++)
831 if (RES.value[i][startc] % RES.value[startr][startc] != 0) {
832 div_rem(TMP1, TMP2, RES.value[i][startc], RES.value[startr][startc]);
833 for (j = startc; j < RES.columns; j++) {
834 LiDIA::multiply(TMP2, RES.value[startr][j], TMP1);
835 LiDIA::subtract(RES.value[i][j], RES.value[i][j], TMP2);
836 }
837 for (j = 0; j < RES.rows; j++) {
838 LiDIA::multiply(TMP2, T1.value[startr][j], TMP1);
839 LiDIA::subtract(T1.value[i][j], T1.value[i][j], TMP2);
840 }
841 TEILBARKEIT = 0;
842 RES.swap_rows(i, startr);
843 T1.swap_rows(i, startr);
844 i = startr;
845 }
846 }
847
848 // row elimination
849 for (i = startc + 1; i < RES.columns; i++) {
850 div_rem(TMP1, TMP2, RES.value[startr][i], RES.value[startr][startc]);
851 for (j = startr; j < RES.rows; j++) {
852 LiDIA::multiply(TMP2, RES.value[j][startc], TMP1);
853 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
854 }
855 for (j = 0; j < RES.columns; j++) {
856 LiDIA::multiply(TMP2, T2.value[j][startc], TMP1);
857 LiDIA::subtract(T2.value[j][i], T2.value[j][i], TMP2);
858 }
859 }
860
861 // column elimination
862 for (i = startr + 1; i < RES.rows; i++) {
863 div_rem(TMP1, TMP2, RES.value[i][startc], RES.value[startr][startc]);
864 for (j = startc; j < RES.columns; j++) {
865 LiDIA::multiply(TMP2, RES.value[startr][j], TMP1);
866 LiDIA::subtract(RES.value[i][j], RES.value[i][j], TMP2);
867 }
868 for (j = 0; j < RES.rows; j++) {
869 LiDIA::multiply(TMP2, T1.value[startr][j], TMP1);
870 LiDIA::subtract(T1.value[i][j], T1.value[i][j], TMP2);
871 }
872 }
873
874 // modulo test
875 for (i = startr + 1; i < RES.rows; i++)
876 for (j = startc + 1; j < RES.columns; j++)
877 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
878 for (z = 0; z < RES.columns; z++)
879 LiDIA::add(RES.value[startr][z], RES.value[startr][z], RES.value[i][z]);
880 for (z = 0; z < RES.rows; z++)
881 LiDIA::add(T1.value[startr][z], T1.value[startr][z], T1.value[i][z]);
882 i = RES.rows;
883 j = RES.columns;
884 startc = 0;
885 startr = 0;
886 }
887 }
888 }
889
890 // diagonal >= 0
891 for (i = 0; i < RES.rows && i < RES.columns; i++)
892 if (RES.value[i][i] < 0) {
893 RES.value[i][i].negate();
894 for (z = 0; z < RES.columns; z++)
895 T2.value[z][i].negate();
896 }
897 }
898
899
900
901 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
902 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_simple(matrix<bigint> & RES) const903 snf_simple(matrix< bigint > &RES) const
904 {
905 debug_handler_l(DMESSAGE, "in member - function "
906 "snf_simple()", DVALUE + 8);
907 bigint PIVOT, TMP1, TMP2;
908 bigint *tmp = NULL, *deltmp;
909
910 matrix< bigint > TR1(RES.rows, RES.rows);
911 matrix< bigint > TR2(RES.columns, RES.columns);
912
913 bigint *REM;
914 register lidia_size_t startr, startc, pivot, i, j, z, TEILBARKEIT;
915
916 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
917
918 // pivot: first non zero
919 pivot = -1;
920 for (i = startr; i < RES.rows; i++)
921 for (j = startc; j < RES.columns; j++)
922 if (RES.value[i][j] != 0) {
923 pivot = i;
924 i = RES.rows;
925 j = RES.columns;
926 }
927
928 if (pivot != -1) {
929 // swap pivot in actual row
930 RES.swap_rows(startr, pivot);
931
932 TEILBARKEIT = 0;
933
934 while (TEILBARKEIT == 0) {
935 TEILBARKEIT = 1;
936
937 // mgcd computation and row elimination
938 deltmp = new bigint[RES.columns];
939 RES.get_row(deltmp, startr);
940 REM = mgcd2(TR2, deltmp, RES.columns);
941 delete[] deltmp;
942
943 delete[] REM;
944 TR2 = TR2.trans();
945 LiDIA::multiply(RES, RES, TR2);
946
947 tmp = RES.value[startr];
948 for (i = 0; tmp[i].is_zero() && i < RES.columns; i++);
949 RES.swap_columns(startc, i);
950
951 // mgcd computation and column elimination
952 deltmp = new bigint[RES.rows];
953 RES.get_column(deltmp, startc);
954 REM = mgcd2(TR1, deltmp, RES.rows);
955 delete[] deltmp;
956
957 delete[] REM;
958 LiDIA::multiply(RES, TR1, RES);
959
960 for (i = 0; RES.value[i][startc].is_zero() && i < RES.rows; i++);
961 RES.swap_rows(startr, i);
962
963 // control: row == 0
964 tmp = RES.value[startr];
965 for (i = startc+1; tmp[i].is_zero() && i < RES.columns; i++);
966 if (i != RES.columns)
967 TEILBARKEIT = 0;
968 }
969
970 // modulo test
971 for (i = startr; i < RES.rows; i++)
972 for (j = startc + 1; j < RES.columns; j++)
973 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
974 if (i != startr)
975 for (z = 0; z < RES.columns; z++)
976 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
977 i = RES.rows;
978 j = RES.columns;
979 startc--;
980 startr--;
981 }
982 }
983 }
984
985 // diagonal >= 0
986 for (i = 0; i < RES.rows && i < RES.columns; i++)
987 if (RES.value[i][i] < 0)
988 RES.value[i][i].negate();
989 }
990
991
992
993 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
994 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_simple(matrix<bigint> & RES,matrix<bigint> & T1,matrix<bigint> & T2) const995 snf_simple(matrix< bigint > &RES, matrix< bigint > & T1, matrix< bigint > & T2) const
996 {
997 debug_handler_l(DMESSAGE, "in member - function "
998 "snf_simple(matrix< bigint > &, matrix< bigint > &)", DVALUE + 8);
999 bigint PIVOT, TMP1, TMP2;
1000 bigint *tmp = NULL, *deltmp;
1001
1002 matrix< bigint > TR1 = T1;
1003 matrix< bigint > TR2 = T2;
1004 bigint *REM;
1005 register lidia_size_t startr, startc, pivot, i, j, z, TEILBARKEIT;
1006
1007 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1008
1009 // pivot: first non zero
1010 pivot = -1;
1011 for (i = startr; i < RES.rows; i++)
1012 for (j = startc; j < RES.columns; j++)
1013 if (RES.value[i][j] != 0) {
1014 pivot = i;
1015 i = RES.rows;
1016 j = RES.columns;
1017 }
1018
1019 if (pivot != -1) {
1020 // swap pivot in actual row
1021 RES.swap_rows(startr, pivot);
1022 T1.swap_rows(startr, pivot);
1023
1024 TEILBARKEIT = 0;
1025
1026 while (TEILBARKEIT == 0) {
1027 TEILBARKEIT = 1;
1028
1029 // mgcd computation and row elimination
1030 deltmp = new bigint[RES.columns];
1031 RES.get_row(deltmp, startr);
1032 REM = mgcd2(TR2, deltmp, RES.columns);
1033 delete[] deltmp;
1034
1035 delete[] REM;
1036 TR2 = TR2.trans();
1037 LiDIA::multiply(RES, RES, TR2);
1038 LiDIA::multiply(T2, T2, TR2);
1039
1040 tmp = RES.value[startr];
1041 for (i = 0; tmp[i].is_zero() && i < RES.columns; i++);
1042 RES.swap_columns(startc, i);
1043 T2.swap_columns(startc, i);
1044
1045 // mgcd computation and column elimination
1046 deltmp = new bigint[RES.rows];
1047 RES.get_column(deltmp, startc);
1048 REM = mgcd2(TR1, deltmp, RES.rows);
1049 delete[] deltmp;
1050
1051 delete[] REM;
1052 LiDIA::multiply(RES, TR1, RES);
1053 LiDIA::multiply(T1, TR1, T1);
1054
1055 for (i = 0; RES.value[i][startc].is_zero() && i < RES.rows; i++);
1056 RES.swap_rows(startr, i);
1057 T1.swap_rows(startr, i);
1058
1059 // control: row == 0
1060 tmp = RES.value[startr];
1061 for (i = startc+1; tmp[i].is_zero() && i < RES.columns; i++);
1062 if (i != RES.columns)
1063 TEILBARKEIT = 0;
1064 }
1065
1066 // modulo test
1067 for (i = startr; i < RES.rows; i++)
1068 for (j = startc + 1; j < RES.columns; j++)
1069 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
1070 if (i != startr) {
1071 for (z = 0; z < RES.columns; z++)
1072 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
1073 for (z = 0; z < RES.rows; z++)
1074 LiDIA::add(T1.value[startr][z], T1.value[startr][z], T1.value[i][z]);
1075 }
1076 i = RES.rows;
1077 j = RES.columns;
1078 startc--;
1079 startr--;
1080 }
1081 }
1082 }
1083
1084 // diagonal >= 0
1085 for (i = 0; i < RES.rows && i < RES.columns; i++)
1086 if (RES.value[i][i] < 0) {
1087 RES.value[i][i].negate();
1088 for (z = 0; z < RES.columns; z++)
1089 T2.value[z][i].negate();
1090 }
1091 }
1092
1093
1094
1095 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
1096 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_havas(matrix<bigint> & RES) const1097 snf_havas(matrix< bigint > &RES) const
1098 {
1099 debug_handler_l(DMESSAGE, "in member - function "
1100 "snf_havas()", DVALUE + 8);
1101
1102 register lidia_size_t i, j, z, index;
1103 bigint PIVOT;
1104 bigint *tmp = NULL;
1105
1106 register lidia_size_t startr, startc, xpivot, ypivot, SW, TEILBARKEIT;
1107 bigint TMP1, TMP2;
1108
1109 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1110 // pivot: first non zero
1111 xpivot = -1;
1112 ypivot = -1;
1113 for (i = startr; i < RES.rows; i++)
1114 for (j = startc; j < RES.columns; j++)
1115 if (!RES.value[i][j].is_zero()) {
1116 xpivot = i;
1117 ypivot = j;
1118 i = RES.rows;
1119 j = RES.columns;
1120 }
1121
1122 if (xpivot != -1) {
1123 // swap to actual row
1124 RES.swap_rows(startr, xpivot);
1125
1126 index = ypivot;
1127
1128 TEILBARKEIT = 0;
1129
1130 while (TEILBARKEIT == 0) {
1131 TEILBARKEIT = 1;
1132
1133 // gcd2(row(startr), columns, TR2);
1134 tmp = RES.value[startr];
1135 do
1136 {
1137 SW = 0;
1138 for (i = 0; i < RES.columns; i++)
1139 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
1140 index = i;
1141
1142 for (i = 0; i < RES.columns; i++)
1143 if (i != index && !tmp[i].is_zero()) {
1144 SW = 1;
1145 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
1146 for (j = 0; j < RES.rows; j++) {
1147 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
1148 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
1149 }
1150 }
1151 }
1152 while (SW == 1);
1153
1154 for (i = 0; RES.value[startr][i].is_zero() && i < RES.columns; i++);
1155 RES.swap_columns(startc, i);
1156
1157 // mgcd2(column(startc), rows, TR1);
1158 index = startr; // no index search
1159 do
1160 {
1161 SW = 0;
1162 for (i = 0; i < RES.rows; i++)
1163 if ((abs_compare(RES.value[index][startc], RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
1164 index = i;
1165
1166 for (i = 0; i < RES.rows; i++)
1167 if ((i != index) && !RES.value[i][startc].is_zero()) {
1168 SW = 1;
1169 tmp = RES.value[i];
1170 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
1171 for (j = 0; j < RES.columns; j++) {
1172 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
1173 LiDIA::subtract(tmp[j], tmp[j], TMP2);
1174 }
1175 }
1176 }
1177 while (SW == 1);
1178
1179 for (i = 0; RES.value[i][startc].is_zero() && i < RES.rows; i++);
1180 RES.swap_rows(startr, i);
1181
1182 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
1183 if (index != RES.columns)
1184 TEILBARKEIT = 0;
1185 index = startc;
1186 }
1187
1188 // modulo test
1189 for (i = startr; i < RES.rows; i++)
1190 for (j = startc + 1; j < RES.columns; j++)
1191 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
1192 if (i != startr)
1193 for (z = 0; z < RES.columns; z++)
1194 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
1195 i = RES.rows;
1196 j = RES.columns;
1197 startc--;
1198 startr--;
1199 }
1200 }
1201 }
1202
1203 // diagonal >= 0
1204 for (i = 0; i < RES.rows && i < RES.columns; i++)
1205 if (RES.value[i][i] < 0)
1206 RES.value[i][i].negate();
1207 }
1208
1209
1210
1211 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
1212 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_havas(matrix<bigint> & RES,matrix<bigint> & T1,matrix<bigint> & T2) const1213 snf_havas(matrix< bigint > &RES, matrix< bigint > & T1, matrix< bigint > & T2) const
1214 {
1215 debug_handler_l(DMESSAGE, "in member - function "
1216 "snf_havas(matrix< bigint > &, matrix< bigint > &)", DVALUE + 8);
1217
1218 register lidia_size_t i, j, z, index;
1219 bigint PIVOT;
1220 bigint *tmp = NULL;
1221
1222 register lidia_size_t startr, startc, xpivot, ypivot, SW, TEILBARKEIT;
1223 bigint TMP1, TMP2;
1224
1225 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1226 // pivot: first non zero
1227 xpivot = -1;
1228 ypivot = -1;
1229 for (i = startr; i < RES.rows; i++)
1230 for (j = startc; j < RES.columns; j++)
1231 if (!RES.value[i][j].is_zero()) {
1232 xpivot = i;
1233 ypivot = j;
1234 i = RES.rows;
1235 j = RES.columns;
1236 }
1237
1238 if (xpivot != -1) {
1239 // swap to actual row
1240 RES.swap_rows(startr, xpivot);
1241 T1.swap_rows(startr, xpivot);
1242
1243 index = ypivot;
1244
1245 TEILBARKEIT = 0;
1246
1247 while (TEILBARKEIT == 0) {
1248 TEILBARKEIT = 1;
1249
1250 // gcd2(row(startr), columns, TR2);
1251 tmp = RES.value[startr];
1252 do
1253 {
1254 SW = 0;
1255 for (i = 0; i < RES.columns; i++)
1256 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
1257 index = i;
1258
1259 for (i = 0; i < RES.columns; i++)
1260 if (i != index && !tmp[i].is_zero()) {
1261 SW = 1;
1262 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
1263 for (j = 0; j < RES.rows; j++) {
1264 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
1265 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
1266 }
1267 for (j = 0; j < RES.columns; j++) {
1268 LiDIA::multiply(TMP2, T2.value[j][index], TMP1);
1269 LiDIA::subtract(T2.value[j][i], T2.value[j][i], TMP2);
1270 }
1271 }
1272 }
1273 while (SW == 1);
1274
1275 for (i = 0; RES.value[startr][i].is_zero() && i < RES.columns; i++);
1276 RES.swap_columns(startc, i);
1277 T2.swap_columns(startc, i);
1278
1279 // mgcd2(column(startc), rows, TR1);
1280 index = startr; // no index search
1281 do
1282 {
1283 SW = 0;
1284 for (i = 0; i < RES.rows; i++)
1285 if ((abs_compare(RES.value[index][startc], RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
1286 index = i;
1287
1288 for (i = 0; i < RES.rows; i++)
1289 if ((i != index) && !RES.value[i][startc].is_zero()) {
1290 SW = 1;
1291 tmp = RES.value[i];
1292 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
1293 for (j = 0; j < RES.columns; j++) {
1294 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
1295 LiDIA::subtract(tmp[j], tmp[j], TMP2);
1296 }
1297 for (j = 0; j < RES.rows; j++) {
1298 LiDIA::multiply(TMP2, T1.value[index][j], TMP1);
1299 LiDIA::subtract(T1.value[i][j], T1.value[i][j], TMP2);
1300 }
1301 }
1302 }
1303 while (SW == 1);
1304
1305 for (i = 0; RES.value[i][startc].is_zero() && i < RES.rows; i++);
1306 RES.swap_rows(startr, i);
1307 T1.swap_rows(startr, i);
1308
1309 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
1310 if (index != RES.columns)
1311 TEILBARKEIT = 0;
1312 index = startc;
1313 }
1314
1315 // modulo test
1316 for (i = startr; i < RES.rows; i++)
1317 for (j = startc + 1; j < RES.columns; j++)
1318 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
1319 if (i != startr) {
1320 for (z = 0; z < RES.columns; z++)
1321 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
1322 for (z = 0; z < RES.rows; z++)
1323 LiDIA::add(T1.value[startr][z], T1.value[startr][z], T1.value[i][z]);
1324 }
1325 i = RES.rows;
1326 j = RES.columns;
1327 startc--;
1328 startr--;
1329 }
1330 }
1331 }
1332
1333 // diagonal >= 0
1334 for (i = 0; i < RES.rows && i < RES.columns; i++)
1335 if (RES.value[i][i] < 0) {
1336 RES.value[i][i].negate();
1337 for (z = 0; z < RES.columns; z++)
1338 T2.value[z][i].negate();
1339 }
1340 }
1341
1342
1343
1344 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
1345 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_mult(matrix<bigint> & RES,long art) const1346 snf_mult(matrix< bigint > &RES, long art) const
1347 {
1348 debug_handler_l(DMESSAGE, "in member - function "
1349 "snf_mult(matrix< bigint > &, matrix< bigint > &)", DVALUE + 8);
1350 register lidia_size_t i, j, z, index, SW;
1351 bigint TMP1, TMP2;
1352 bigint *tmp = NULL;
1353
1354 register lidia_size_t startr, startc, xpivot, ypivot, TEILBARKEIT;
1355 bigint ROW, COLUMN, PIVOT, NORM;
1356
1357 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1358
1359 // pivotselection: minimale C * R norm
1360 xpivot = -1;
1361 ypivot = -1;
1362 PIVOT.assign_zero();
1363 for (i = startr; i < RES.rows; i++)
1364 for (j = startc; j < RES.columns; j++) {
1365 if (PIVOT == abs(RES.value[i][j])) {
1366 RES.row_norm(ROW, i, art);
1367 RES.column_norm(COLUMN, j, art);
1368 LiDIA::multiply(TMP1, ROW, COLUMN);
1369 if (TMP1 < NORM) {
1370 NORM.assign(TMP1);
1371 PIVOT.assign(abs(RES.value[i][j]));
1372 xpivot = i;
1373 ypivot = j;
1374 }
1375 }
1376
1377 if (((abs_compare(PIVOT, RES.value[i][j]) > 0) && !RES.value[i][j].is_zero()) || PIVOT.is_zero()) {
1378 PIVOT.assign(abs(RES.value[i][j]));
1379 RES.row_norm(ROW, i, art);
1380 RES.column_norm(COLUMN, j, art);
1381 LiDIA::multiply(NORM, ROW, COLUMN);
1382 xpivot = i;
1383 ypivot = j;
1384 }
1385 }
1386
1387 if (!PIVOT.is_zero()) {
1388
1389 // swap to actual row
1390 RES.swap_rows(startr, xpivot);
1391
1392 index = ypivot;
1393
1394 TEILBARKEIT = 0;
1395
1396 while (TEILBARKEIT == 0) {
1397 TEILBARKEIT = 1;
1398
1399 // gcd2(row(startr), columns, TR2);
1400 tmp = RES.value[startr];
1401 do
1402 {
1403 SW = 0;
1404 for (i = 0; i < RES.columns; i++)
1405 if ((i != index) && !tmp[i].is_zero()) {
1406 SW = 1;
1407 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
1408 for (j = 0; j < RES.rows; j++) {
1409 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
1410 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
1411 }
1412 }
1413
1414 for (i = 0; i < RES.columns; i++)
1415 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
1416 index = i;
1417 }
1418 while (SW == 1);
1419
1420 for (i = 0; RES.value[startr][i].is_zero(); i++);
1421 RES.swap_columns(startc, i);
1422
1423 // mgcd2(column(startc), rows, TR1);
1424 index = startr;
1425 do
1426 {
1427 SW = 0;
1428 for (i = 0; i < RES.rows; i++)
1429 if ((abs_compare(RES.value[index][startc], RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
1430 index = i;
1431
1432 for (i = 0; i < RES.rows; i++)
1433 if ((i != index) && !RES.value[i][startc].is_zero()) {
1434 SW = 1;
1435 tmp = RES.value[i];
1436 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
1437 for (j = 0; j < RES.columns; j++) {
1438 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
1439 LiDIA::subtract(tmp[j], tmp[j], TMP2);
1440 }
1441 }
1442 }
1443 while (SW == 1);
1444
1445 for (i = 0; RES.value[i][startc].is_zero(); i++);
1446 RES.swap_rows(startr, i);
1447
1448 tmp = RES.value[startr];
1449 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
1450 if (index != RES.columns)
1451 TEILBARKEIT = 0;
1452
1453 index = startr;
1454 }
1455
1456 // modulo test
1457 for (i = startr; i < RES.rows; i++)
1458 for (j = startc + 1; j < RES.columns; j++)
1459 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
1460 if (i != startr)
1461 for (z = 0; z < RES.columns; z++)
1462 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
1463 i = RES.rows;
1464 j = RES.columns;
1465 startc--;
1466 startr--;
1467 }
1468 }
1469 }
1470
1471 // diagonal >= 0
1472 for (i = 0; i < RES.rows && i < RES.columns; i++)
1473 if (RES.value[i][i] < 0)
1474 RES.value[i][i].negate();
1475 }
1476
1477
1478
1479 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
1480 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_mult(matrix<bigint> & RES,matrix<bigint> & T1,matrix<bigint> & T2,long art) const1481 snf_mult(matrix< bigint > &RES, matrix< bigint > & T1, matrix< bigint > & T2, long art) const
1482 {
1483 debug_handler_l(DMESSAGE, "in member - function "
1484 "snf_mult(matrix< bigint > &, matrix< bigint > &)", DVALUE + 8);
1485
1486 register lidia_size_t i, j, z, index, SW;
1487 bigint TMP1, TMP2;
1488 bigint *tmp = NULL;
1489
1490 register lidia_size_t startr, startc, xpivot, ypivot, TEILBARKEIT;
1491 bigint ROW, COLUMN, PIVOT, NORM;
1492
1493 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1494
1495 // pivotselection: minimale C * R norm
1496 xpivot = -1;
1497 ypivot = -1;
1498 PIVOT.assign_zero();
1499 for (i = startr; i < RES.rows; i++)
1500 for (j = startc; j < RES.columns; j++) {
1501 if (PIVOT == abs(RES.value[i][j])) {
1502 RES.row_norm(ROW, i, art);
1503 RES.column_norm(COLUMN, j, art);
1504 LiDIA::multiply(TMP1, ROW, COLUMN);
1505 if (TMP1 < NORM) {
1506 NORM.assign(TMP1);
1507 PIVOT.assign(abs(RES.value[i][j]));
1508 xpivot = i;
1509 ypivot = j;
1510 }
1511 }
1512
1513 if (((abs_compare(PIVOT, RES.value[i][j]) > 0) && !RES.value[i][j].is_zero()) || PIVOT.is_zero()) {
1514 PIVOT.assign(abs(RES.value[i][j]));
1515 RES.row_norm(ROW, i, art);
1516 RES.column_norm(COLUMN, j, art);
1517 LiDIA::multiply(NORM, ROW, COLUMN);
1518 xpivot = i;
1519 ypivot = j;
1520 }
1521 }
1522
1523 if (!PIVOT.is_zero()) {
1524
1525 // swap to actual row
1526 RES.swap_rows(startr, xpivot);
1527 T1.swap_rows(startr, xpivot);
1528
1529 index = ypivot;
1530
1531 TEILBARKEIT = 0;
1532
1533 while (TEILBARKEIT == 0) {
1534 TEILBARKEIT = 1;
1535
1536 // gcd2(row(startr), columns, TR2);
1537 tmp = RES.value[startr];
1538 do
1539 {
1540 SW = 0;
1541 for (i = 0; i < RES.columns; i++)
1542 if ((i != index) && !tmp[i].is_zero()) {
1543 SW = 1;
1544 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
1545 for (j = 0; j < RES.rows; j++) {
1546 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
1547 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
1548 }
1549 for (j = 0; j < RES.columns; j++) {
1550 LiDIA::multiply(TMP2, T2.value[j][index], TMP1);
1551 LiDIA::subtract(T2.value[j][i], T2.value[j][i], TMP2);
1552 }
1553 }
1554
1555 for (i = 0; i < RES.columns; i++)
1556 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
1557 index = i;
1558 }
1559 while (SW == 1);
1560
1561 for (i = 0; RES.value[startr][i].is_zero(); i++);
1562 RES.swap_columns(startc, i);
1563 T2.swap_columns(startc, i);
1564
1565 // mgcd2(column(startc), rows, TR1);
1566 index = startr;
1567 do
1568 {
1569 SW = 0;
1570 for (i = 0; i < RES.rows; i++)
1571 if ((abs_compare(RES.value[index][startc], RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
1572 index = i;
1573
1574 for (i = 0; i < RES.rows; i++)
1575 if ((i != index) && !RES.value[i][startc].is_zero()) {
1576 SW = 1;
1577 tmp = RES.value[i];
1578 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
1579 for (j = 0; j < RES.columns; j++) {
1580 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
1581 LiDIA::subtract(tmp[j], tmp[j], TMP2);
1582 }
1583 for (j = 0; j < RES.rows; j++) {
1584 LiDIA::multiply(TMP2, T1.value[index][j], TMP1);
1585 LiDIA::subtract(T1.value[i][j], T1.value[i][j], TMP2);
1586 }
1587 }
1588 }
1589 while (SW == 1);
1590
1591 for (i = 0; RES.value[i][startc].is_zero(); i++);
1592 RES.swap_rows(startr, i);
1593 T1.swap_rows(startr, i);
1594
1595 tmp = RES.value[startr];
1596 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
1597 if (index != RES.columns)
1598 TEILBARKEIT = 0;
1599
1600 index = startr;
1601 }
1602
1603 // modulo test
1604 for (i = startr; i < RES.rows; i++)
1605 for (j = startc + 1; j < RES.columns; j++)
1606 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
1607 if (i != startr) {
1608 for (z = 0; z < RES.columns; z++)
1609 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
1610 for (z = 0; z < RES.rows; z++)
1611 LiDIA::add(T1.value[startr][z], T1.value[startr][z], T1.value[i][z]);
1612 }
1613 i = RES.rows;
1614 j = RES.columns;
1615 startc--;
1616 startr--;
1617 }
1618 }
1619 }
1620
1621 // diagonal >= 0
1622 for (i = 0; i < RES.rows && i < RES.columns; i++)
1623 if (RES.value[i][i] < 0) {
1624 RES.value[i][i].negate();
1625 for (z = 0; z < RES.columns; z++)
1626 T2.value[z][i].negate();
1627 }
1628 }
1629
1630
1631
1632 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
1633 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_add(matrix<bigint> & RES,long art) const1634 snf_add(matrix< bigint > &RES, long art) const
1635 {
1636 debug_handler_l(DMESSAGE, "in member - function "
1637 "snf_add(long)", DVALUE + 8);
1638
1639 register lidia_size_t i, j, z, index, SW;
1640 bigint TMP1, TMP2;
1641 bigint *tmp = NULL;
1642
1643 register lidia_size_t startr, startc, xpivot, ypivot, TEILBARKEIT;
1644 bigint ROW, COLUMN, PIVOT, NORM;
1645
1646 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1647
1648 // pivotselection: minimale C * R norm
1649 xpivot = -1;
1650 ypivot = -1;
1651 PIVOT.assign_zero();
1652 for (i = startr; i < RES.rows; i++)
1653 for (j = startc; j < RES.columns; j++) {
1654 if (PIVOT == abs(RES.value[i][j])) {
1655 RES.row_norm(ROW, i, art);
1656 RES.column_norm(COLUMN, j, art);
1657 LiDIA::add(TMP1, ROW, COLUMN);
1658 if (TMP1 < NORM) {
1659 NORM.assign(TMP1);
1660 PIVOT.assign(abs(RES.value[i][j]));
1661 xpivot = i;
1662 ypivot = j;
1663 }
1664 }
1665
1666 if (((abs_compare(PIVOT, RES.value[i][j]) > 0) && !RES.value[i][j].is_zero()) || PIVOT.is_zero()) {
1667 PIVOT.assign(abs(RES.value[i][j]));
1668 RES.row_norm(ROW, i, art);
1669 RES.column_norm(COLUMN, j, art);
1670 LiDIA::add(NORM, ROW, COLUMN);
1671 xpivot = i;
1672 ypivot = j;
1673 }
1674 }
1675
1676 if (!PIVOT.is_zero()) {
1677
1678 // swap to actual row
1679 RES.swap_rows(startr, xpivot);
1680
1681 index = ypivot;
1682
1683 TEILBARKEIT = 0;
1684
1685 while (TEILBARKEIT == 0) {
1686 TEILBARKEIT = 1;
1687
1688 // gcd2(row(startr), columns, TR2);
1689 tmp = RES.value[startr];
1690 do
1691 {
1692 SW = 0;
1693 for (i = 0; i < RES.columns; i++)
1694 if ((i != index) && !tmp[i].is_zero()) {
1695 SW = 1;
1696 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
1697 for (j = 0; j < RES.rows; j++) {
1698 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
1699 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
1700 }
1701 }
1702
1703 for (i = 0; i < RES.columns; i++)
1704 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
1705 index = i;
1706 }
1707 while (SW == 1);
1708
1709 for (i = 0; RES.value[startr][i].is_zero(); i++);
1710 RES.swap_columns(startc, i);
1711
1712 // mgcd2(column(startc), rows, TR1);
1713 index = startr;
1714 do
1715 {
1716 SW = 0;
1717 for (i = 0; i < RES.rows; i++)
1718 if ((abs_compare(RES.value[index][startc], RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
1719 index = i;
1720
1721 for (i = 0; i < RES.rows; i++)
1722 if ((i != index) && !RES.value[i][startc].is_zero()) {
1723 SW = 1;
1724 tmp = RES.value[i];
1725 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
1726 for (j = 0; j < RES.columns; j++) {
1727 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
1728 LiDIA::subtract(tmp[j], tmp[j], TMP2);
1729 }
1730 }
1731 }
1732 while (SW == 1);
1733
1734 for (i = 0; RES.value[i][startc].is_zero(); i++);
1735 RES.swap_rows(startr, i);
1736
1737 tmp = RES.value[startr];
1738 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
1739 if (index != RES.columns)
1740 TEILBARKEIT = 0;
1741
1742 index = startr;
1743 }
1744
1745 // modulo test
1746 for (i = startr; i < RES.rows; i++)
1747 for (j = startc + 1; j < RES.columns; j++)
1748 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
1749 if (i != startr)
1750 for (z = 0; z < RES.columns; z++)
1751 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
1752 i = RES.rows;
1753 j = RES.columns;
1754 startc--;
1755 startr--;
1756 }
1757 }
1758 }
1759
1760 // diagonal >= 0
1761 for (i = 0; i < RES.rows && i < RES.columns; i++)
1762 if (RES.value[i][i] < 0)
1763 RES.value[i][i].negate();
1764 }
1765
1766
1767
1768 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
1769 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_add(matrix<bigint> & RES,matrix<bigint> & T1,matrix<bigint> & T2,long art) const1770 snf_add(matrix< bigint > &RES, matrix< bigint > & T1, matrix< bigint > & T2, long art) const
1771 {
1772 debug_handler_l(DMESSAGE, "in member - function "
1773 "snf_add(matrix< bigint > &, matrix< bigint > &, long)", DVALUE + 8);
1774
1775 register lidia_size_t i, j, z, index, SW;
1776 bigint TMP1, TMP2;
1777 bigint *tmp = NULL;
1778
1779 register lidia_size_t startr, startc, xpivot, ypivot, TEILBARKEIT;
1780 bigint ROW, COLUMN, PIVOT, NORM;
1781
1782 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1783
1784 // pivotselection: minimale C * R norm
1785 xpivot = -1;
1786 ypivot = -1;
1787 PIVOT.assign_zero();
1788 for (i = startr; i < RES.rows; i++)
1789 for (j = startc; j < RES.columns; j++) {
1790 if (PIVOT == abs(RES.value[i][j])) {
1791 RES.row_norm(ROW, i, art);
1792 RES.column_norm(COLUMN, j, art);
1793 LiDIA::add(TMP1, ROW, COLUMN);
1794 if (TMP1 < NORM) {
1795 NORM.assign(TMP1);
1796 PIVOT.assign(abs(RES.value[i][j]));
1797 xpivot = i;
1798 ypivot = j;
1799 }
1800 }
1801
1802 if (((abs_compare(PIVOT, RES.value[i][j]) > 0) && !RES.value[i][j].is_zero()) || PIVOT.is_zero()) {
1803 PIVOT.assign(abs(RES.value[i][j]));
1804 RES.row_norm(ROW, i, art);
1805 RES.column_norm(COLUMN, j, art);
1806 LiDIA::add(NORM, ROW, COLUMN);
1807 xpivot = i;
1808 ypivot = j;
1809 }
1810 }
1811
1812 if (!PIVOT.is_zero()) {
1813
1814 // swap to actual row
1815 RES.swap_rows(startr, xpivot);
1816 T1.swap_rows(startr, xpivot);
1817
1818 index = ypivot;
1819
1820 TEILBARKEIT = 0;
1821
1822 while (TEILBARKEIT == 0) {
1823 TEILBARKEIT = 1;
1824
1825 // gcd2(row(startr), columns, TR2);
1826 tmp = RES.value[startr];
1827 do
1828 {
1829 SW = 0;
1830 for (i = 0; i < RES.columns; i++)
1831 if ((i != index) && !tmp[i].is_zero()) {
1832 SW = 1;
1833 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
1834 for (j = 0; j < RES.rows; j++) {
1835 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
1836 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
1837 }
1838 for (j = 0; j < RES.columns; j++) {
1839 LiDIA::multiply(TMP2, T2.value[j][index], TMP1);
1840 LiDIA::subtract(T2.value[j][i], T2.value[j][i], TMP2);
1841 }
1842 }
1843
1844 for (i = 0; i < RES.columns; i++)
1845 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
1846 index = i;
1847 }
1848 while (SW == 1);
1849
1850 for (i = 0; RES.value[startr][i].is_zero(); i++);
1851 RES.swap_columns(startc, i);
1852 T2.swap_columns(startc, i);
1853
1854 // mgcd2(column(startc), rows, TR1);
1855 index = startr;
1856 do
1857 {
1858 SW = 0;
1859 for (i = 0; i < RES.rows; i++)
1860 if ((abs_compare(RES.value[index][startc], RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
1861 index = i;
1862
1863 for (i = 0; i < RES.rows; i++)
1864 if ((i != index) && !RES.value[i][startc].is_zero()) {
1865 SW = 1;
1866 tmp = RES.value[i];
1867 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
1868 for (j = 0; j < RES.columns; j++) {
1869 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
1870 LiDIA::subtract(tmp[j], tmp[j], TMP2);
1871 }
1872 for (j = 0; j < RES.rows; j++) {
1873 LiDIA::multiply(TMP2, T1.value[index][j], TMP1);
1874 LiDIA::subtract(T1.value[i][j], T1.value[i][j], TMP2);
1875 }
1876 }
1877 }
1878 while (SW == 1);
1879
1880 for (i = 0; RES.value[i][startc].is_zero(); i++);
1881 RES.swap_rows(startr, i);
1882 T1.swap_rows(startr, i);
1883
1884 tmp = RES.value[startr];
1885 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
1886 if (index != RES.columns)
1887 TEILBARKEIT = 0;
1888
1889 index = startr;
1890 }
1891
1892 // modulo test
1893 for (i = startr; i < RES.rows; i++)
1894 for (j = startc + 1; j < RES.columns; j++)
1895 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
1896 if (i != startr) {
1897 for (z = 0; z < RES.columns; z++)
1898 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
1899 for (z = 0; z < RES.rows; z++)
1900 LiDIA::add(T1.value[startr][z], T1.value[startr][z], T1.value[i][z]);
1901 }
1902 i = RES.rows;
1903 j = RES.columns;
1904 startc--;
1905 startr--;
1906 }
1907 }
1908 }
1909
1910 // diagonal >= 0
1911 for (i = 0; i < RES.rows && i < RES.columns; i++)
1912 if (RES.value[i][i] < 0) {
1913 RES.value[i][i].negate();
1914 for (z = 0; z < RES.columns; z++)
1915 T2.value[z][i].negate();
1916 }
1917 }
1918
1919
1920
1921 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
1922 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_new(matrix<bigint> & RES,long art) const1923 snf_new(matrix< bigint > &RES, long art) const
1924 {
1925 debug_handler_l(DMESSAGE, "in member - function "
1926 "snf_new(long)", DVALUE + 8);
1927
1928 register lidia_size_t i, j, z, index, SW;
1929 bigint TMP1, TMP2;
1930 bigint *tmp = NULL;
1931
1932 register lidia_size_t startr, startc, xpivot, ypivot, TEILBARKEIT;
1933 bigint ROW, COLUMN, PIVOT, NORM;
1934
1935 bigint *RO = new bigint[RES.rows];
1936 bigint *CO = new bigint[RES.columns];
1937
1938 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
1939 // norm computation
1940 for (i = 0; i < RES.rows; i++)
1941 RES.row_norm(RO[i], i, art);
1942 for (i = 0; i < RES.columns; i++)
1943 RES.column_norm(CO[i], i, art);
1944
1945 // pivotselection: new minimale C * R norm
1946 xpivot = -1;
1947 ypivot = -1;
1948 PIVOT.assign_zero();
1949 NORM.assign_zero();
1950 for (i = startr; i < RES.rows; i++)
1951 for (j = startc; j < RES.columns; j++) {
1952 LiDIA::multiply(TMP1, RO[i], CO[j]);
1953
1954 if (!RES.value[i][j].is_zero() && (NORM > TMP1 || PIVOT.is_zero())) {
1955 NORM.assign(TMP1);
1956 PIVOT.assign(abs(RES.value[i][j]));
1957 xpivot = i;
1958 ypivot = j;
1959 }
1960
1961 if (NORM == TMP1 && !PIVOT.is_zero() && !RES.value[i][j].is_zero()) {
1962 if (abs_compare(PIVOT, RES.value[i][j]) > 0) {
1963 PIVOT.assign(abs(RES.value[i][j]));
1964 NORM.assign(TMP1);
1965 xpivot = i;
1966 ypivot = j;
1967 }
1968 }
1969 }
1970
1971 if (!PIVOT.is_zero()) {
1972
1973 // swap to actual row
1974 RES.swap_rows(startr, xpivot);
1975
1976 index = ypivot;
1977
1978 TEILBARKEIT = 0;
1979
1980 while (TEILBARKEIT == 0) {
1981 TEILBARKEIT = 1;
1982
1983 // gcd2(row(startr), columns, TR2);
1984 tmp = RES.value[startr];
1985 do
1986 {
1987 SW = 0;
1988 for (i = 0; i < RES.columns; i++)
1989 if ((i != index) && !tmp[i].is_zero()) {
1990 SW = 1;
1991 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
1992 for (j = 0; j < RES.rows; j++) {
1993 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
1994 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
1995 }
1996 }
1997
1998 for (i = 0; i < RES.columns; i++)
1999 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
2000 index = i;
2001 }
2002 while (SW == 1);
2003
2004 for (i = 0; RES.value[startr][i].is_zero(); i++);
2005 RES.swap_columns(startc, i);
2006
2007 // mgcd2(column(startc), rows, TR1);
2008 index = startr;
2009 do
2010 {
2011 SW = 0;
2012 for (i = 0; i < RES.rows; i++)
2013 if ((abs_compare(RES.value[index][startc] , RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
2014 index = i;
2015
2016 for (i = 0; i < RES.rows; i++)
2017 if ((i != index) && !RES.value[i][startc].is_zero()) {
2018 SW = 1;
2019 tmp = RES.value[i];
2020 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
2021 for (j = 0; j < RES.columns; j++) {
2022 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
2023 LiDIA::subtract(tmp[j], tmp[j], TMP2);
2024 }
2025 }
2026 }
2027 while (SW == 1);
2028
2029 for (i = 0; RES.value[i][startc].is_zero(); i++);
2030 RES.swap_rows(startr, i);
2031
2032 tmp = RES.value[startr];
2033 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
2034 if (index != RES.columns)
2035 TEILBARKEIT = 0;
2036
2037 index = startr;
2038 }
2039
2040 // modulo test
2041 for (i = startr; i < RES.rows; i++)
2042 for (j = startc + 1; j < RES.columns; j++)
2043 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
2044 if (i != startr)
2045 for (z = 0; z < RES.columns; z++)
2046 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
2047 i = RES.rows;
2048 j = RES.columns;
2049 startc--;
2050 startr--;
2051 }
2052 }
2053 }
2054
2055 // diagonal >= 0
2056 for (i = 0; i < RES.rows && i < RES.columns; i++)
2057 if (RES.value[i][i] < 0)
2058 RES.value[i][i].negate();
2059 }
2060
2061
2062
2063 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
2064 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snf_new(matrix<bigint> & RES,matrix<bigint> & T1,matrix<bigint> & T2,long art) const2065 snf_new(matrix< bigint > &RES, matrix< bigint > & T1, matrix< bigint > & T2, long art) const
2066 {
2067 debug_handler_l(DMESSAGE, "in member - function "
2068 "snf_new(matrix< bigint > &, matrix< bigint > &, long)", DVALUE + 8);
2069 register lidia_size_t i, j, z, index, SW;
2070 bigint TMP1, TMP2;
2071 bigint *tmp = NULL;
2072
2073 register lidia_size_t startr, startc, xpivot, ypivot, TEILBARKEIT;
2074 bigint ROW, COLUMN, PIVOT, NORM;
2075
2076 bigint *RO = new bigint[RES.rows];
2077 bigint *CO = new bigint[RES.columns];
2078
2079 for (startc = 0, startr = 0; startr < RES.rows && startc < RES.columns; startr++, startc++) {
2080 // norm computation
2081 for (i = 0; i < RES.rows; i++)
2082 RES.row_norm(RO[i], i, art);
2083 for (i = 0; i < RES.columns; i++)
2084 RES.column_norm(CO[i], i, art);
2085
2086 // pivotselection: new minimale C * R norm
2087 xpivot = -1;
2088 ypivot = -1;
2089 PIVOT.assign_zero();
2090 NORM.assign_zero();
2091 for (i = startr; i < RES.rows; i++)
2092 for (j = startc; j < RES.columns; j++) {
2093 //row_norm(ROW, i, art);
2094 //column_norm(RES, COLUMN, j, art);
2095 LiDIA::multiply(TMP1, RO[i], CO[j]);
2096
2097 if (!RES.value[i][j].is_zero() && (NORM > TMP1 || PIVOT.is_zero())) {
2098 NORM.assign(TMP1);
2099 PIVOT.assign(abs(RES.value[i][j]));
2100 xpivot = i;
2101 ypivot = j;
2102 }
2103
2104 if (NORM == TMP1 && !PIVOT.is_zero() && !RES.value[i][j].is_zero()) {
2105 if (abs_compare(PIVOT, RES.value[i][j]) > 0) {
2106 PIVOT.assign(abs(RES.value[i][j]));
2107 NORM.assign(TMP1);
2108 xpivot = i;
2109 ypivot = j;
2110 }
2111 }
2112 }
2113
2114 if (!PIVOT.is_zero()) {
2115
2116 // swap to actual row
2117 RES.swap_rows(startr, xpivot);
2118 T1.swap_rows(startr, xpivot);
2119
2120 index = ypivot;
2121
2122 TEILBARKEIT = 0;
2123
2124 while (TEILBARKEIT == 0) {
2125 TEILBARKEIT = 1;
2126
2127 // gcd2(row(startr), columns, TR2);
2128 tmp = RES.value[startr];
2129 do
2130 {
2131 SW = 0;
2132 for (i = 0; i < RES.columns; i++)
2133 if ((i != index) && !tmp[i].is_zero()) {
2134 SW = 1;
2135 div_rem(TMP1, TMP2, tmp[i], tmp[index]);
2136 for (j = 0; j < RES.rows; j++) {
2137 LiDIA::multiply(TMP2, RES.value[j][index], TMP1);
2138 LiDIA::subtract(RES.value[j][i], RES.value[j][i], TMP2);
2139 }
2140 for (j = 0; j < RES.columns; j++) {
2141 LiDIA::multiply(TMP2, T2.value[j][index], TMP1);
2142 LiDIA::subtract(T2.value[j][i], T2.value[j][i], TMP2);
2143 }
2144 }
2145
2146 for (i = 0; i < RES.columns; i++)
2147 if ((abs_compare(tmp[index], tmp[i]) > 0) && !tmp[i].is_zero())
2148 index = i;
2149 }
2150 while (SW == 1);
2151
2152 for (i = 0; RES.value[startr][i].is_zero(); i++);
2153 RES.swap_columns(startc, i);
2154 T2.swap_columns(startc, i);
2155
2156 // mgcd2(column(startc), rows, TR1);
2157 index = startr;
2158 do
2159 {
2160 SW = 0;
2161 for (i = 0; i < RES.rows; i++)
2162 if ((abs_compare(RES.value[index][startc], RES.value[i][startc]) > 0) && !RES.value[i][startc].is_zero())
2163 index = i;
2164
2165 for (i = 0; i < RES.rows; i++)
2166 if ((i != index) && !RES.value[i][startc].is_zero()) {
2167 SW = 1;
2168 tmp = RES.value[i];
2169 div_rem(TMP1, TMP2, tmp[startc], RES.value[index][startc]);
2170 for (j = 0; j < RES.columns; j++) {
2171 LiDIA::multiply(TMP2, RES.value[index][j], TMP1);
2172 LiDIA::subtract(tmp[j], tmp[j], TMP2);
2173 }
2174 for (j = 0; j < RES.rows; j++) {
2175 LiDIA::multiply(TMP2, T1.value[index][j], TMP1);
2176 LiDIA::subtract(T1.value[i][j], T1.value[i][j], TMP2);
2177 }
2178 }
2179 }
2180 while (SW == 1);
2181
2182 for (i = 0; RES.value[i][startc].is_zero(); i++);
2183 RES.swap_rows(startr, i);
2184 T1.swap_rows(startr, i);
2185
2186 tmp = RES.value[startr];
2187 for (index = startc+1; index < RES.columns && tmp[index].is_zero(); index++);
2188 if (index != RES.columns)
2189 TEILBARKEIT = 0;
2190
2191 index = startr;
2192 }
2193
2194 // modulo test
2195 for (i = startr; i < RES.rows; i++)
2196 for (j = startc + 1; j < RES.columns; j++)
2197 if (RES.value[i][j] % RES.value[startr][startc] != 0) {
2198 if (i != startr) {
2199 for (z = 0; z < RES.columns; z++)
2200 LiDIA::add(tmp[z], tmp[z], RES.value[i][z]);
2201 for (z = 0; z < RES.rows; z++)
2202 LiDIA::add(T1.value[startr][z], T1.value[startr][z], T1.value[i][z]);
2203 }
2204 i = RES.rows;
2205 j = RES.columns;
2206 startc--;
2207 startr--;
2208 }
2209 }
2210 }
2211
2212 // diagonal >= 0
2213 for (i = 0; i < RES.rows && i < RES.columns; i++)
2214 if (RES.value[i][i] < 0) {
2215 RES.value[i][i].negate();
2216 for (z = 0; z < RES.columns; z++)
2217 T2.value[z][i].negate();
2218 }
2219 }
2220
2221
2222
2223 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
2224 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snfmod_dkt(matrix<bigint> & RES,const bigint & mod) const2225 snfmod_dkt(matrix< bigint > &RES, const bigint &mod) const
2226 {
2227 debug_handler_l(DMESSAGE, "in member - function "
2228 "snfmod_dkt(const bigint &)", DVALUE + 8);
2229
2230 register lidia_size_t diagindex, j, z, l;
2231
2232 bigint RES0, RES1, RES2, RES3; // 0 = lggT, 1 = rggt, 2 = ggt
2233 bigint x, y;
2234 bigint TMP, TMP1, TMP2, TMP3;
2235 bigint *Atmp, *Atmp1 = NULL;
2236
2237 // A = A % mod
2238 for (z = 0; z < RES.rows; z++) {
2239 Atmp = RES.value[z];
2240 for (j = 0; j < RES.columns; j++)
2241 best_remainder(Atmp[j], Atmp[j], mod);
2242 }
2243
2244 // Step 3 - 5
2245 for (diagindex = 0; diagindex < RES.rows; diagindex++) {
2246 Atmp = RES.value[diagindex];
2247
2248 // if diagonalelement == 0 then diagonalelement = mod
2249 if (Atmp[diagindex].is_zero())
2250 Atmp[diagindex].assign(mod);
2251
2252 // Step 6 -8
2253 for (j = diagindex+1; j < RES.columns; j++)
2254 if (!Atmp[j].is_zero()) {
2255 // columns reduction
2256 RES2 = xgcd(RES0, RES1, Atmp[j], Atmp[diagindex]);
2257 div_rem(x, RES3, Atmp[diagindex], RES2);
2258 div_rem(y, RES3, Atmp[j], RES2);
2259
2260 for (z = 0; z < RES.rows; z++) {
2261 Atmp1 = RES.value[z];
2262 TMP.assign(Atmp1[j]);
2263 TMP1.assign(Atmp1[diagindex]);
2264
2265 mult_mod(TMP2, TMP, x, mod);
2266 mult_mod(TMP3, TMP1, y, mod);
2267 sub_mod(Atmp1[j], TMP2, TMP3, mod);
2268
2269 mult_mod(TMP2, TMP, RES0, mod);
2270 mult_mod(TMP3, TMP1, RES1, mod);
2271 add_mod(Atmp1[diagindex], TMP2, TMP3, mod);
2272 }
2273 }
2274
2275 for (j = diagindex+1; j < RES.rows; j++)
2276 if (!RES.value[j][diagindex].is_zero()) {
2277 // row reduction
2278 RES2 = xgcd(RES0, RES1, RES.value[j][diagindex], Atmp[diagindex]);
2279 div_rem(x, RES3, Atmp[diagindex], RES2);
2280 div_rem(y, RES3, RES.value[j][diagindex], RES2);
2281
2282 for (z = 0; z < RES.columns; z++) {
2283 TMP.assign(RES.value[j][z]);
2284 TMP1.assign(RES.value[diagindex][z]);
2285
2286 mult_mod(TMP2, TMP, x, mod);
2287 mult_mod(TMP3, TMP1, y, mod);
2288 sub_mod(RES.value[j][z], TMP2, TMP3, mod);
2289
2290 mult_mod(TMP2, TMP, RES0, mod);
2291 mult_mod(TMP3, TMP1, RES1, mod);
2292 add_mod(RES.value[diagindex][z], TMP2, TMP3, mod);
2293 }
2294 }
2295
2296 // value[diagindex][diagindex] | value[i][j] ???
2297 TMP = Atmp[diagindex];
2298 for (j = diagindex+1; j < RES.rows; j++)
2299 for (z = diagindex+1; z < RES.columns; z++) {
2300 if (RES.value[j][z] % TMP != 0) {
2301 if (j != diagindex)
2302 for (l = diagindex; l < RES.columns; l++)
2303 add_mod(Atmp[l], Atmp[l], RES.value[j][l], mod);
2304 j = RES.rows;
2305 z = RES.columns;
2306 }
2307 }
2308
2309 for (z = diagindex+1; z < RES.columns && Atmp[z].is_zero(); z++);
2310 if (z != RES.columns)
2311 diagindex--;
2312 }
2313
2314 // Step 29 - 32
2315 bigint D = mod;
2316 for (j = 0; j < RES.rows; j++) {
2317 Atmp = RES.value[j];
2318 Atmp[j].assign(gcd(Atmp[j], D));
2319 div_rem(D, TMP, D, Atmp[j]);
2320 }
2321 }
2322
2323
2324
2325 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
2326 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
snfmod_cohen(matrix<bigint> & RES,const bigint & mod) const2327 snfmod_cohen(matrix< bigint > &RES, const bigint & mod) const
2328 {
2329 debug_handler_l(DMESSAGE, "in member - function "
2330 "snfmod_cohen(const bigint &)", DVALUE + 8);
2331
2332 register lidia_size_t diagindex, j, z, l;
2333
2334 bigint RES0, RES1, RES2, RES3; // 0 = lggT, 1 = rggt, 2 = ggt
2335 bigint x, y;
2336 bigint TMP, TMP1, TMP2, TMP3;
2337 bigint *Atmp, *Atmp1 = NULL;
2338 bigint D = mod;
2339
2340 // Step 1
2341 for (diagindex = 0; diagindex < RES.rows; diagindex++) {
2342 Atmp = RES.value[diagindex];
2343
2344 if (Atmp[diagindex].is_zero())
2345 Atmp[diagindex].assign(mod);
2346
2347 // Step 2 - 4
2348 for (j = diagindex+1; j < RES.columns; j++)
2349 if (!Atmp[j].is_zero()) {
2350 // columns reduction
2351 RES2 = xgcd(RES0, RES1, Atmp[j], Atmp[diagindex]);
2352 div_rem(x, RES3, Atmp[diagindex], RES2);
2353 div_rem(y, RES3, Atmp[j], RES2);
2354
2355 for (z = 0; z < RES.rows; z++) {
2356 Atmp1 = RES.value[z];
2357 TMP.assign(Atmp1[j]);
2358 TMP1.assign(Atmp1[diagindex]);
2359
2360 mult_mod(TMP2, TMP, x, mod);
2361 mult_mod(TMP3, TMP1, y, mod);
2362 sub_mod(Atmp1[j], TMP2, TMP3, mod);
2363
2364 mult_mod(TMP2, TMP, RES0, mod);
2365 mult_mod(TMP3, TMP1, RES1, mod);
2366 add_mod(Atmp1[diagindex], TMP2, TMP3, mod);
2367 }
2368 }
2369
2370 // Step 5 - 7
2371 for (j = diagindex+1; j < RES.rows; j++)
2372 if (!RES.value[j][diagindex].is_zero()) {
2373 // row reduction
2374 RES2 = xgcd(RES0, RES1, RES.value[j][diagindex], Atmp[diagindex]);
2375 div_rem(x, RES3, Atmp[diagindex], RES2);
2376 div_rem(y, RES3, RES.value[j][diagindex], RES2);
2377
2378 for (z = 0; z < RES.columns; z++) {
2379 TMP.assign(RES.value[j][z]);
2380 TMP1.assign(RES.value[diagindex][z]);
2381
2382 mult_mod(TMP2, TMP, x, mod);
2383 mult_mod(TMP3, TMP1, y, mod);
2384 sub_mod(RES.value[j][z], TMP2, TMP3, mod);
2385
2386 mult_mod(TMP2, TMP, RES0, mod);
2387 mult_mod(TMP3, TMP1, RES1, mod);
2388 add_mod(RES.value[diagindex][z], TMP2, TMP3, mod);
2389 }
2390 }
2391
2392 // Step 8, 9
2393 TMP = Atmp[diagindex];
2394
2395 for (j = diagindex+1; j < RES.rows; j++)
2396 for (z = diagindex+1; z < RES.columns; z++) {
2397 if (RES.value[j][z] % TMP != 0) {
2398 if (j != diagindex)
2399 for (l = diagindex; l < RES.columns; l++)
2400 add_mod(Atmp[l], Atmp[l], RES.value[j][l], mod);
2401 j = RES.rows;
2402 z = RES.columns;
2403 }
2404 }
2405
2406 for (z = diagindex+1; z < RES.columns && Atmp[z].is_zero(); z++);
2407 if (z != RES.columns)
2408 diagindex--;
2409 else {
2410 // Step 10
2411 Atmp[diagindex] = xgcd(TMP1, TMP2, TMP, D);
2412 div_rem(D, TMP1, D, Atmp[diagindex]);
2413 }
2414 }
2415 }
2416
2417
2418
2419 //
2420 // END: Linear algebra
2421 // PART 2
2422 //
2423
2424 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
2425 inline void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
gauss(matrix<bigint> & RES) const2426 gauss(matrix< bigint > &RES) const
2427 {
2428 debug_handler(DMESSAGE, "in member - function gauss()");
2429
2430 matrix< bigint > TR(RES.columns, RES.columns);
2431 bigint *REM = NULL;
2432 register lidia_size_t startr = 0, startc = 0, i;
2433
2434 for (startc = RES.columns - 1, startr = RES.rows - 1; startr >= 0 && startc >= 0; startr--, startc--) {
2435
2436 bigint *ZU = new bigint[RES.columns];
2437 RES.get_row(ZU, startr);
2438 for (i = startc + 1; i < RES.columns; ZU[i].assign_zero(), i++);
2439 REM = TR.mgcd1(ZU, RES.columns);
2440 delete[] REM;
2441 delete[] ZU;
2442
2443 TR = TR.trans();
2444 LiDIA::multiply(RES, RES, TR);
2445 for (i = 0; RES.value[startr][i].is_zero() && i <= startc; i++);
2446 if (i > startc)
2447 return;
2448 RES.swap_columns(startc, i);
2449
2450 }
2451 }
2452
2453
2454
2455 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
2456 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
basis_completion(matrix<bigint> & M,bigint * x,lidia_size_t n) const2457 basis_completion(matrix< bigint > &M, bigint *x, lidia_size_t n) const
2458 {
2459 bigint *y = new bigint[n+1];
2460
2461 lidia_size_t i, j;
2462 for (i = 1; i <= n; i++) {
2463 y[i] = x[i-1];
2464 }
2465 y[0] = -1;
2466
2467 if (M.rows != n)
2468 M.set_no_of_rows(n);
2469 if (M.columns != n+1)
2470 M.set_no_of_columns(n+1);
2471
2472 // diag
2473 for (i = 0; i < n; i++) {
2474 M.value[i][i+1].assign_one();
2475 M.value[i][0].assign(x[i]);
2476 }
2477
2478 bigint a, b, g, c, d, gold = y[1], TMP;
2479
2480 for (i = 1; i < n - 1; i++) {
2481 if (y[i] == 1) {
2482 for (j = i; j < n; j++)
2483 M.swap_columns(j, j+1);
2484 M.columns--;
2485 return;
2486 }
2487
2488 if (y[i + 1] == 0) {
2489 continue;
2490 }
2491 g = xgcd(a, b, y[i + 1], gold);
2492 b.negate();
2493 c = gold / g;
2494 d = y[i+1] /g;
2495
2496 for (j = 0; j < n; j++) {
2497
2498 TMP = M.value[j][i];
2499 M.value[j][i].assign(TMP * a + M.value[j][i+1] * b);
2500 M.value[j][i+1].assign(TMP * c + M.value[j][i+1] * d);
2501 }
2502 y[i+1] = g;
2503 y[i] = 0;
2504
2505 gold = g;
2506 }
2507 }
2508
2509
2510
2511 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
2512 void modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
cond_matrix(matrix<bigint> & M,bigint * v,lidia_size_t n) const2513 cond_matrix(matrix< bigint > &M, bigint *v, lidia_size_t n) const
2514 {
2515 if (M.rows != n)
2516 M.set_no_of_rows(n);
2517 if (M.columns != n)
2518 M.set_no_of_columns(n);
2519
2520 M.diag(1, 0);
2521
2522 lidia_size_t i;
2523
2524 bigint *b = new bigint[n];
2525 for (i = 0; i < n; i++)
2526 b[i] = v[i];
2527
2528 bigint g, N = b[0], a = b[1], qa, ra, qb, rb, t = 0;
2529
2530 for (i = 2; i < n && g != 1; i++) {
2531 if (b[i] != 0) {
2532 g = gcd(b[i], a);
2533
2534 div_rem(qa, ra, a/g, N);
2535 div_rem(qb, rb, b[i]/g, N);
2536
2537 for (t = 0; gcd((ra + t*rb), N) != 1; t++);
2538
2539 a = a + t * b[i]; // NEW
2540
2541 M.sto(1, i, t);
2542 }
2543 }
2544 }
2545
2546
2547
2548 //
2549 // mgcd2
2550 //
2551
2552 template< class REP, class SINGLE_MODUL, class MULTI_MODUL >
2553 bigint *modular_bigint_matrix_algorithms< REP, SINGLE_MODUL, MULTI_MODUL >::
mgcd2(matrix<bigint> & RES,const bigint * aconst,lidia_size_t n) const2554 mgcd2(matrix< bigint > &RES, const bigint * aconst, lidia_size_t n) const
2555 {
2556 //
2557 // DESCRIPTION: RES = T.mgcd2(a, n);
2558 // = > RES[0] = RES[1]*a[0] + ... + RES[n]*a[n-1]
2559 // = > RES[0] = gcd(a[0], ..., a[n-1])
2560 // = > T*a = RES
2561 // ALGORITHM: Blankinship
2562 // IMPROVEMENTS: Havas, Majewski, reduction of all elements, MIN assignments
2563 // PAPER: Hermite normal form computation for integer matrices, Havas
2564 // VERSION: 1.8
2565 //
2566
2567 debug_handler("multiple_gcd", "in member - function "
2568 "mgcd2(const bigint *, lidia_size_t)");
2569
2570 register lidia_size_t i, j, index, bound, SW;
2571 bigint MIN, TMP, q, r, *Ttmp1, *Ttmp2 = NULL;
2572
2573 if (RES.columns != n)
2574 RES.set_no_of_columns(n);
2575 if (RES.rows != n)
2576 RES.set_no_of_rows(n);
2577 RES.diag(1, 0);
2578
2579 bigint *a = new bigint[n + 1];
2580 memory_handler(a, "multiple_gcd", "mgcd2 :: "
2581 "Error in memory allocation (a)");
2582
2583 for (i = 0; i < n; i++)
2584 a[i].assign(aconst[i]);
2585
2586 // init
2587 for (index = 0; index < n && a[index].is_zero(); index++);
2588
2589 if (index == n) {
2590 delete[] a;
2591 return new bigint[n];
2592 }
2593 else
2594 bound = index;
2595
2596 do
2597 {
2598 MIN.assign(a[index]);
2599
2600 // Pivot search: MINIMUM
2601 for (i = bound; i < n; i++)
2602 if ((abs_compare(MIN, a[i]) > 0) && !a[i].is_zero()) {
2603 MIN.assign(a[i]);
2604 index = i;
2605 }
2606
2607 // all elements
2608 SW = 0;
2609 Ttmp2 = RES.value[index];
2610 for (i = bound; i < n; i++)
2611 if ((i != index) && !a[i].is_zero()) {
2612 SW = 1;
2613 Ttmp1 = RES.value[i];
2614 div_rem(q, r, a[i], MIN);
2615 a[i].assign(r);
2616 for (j = 0; j < n; j++) {
2617 LiDIA::multiply(TMP, q, Ttmp2[j]);
2618 LiDIA::subtract(Ttmp1[j], Ttmp1[j], TMP);
2619 }
2620 }
2621 }
2622 while (SW == 1);
2623
2624 Ttmp2 = RES.value[index];
2625
2626 // gcd < 0 ?
2627 if (a[index] < 0) {
2628 a[index].negate();
2629 for (i = 0; i < n; i++)
2630 Ttmp2[i].negate();
2631 }
2632
2633 if (index != 0)
2634 a[0].assign(a[index]);
2635 for (i = 1; i <= n; i++)
2636 a[i].assign(Ttmp2[i - 1]);
2637
2638 return a;
2639 }
2640
2641
2642
2643 #undef DVALUE
2644 #undef DMESSAGE
2645 #undef EMESSAGE
2646
2647
2648
2649 #ifdef LIDIA_NAMESPACE
2650 # ifndef IN_NAMESPACE_LIDIA
2651 } // end of namespace LiDIA
2652 # endif
2653 #endif
2654
2655
2656
2657 #endif // LIDIA_BIGINT_MATRIX_ALGORITHMS_CC_GUARD_
2658