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 : Thorsten Lauer (TL)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22 #include "LiDIA/bigfloat_matrix.h"
23 #include "LiDIA/bigfloat_int.h"
24
25
26
27 #ifdef LIDIA_NAMESPACE
28 namespace LiDIA {
29 #endif
30
31
32
33 void
34 bigfloat_matrix::
gauss_jordan(base_vector<bigfloat> & x,const base_vector<bigfloat> & b) const35 gauss_jordan(base_vector< bigfloat > & x, const base_vector< bigfloat > & b) const
36 {
37 debug_handler("bigfloat_matrix", "gauss_jordan_pc(x, b)");
38 if (rows != b.size())
39 lidia_error_handler("bigfloat_lattice", "gauss_jordan_pc(x, b) :: illegal vector size");
40 if (rows > columns)
41 lidia_error_handler("bigfloat_lattice", "gauss_jordan_pc(x, b) :: "
42 "rows have to be less or equal columns");
43
44 long i, j, k, r = 0, c = 0, *index = new long[columns]; // MM
45 bigfloat_int factor, tmp, *help0, *help1, *helpr, *l, **p;
46 bigfloat max, temp;
47
48 x.set_capacity(rows);
49
50 help1 = new bigfloat_int[columns];
51 memory_handler(help1, "bigfloat_lattice", "gauss_jordan_pc(x, b) :: not enough memory !");
52
53 l = new bigfloat_int[rows];
54 memory_handler(l, "bigfloat_lattice", "gauss_jordan_pc(x, b) :: not enough memory !");
55
56 for (i = 0; i < rows; i++)
57 l[i].assign(b[i]);
58 for (i = 0; i < columns; i++)
59 index[i] = i;
60
61 p = new bigfloat_int*[rows];
62 memory_handler(p, "bigfloat_matrix", "gauss_jordan(b) :: not enough memory");
63 for (i = 0; i < rows; i++) {
64 p[i] = new bigfloat_int[columns];
65 memory_handler(p[i], "bigfloat_matrix", "gauss_jordan(b) :: not enough memory");
66 }
67
68 for (i = 0; i < rows; i++)
69 for (j = 0; j < columns; j++)
70 p[i][j].assign(value[i][j]);
71
72 // Begin
73
74 for (i = 0; i < rows; i++) {
75 // search for the biggest absolute value in column
76
77 max.assign_zero();
78 k = i;
79 for (j = i; j < rows; j++)
80 for (k = i; k < rows; k++) {
81 p[j][k].approx(temp);
82 temp = abs(temp);
83 if (max < temp) {
84 r = j;
85 c = k;
86 max.assign(temp);
87 }
88 }
89 if (r != i) {
90 helpr = p[i];
91 p[i] = p[r];
92 p[r] = helpr;
93 LiDIA::swap(l[i], l[r]);
94 }
95 if (c != i) {
96 for (j = 0; j < rows; j++) {
97 factor.assign(p[j][c]);
98 p[j][c].assign(p[j][i]);
99 p[j][i].assign(factor);
100 }
101 k = index[i];
102 index[i] = index[c];
103 index[c] = k;
104 }
105
106 // is it zero ? then finish
107
108 p[i][i].approx(temp);
109 temp = abs(temp);
110 if (temp.is_approx_zero())
111 break;
112
113 // else continue
114
115 helpr = p[i];
116 for (j = i+1; j < rows; j++) {
117 LiDIA::divide(factor, p[j][i], p[i][i]);
118 help0 = p[j];
119 for (k = i; k < columns; k++) {
120 LiDIA::multiply(help1[k], helpr[k], factor);
121 LiDIA::subtract(help0[k], help0[k], help1[k]);
122 }
123 LiDIA::multiply(tmp, l[i], factor);
124 LiDIA::subtract(l[j], l[j], tmp);
125 }
126 for (j = 0; j < i; j++) {
127 LiDIA::divide(factor, p[j][i], p[i][i]);
128 help0 = p[j];
129 for (k = i; k < columns; k++) {
130 LiDIA::multiply(help1[k], helpr[k], factor);
131 LiDIA::subtract(help0[k], help0[k], help1[k]);
132 }
133 LiDIA::multiply(tmp, l[i], factor);
134 LiDIA::subtract(l[j], l[j], tmp);
135 }
136 }
137
138 // check if the equation-system is solvable
139
140 for (j = i; j < rows; j++) {
141 l[j].approx(temp);
142 if (!temp.is_approx_zero()) {
143 for (i = 0; i < rows; i++)
144 x[i].assign_zero();
145 delete[] index; // MM
146 return; // no solution, return zero-vector
147 }
148 }
149
150 // else compute x in Ax=b
151
152 for (j = 0; j < i; j++) {
153 LiDIA::divide(factor, l[j], p[j][j]);
154 factor.approx(x[index[j]]);
155 }
156
157 for (j = i; j < columns; j++)
158 x[index[j]].assign_zero();
159
160 for (i = 0; i < rows; i++)
161 delete[] p[i];
162 delete[] p;
163
164 delete[] help1;
165 delete[] l;
166 delete[] index; // MM
167 }
168
169
170
171 base_vector< bigfloat >
gauss_jordan(const base_vector<bigfloat> & b) const172 bigfloat_matrix::gauss_jordan(const base_vector< bigfloat > & b) const
173 {
174 debug_handler("bigfloat_matrix", "gauss_jordan_pc(b)");
175 if (rows != b.size())
176 lidia_error_handler("bigfloat_lattice", "gauss_jordan_pc(b) :: illegal vector size");
177 if (rows > columns)
178 lidia_error_handler("bigfloat_lattice", "gauss_jordan_pc(b) :: "
179 "rows have to be less or equal columns");
180
181 base_vector< bigfloat > x(columns, vector_flags(vector_flags::fixed));
182 long i, j, k, r = 0, c = 0, *index = new long[columns];
183 bigfloat_int factor, tmp, *help0, *help1, *helpr, *l, **p;
184 bigfloat max, temp;
185
186 help1 = new bigfloat_int[columns];
187 memory_handler(help1, "bigfloat_lattice", "gauss_jordan_pc(b) :: not enough memory !");
188
189 l = new bigfloat_int[rows];
190 memory_handler(l, "bigfloat_lattice", "gauss_jordan_pc(b) :: not enough memory !");
191
192 for (i = 0; i < rows; i++)
193 l[i].assign(b[i]);
194 for (i = 0; i < columns; i++)
195 index[i] = i;
196
197 p = new bigfloat_int*[rows];
198 memory_handler(p, "bigfloat_matrix", "gauss_jordan(b) :: not enough memory");
199 for (i = 0; i < rows; i++) {
200 p[i] = new bigfloat_int[columns];
201 memory_handler(p[i], "bigfloat_matrix", "gauss_jordan(b) :: not enough memory");
202 }
203
204 for (i = 0; i < rows; i++)
205 for (j = 0; j < columns; j++)
206 p[i][j].assign(value[i][j]);
207
208 // Begin
209
210 for (i = 0; i < rows; i++) {
211 // search for the biggest absolute value in column
212
213 max.assign_zero();
214 k = i;
215 for (j = i; j < rows; j++)
216 for (k = i; k < rows; k++) {
217 p[j][k].approx(temp);
218 temp = abs(temp);
219 if (max < temp) {
220 r = j;
221 c = k;
222 max.assign(temp);
223 }
224 }
225 if (r != i) {
226 helpr = p[i];
227 p[i] = p[r];
228 p[r] = helpr;
229 LiDIA::swap(l[i], l[r]);
230 }
231 if (c != i) {
232 for (j = 0; j < rows; j++) {
233 factor.assign(p[j][c]);
234 p[j][c].assign(p[j][i]);
235 p[j][i].assign(factor);
236 }
237 k = index[i];
238 index[i] = index[c];
239 index[c] = k;
240 }
241
242 // is it zero ? then finish
243
244 p[i][i].approx(temp);
245 temp = abs(temp);
246 if (temp.is_approx_zero())
247 break;
248
249 // else continue
250
251 helpr = p[i];
252 for (j = i+1; j < rows; j++) {
253 LiDIA::divide(factor, p[j][i], p[i][i]);
254 help0 = p[j];
255 for (k = i; k < columns; k++) {
256 LiDIA::multiply(help1[k], helpr[k], factor);
257 LiDIA::subtract(help0[k], help0[k], help1[k]);
258 }
259 LiDIA::multiply(tmp, l[i], factor);
260 LiDIA::subtract(l[j], l[j], tmp);
261 }
262 for (j = 0; j < i; j++) {
263 LiDIA::divide(factor, p[j][i], p[i][i]);
264 help0 = p[j];
265 for (k = i; k < columns; k++) {
266 LiDIA::multiply(help1[k], helpr[k], factor);
267 LiDIA::subtract(help0[k], help0[k], help1[k]);
268 }
269 LiDIA::multiply(tmp, l[i], factor);
270 LiDIA::subtract(l[j], l[j], tmp);
271 }
272 }
273
274 // check if the equation-system is solvable
275
276 for (j = i; j < rows; j++) {
277 l[j].approx(temp);
278 if (!temp.is_approx_zero()) {
279 delete [] index; // MM
280 return x; // no solution, return zero-vector
281 }
282 }
283
284 // else compute x in Ax=b
285
286 for (j = 0; j < i; j++) {
287 LiDIA::divide(factor, l[j], p[j][j]);
288 factor.approx(x[index[j]]);
289 }
290
291 for (j = i; j < columns; j++)
292 x[index[j]].assign_zero();
293
294 for (i = 0; i < rows; i++)
295 delete[] p[i];
296 delete[] p;
297
298 delete[] help1;
299 delete[] l;
300 delete[] index; // MM
301 // return it
302
303 return x;
304 }
305
306
307
308 void
cholesky(bigfloat_matrix & A) const309 bigfloat_matrix::cholesky(bigfloat_matrix & A) const
310 {
311 debug_handler("bigfloat_matrix", "cholesky()");
312 if (rows != columns)
313 lidia_error_handler("bigfloat_matrix", "cholesky() :: matrix is not symmetrical");
314
315 bigfloat_int l, **p, *helpr;
316 bigfloat max, temp;
317 long n, k, i, r, c, j, *index = new long[columns];
318
319 p = new bigfloat_int*[rows];
320 memory_handler(p, "bigfloat_matrix", "cholesky() :: not enough memory");
321 for (n = 0; n < rows; n++) {
322 p[n] = new bigfloat_int[columns];
323 memory_handler(p[n], "bigfloat_matrix", "cholesky() :: not enough memory");
324 index[n] = n;
325 }
326
327 for (n = 0; n < rows; n++) {
328 max.assign_zero();
329 r = c = n;
330 for (j = n; j < rows; j++)
331 for (k = n; k < columns; k++) {
332 p[j][k].approx(temp);
333 temp = abs(temp);
334 if (max < temp) {
335 r = j;
336 c = k;
337 max.assign(temp);
338 }
339 }
340 if (c != n) {
341 helpr = p[n];
342 p[n] = p[c];
343 p[c] = helpr;
344 }
345 if (r != n) {
346 for (i = 0; i < rows; i++) {
347 l = p[i][r];
348 p[i][r] = p[i][n];
349 p[i][n] = l;
350 }
351 k = index[n];
352 index[n] = index[r];
353 index[r] = k;
354 }
355
356
357 (p[n][n]).assign(value[n][n]);
358 for (k = 0; k < n; k++) {
359 LiDIA::square(l, p[n][k]);
360 LiDIA::subtract(p[n][n], p[n][n], l);
361 }
362 LiDIA::sqrt(p[n][n], p[n][n]);
363 for (i = n+1; i < rows; i++) {
364 (p[i][n]).assign(value[i][n]);
365 for (k = 0; k < n; k++) {
366 LiDIA::multiply(l, p[i][k], p[n][k]);
367 LiDIA::subtract(p[i][n], p[i][n], l);
368 }
369 LiDIA::divide(p[i][n], p[i][n], p[n][n]);
370 }
371 }
372
373 A.resize(rows, columns);
374 for (n = 0; n < rows; n++)
375 for (k = 0; k < rows; k++)
376 p[n][index[k]].approx(A.value[n][k]);
377
378 for (i = 0; i < rows; i++)
379 delete[] p[i];
380 delete[] p;
381 delete[] index; // MM
382 }
383
384
385
386 void
mod_cholesky(bigfloat_matrix & Q) const387 bigfloat_matrix::mod_cholesky(bigfloat_matrix & Q) const
388 {
389 debug_handler("bigfloat_matrix", "mod_cholesky(Q)");
390 if (rows != columns)
391 lidia_error_handler("bigfloat_lattice", "mod_cholesky(Q) :: rows != columns");
392
393 long i, j, s, l;
394 bigfloat_int temp, **p;
395
396 Q.resize(rows, columns);
397
398 p = new bigfloat_int*[rows];
399 for (i = 0; i < rows; i++)
400 p[i] = new bigfloat_int[columns];
401
402 for (i = 0; i < rows; i++)
403 for (j = i; j < columns; j++)
404 p[i][j].assign(value[i][j]);
405
406 for (i = 0; i < rows; i++) {
407 for (j = i+1; j < rows; j++) {
408 p[j][i].assign(p[i][j]);
409 LiDIA::divide(p[i][j], p[i][j], p[i][i]);
410 }
411
412 for (s = i+1; s < rows; s++) {
413 for (l = s; l < columns; l++) {
414 LiDIA::multiply(temp, p[s][i], p[i][l]);
415 LiDIA::subtract(p[s][l], p[s][l], temp);
416 }
417 }
418 }
419
420 for (i = 0; i < rows; i++)
421 for (j = 0; j < columns; j++)
422 (p[i][j]).approx(Q.value[i][j]);
423
424 for (i = 0; i < rows; i++)
425 delete[] p[i];
426 delete[] p;
427 }
428
429
430
431 void
invert(bigfloat_matrix & A) const432 bigfloat_matrix::invert(bigfloat_matrix & A) const
433 {
434 debug_handler("bigfloat_matrix", "invert(A)");
435 if (rows != columns)
436 lidia_error_handler("bigfloat_lattice", "invert(A) :: rows != columns");
437
438 long i, j, k, c = 0, r = 0, *index = new long[columns];
439 bigfloat temp, max;
440 bigfloat_int tempi, **p, *helpr, *help0, factor;
441 bigfloat_matrix I(rows, columns), M(rows, columns+columns);
442
443 p = new bigfloat_int*[rows];
444 for (i = 0; i < rows; i++)
445 p[i] = new bigfloat_int[columns+columns];
446
447 for (i = 0; i < rows; i++) {
448 for (j = 0; j < columns; j++)
449 p[i][j].assign(value[i][j]);
450 (p[i][i+columns]).assign_one();
451 index[i] = i;
452 }
453
454 // Begin
455
456 for (i = 0; i < rows; i++) {
457 // search for the biggest absolute value
458
459 max.assign_zero();
460 k = i;
461 for (j = i; j < rows; j++)
462 for (k = i; k < rows; k++) {
463 p[j][k].approx(temp);
464 temp = abs(temp);
465 if (max < temp) {
466 r = j;
467 c = k;
468 max.assign(temp);
469 }
470 }
471 if (r != i) {
472 helpr = p[i];
473 p[i] = p[r];
474 p[r] = helpr;
475 }
476 if (c != i) {
477 for (k = 0; k < rows; k++) {
478 LiDIA::swap(p[k][i], p[k][c]);
479 }
480 k = index[i];
481 index[i] = index[c];
482 index[c] = k;
483 }
484
485 // is it zero ? then finish
486
487 p[i][i].approx(temp);
488 if (temp.is_approx_zero())
489 break;
490
491 // else continue
492
493 p[i][i].Approx();
494 helpr = p[i];
495 for (j = i+1; j < rows; j++) {
496 LiDIA::divide(factor, p[j][i], p[i][i]);
497 help0 = p[j];
498 for (k = i; k < columns+columns; k++) {
499 LiDIA::multiply(tempi, helpr[k], factor);
500 LiDIA::subtract(help0[k], help0[k], tempi);
501 }
502 }
503 for (j = 0; j < i; j++) {
504 LiDIA::divide(factor, p[j][i], p[i][i]);
505 help0 = p[j];
506 for (k = i; k < columns+columns; k++) {
507 LiDIA::multiply(tempi, helpr[k], factor);
508 LiDIA::subtract(help0[k], help0[k], tempi);
509 }
510 }
511 }
512
513 for (j = 0; j < rows; j++)
514 for (i = 0; i < columns; i++)
515 LiDIA::divide(p[j][i+columns], p[j][i+columns], p[j][j]);
516
517 for (j = 0; j < rows; j++) {
518 while (index[j] != j) {
519 helpr = p[j];
520 p[j] = p[index[j]];
521 p[index[j]] = helpr;
522 i = index[index[j]];
523 index[index[j]] = index[j];
524 index[j] = i;
525 }
526 }
527
528 A.resize(rows, columns);
529 for (i = 0; i < rows; i++)
530 for (j = 0; j < columns; j++)
531 p[i][columns+j].approx(A.value[i][j]);
532
533 for (i = 0; i < rows; i++)
534 delete[] p[i];
535 delete[] p;
536 delete[] index; // MM
537 }
538
539
540 void
LR(bigfloat_matrix & L,bigfloat_matrix & R) const541 bigfloat_matrix::LR(bigfloat_matrix & L, bigfloat_matrix & R) const
542 {
543 debug_handler("bigfloat_matrix", "LR(L, R)");
544 if (rows != columns)
545 lidia_error_handler("bigfloat_matrix", "LR(L, R) :: rows != columns");
546 lidia_size_t i, j, k;
547 bigfloat_int t, **l, **r, **a;
548 bigfloat temp;
549
550 L.resize(rows, columns);
551 R.resize(rows, columns);
552
553 l = new bigfloat_int*[rows];
554 for (i = 0; i < rows; i++)
555 l[i] = new bigfloat_int[columns];
556
557 r = new bigfloat_int*[rows];
558 for (i = 0; i < rows; i++)
559 r[i] = new bigfloat_int[columns];
560
561 a = new bigfloat_int*[rows];
562 for (i = 0; i < rows; i++)
563 a[i] = new bigfloat_int[columns];
564
565 for (i = 0; i < rows; i++)
566 for (j = 0; j < columns; j++) {
567 l[i][j].assign_zero();
568 r[i][j].assign_zero();
569 a[i][j].assign(value[i][j]);
570 }
571
572 for (j = 0; j < columns; j++) {
573 for (k = 0; k < j; k++)
574 for (i = k+1; i <= j; i++) {
575 LiDIA::multiply(t, a[i][k], a[k][j]);
576 LiDIA::subtract(a[i][j], a[i][j], t);
577 }
578
579 for (k = 0; k < j; k++)
580 for (i = j+1; i < rows; i++) {
581 LiDIA::multiply(t, a[i][k], a[k][j]);
582 LiDIA::subtract(a[i][j], a[i][j], t);
583 }
584
585 for (i = j+1; i < rows; i++)
586 LiDIA::divide(a[i][j], a[i][j], a[j][j]);
587 }
588
589 for (i = 0; i < rows; i++) {
590 for (j = 0; j < i; j++)
591 (a[i][j]).approx(L.value[i][j]);
592 for (; j < columns; j++)
593 (a[i][j]).approx(R.value[i][j]);
594 (L.value[i][i]).assign_one();
595 }
596
597 for (i = 0; i < rows; i++) {
598 delete[] l[i];
599 delete[] r[i];
600 delete[] a[i];
601 }
602 delete[] l;
603 delete[] r;
604 delete[] a;
605 }
606
607
608 void
QR(bigfloat_matrix & Q,bigfloat_matrix & R) const609 bigfloat_matrix::QR(bigfloat_matrix & Q, bigfloat_matrix & R) const
610 {
611 debug_handler("bigfloat_matrix", "QR(Q, R)");
612 if (rows != columns)
613 lidia_error_handler("bigfloat_matrix", "QR(Q, R) :: rows != columns");
614 lidia_size_t i, j, k;
615 bigfloat_int a, a2, c, **q, **r;
616 bigfloat temp;
617 base_vector< bigfloat_int > u(rows, vector_flags(vector_flags::fixed));
618 base_vector< bigfloat_int > v(rows, vector_flags(vector_flags::fixed));
619 base_vector< bigfloat_int > w(columns, vector_flags(vector_flags::fixed));
620 base_vector< bigfloat_int > x(columns, vector_flags(vector_flags::fixed));
621
622 Q.resize(rows, columns);
623 R.resize(rows, columns);
624
625 q = new bigfloat_int*[rows];
626 for (i = 0; i < rows; i++)
627 q[i] = new bigfloat_int[columns];
628
629 for (i = 0; i < rows; i++) {
630 for (j = 0; j < columns; j++)
631 q[i][j].assign_zero();
632 q[i][i].assign_one();
633 }
634
635 r = new bigfloat_int*[rows];
636 for (i = 0; i < rows; i++)
637 r[i] = new bigfloat_int[columns];
638
639 for (i = 0; i < rows; i++)
640 for (j = 0; j < columns; j++)
641 r[i][j].assign(value[i][j]);
642
643 for (j = 0; j < columns-1; j++) {
644 a2.assign_zero();
645
646 for (k = j; k < rows; k++) {
647 square(a, r[k][j]);
648 LiDIA::add(a2, a2, a);
649 }
650
651 a2.approx(temp);
652 if (temp.is_approx_zero())
653 lidia_error_handler("bigfloat_matrix", "QR(Q, R) :: Matrix was singular");
654 a = sqrt(a2);
655
656 LiDIA::multiply(c, a, abs(r[j][j]));
657 LiDIA::add(c, c, a2);
658 c.invert();
659
660 if (r[j][j].is_lt_zero())
661 a.negate();
662
663 for (k = j; k < rows; k++)
664 v[k].assign(r[k][j]);
665
666 LiDIA::add(v[j], v[j], a);
667
668 for (k = j; k < rows; k++)
669 LiDIA::multiply(u[k], v[k], c);
670
671 for (k = j; k < columns; k++) {
672 w[k].assign_zero();
673 x[k].assign_zero();
674 for (i = j; i < rows; i++) {
675 LiDIA::multiply(c, u[i], r[i][k]);
676 LiDIA::add(w[k], w[k], c);
677 LiDIA::multiply(c, u[i], q[i][k]);
678 LiDIA::add(x[k], x[k], c);
679 }
680 }
681
682 for (k = j; k < columns; k++)
683 for (i = j; i < rows; i++) {
684 LiDIA::multiply(c, v[i], w[k]);
685 LiDIA::subtract(r[i][k], r[i][k], c);
686 LiDIA::multiply(c, v[i], x[k]);
687 LiDIA::subtract(q[i][k], q[i][k], c);
688 }
689
690 for (k = 0; k < j; k++) {
691 x[k].assign_zero();
692 for (i = j; i < rows; i++) {
693 LiDIA::multiply(c, u[i], q[i][k]);
694 LiDIA::add(x[k], x[k], c);
695 }
696 }
697
698 for (k = 0; k < j; k++)
699 for (i = j; i < rows; i++) {
700 LiDIA::multiply(c, v[i], x[k]);
701 LiDIA::subtract(q[i][k], q[i][k], c);
702 }
703 }
704
705 for (i = 1; i < rows; i++)
706 for (j = 0; j < i; j++)
707 r[i][j].assign_zero();
708
709 for (i = 0; i < rows; i++)
710 for (j = 0; j < columns; j++)
711 (r[i][j]).approx(R.value[i][j]);
712
713 for (i = 0; i < rows; i++)
714 for (j = 0; j < columns; j++)
715 (q[i][j]).approx(Q.value[j][i]);
716
717 for (i = 0; i < rows; i++) {
718 delete[] q[i];
719 delete[] r[i];
720 }
721 delete[] q;
722 delete[] r;
723 }
724
725 int
equal(const bigfloat_matrix & A) const726 bigfloat_matrix::equal(const bigfloat_matrix & A) const
727 {
728 debug_handler("bigfloat_matrix", "equal(A)");
729 if (rows != A.rows || columns != A.columns) {
730 return 0;
731 }
732 lidia_size_t i = 0, r, c;
733 bigfloat erg;
734 for (r = 0; r < rows; r++)
735 for (c = 0; c < columns; c++) {
736 LiDIA::subtract(erg, value[r][c], A.value[r][c]);
737 erg = abs(erg);
738 if (erg.is_approx_zero())
739 i++;
740 }
741 return (i == rows*columns);
742 }
743
744
745 void
randomize()746 bigfloat_matrix::randomize()
747 {
748 debug_handler("bigfloat_matrix", "randomize()");
749 for (int l = 0; l < rows; l++)
750 for (int k = 0; k < columns; k++) {
751 value[k][l].assign(bigfloat(LiDIA::randomize(bigint(10000000))));
752 }
753 }
754
755
756
757 #ifdef LIDIA_NAMESPACE
758 } // end of namespace LiDIA
759 #endif
760