# source:git/libpolys/coeffs/bigintmat.cc@510dbc

spielwiese
Last change on this file since 510dbc was 510dbc, checked in by Hans Schoenemann <hannes@…>, 11 years ago
• Property mode set to `100644`
File size: 14.4 KB
Line
1/*****************************************
2 *  Computer Algebra System SINGULAR      *
3 *****************************************/
4/*
5 * ABSTRACT: class bigintmat: matrizes of big integers
6 */
7#include "bigintmat.h"
8
9#include <misc/intvec.h>
10
11//#include <kernel/mod2.h>
12//#include <kernel/options.h>
13
14#include <math.h>
15#include <string.h>
16
17#define BIMATELEM(M,I,J) (M)[ (M).index(I,J) ]
18
19bigintmat * bigintmat::transpose()
20{
21  bigintmat * t = new bigintmat(col, row, basecoeffs());
22  for (int i=1; i<=row; i++)
23  {
24    for (int j=1; j<=col; j++)
25    {
26      t->set(j, i, v[(i-1)*col+(j-1)]);
27    }
28  }
29  return t;
30}
31
32// Beginnt bei [0]
33void bigintmat::set(int i, number n, const coeffs C)
34{
35  assume (C == NULL || C == basecoeffs());
36
37  rawset(i, n_Copy(n, basecoeffs()), basecoeffs());
38}
39
40// Beginnt bei [1,1]
41void bigintmat::set(int i, int j, number n, const coeffs C)
42{
43  assume (C == NULL || C == basecoeffs());
44  assume (i > 0 && j > 0);
45  assume (i <= rows() && j <= cols());
46
47  set(index(i, j), n, C);
48}
49
50number bigintmat::get(int i) const
51{
52  assume (i >= 0);
53  assume (i<rows()*cols());
54
55  return n_Copy(v[i], basecoeffs());
56}
57
58number bigintmat::get(int i, int j) const
59{
60  assume (i > 0 && j > 0);
61  assume (i <= rows() && j <= cols());
62
63  return get(index(i, j));
64}
65
66// Ãberladener *=-Operator (fÃŒr int und bigint)
67// Frage hier: *= verwenden oder lieber = und * einzeln?
68void bigintmat::operator*=(int intop)
69{
70  number iop = n_Init(intop, basecoeffs());
71
72  inpMult(iop, basecoeffs());
73
74  n_Delete(&iop, basecoeffs());
75}
76
77void bigintmat::inpMult(number bintop, const coeffs C)
78{
79  assume (C == NULL || C == basecoeffs());
80
81  const int l = rows() * cols();
82
83  for (int i=0; i < l; i++)
84    n_InpMult(v[i], bintop, basecoeffs());
85}
86
87// Stimmen Parameter?
88// Welche der beiden Methoden?
89// Oder lieber eine comp-Funktion?
90
91bool operator==(const bigintmat & lhr, const bigintmat & rhr)
92{
93  if (&lhr == &rhr) { return true; }
94  if (lhr.cols() != rhr.cols()) { return false; }
95  if (lhr.rows() != rhr.rows()) { return false; }
96  if (lhr.basecoeffs() != rhr.basecoeffs()) { return false; }
97
98  const int l = (lhr.rows())*(lhr.cols());
99
100  for (int i=0; i < l; i++)
101  {
102    if (!n_Equal(lhr[i], rhr[i], lhr.basecoeffs())) { return false; }
103  }
104
105  return true;
106}
107
108bool operator!=(const bigintmat & lhr, const bigintmat & rhr)
109{
110  return !(lhr==rhr);
111}
112
113// Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ?
114bigintmat * bimAdd(bigintmat * a, bigintmat * b)
115{
116  if (a->cols() != b->cols()) return NULL;
117  if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
118
119  const int mn = si_min(a->rows(),b->rows());
120  const int ma = si_max(a->rows(),b->rows());
121
122  const coeffs basecoeffs = a->basecoeffs();
123
124  int i;
125
126  if (a->cols() == 1)
127  {
128    bigintmat * bim = new bigintmat(ma, 1, basecoeffs);
129
130    for (i=0; i<mn; i++)
131      bim->rawset(i, n_Add((*a)[i], (*b)[i], basecoeffs), basecoeffs);
132
133    if (ma > mn)
134    {
135      if (ma == a->rows())
136        for(i=mn; i<ma; i++)
137          bim->set(i, (*a)[i], basecoeffs);
138      else
139        for(i=mn; i<ma; i++)
140          bim->set(i, (*b)[i], basecoeffs);
141    }
142    return bim;
143  }
144
145  if (mn != ma) return NULL;
146
147  bigintmat * bim = new bigintmat(mn, a->cols(), basecoeffs);
148
149  for (i=0; i<mn*a->cols(); i++)
150    bim->rawset(i, n_Add((*a)[i], (*b)[i], basecoeffs), basecoeffs);
151
152  return bim;
153}
154bigintmat * bimAdd(bigintmat * a, int b)
155{
156
157  const int mn = a->rows()*a->cols();
158
159  const coeffs basecoeffs = a->basecoeffs();
160  number bb=n_Init(b,basecoeffs);
161
162  int i;
163
164  bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
165
166  for (i=0; i<mn; i++)
167    bim->rawset(i, n_Add((*a)[i], bb, basecoeffs), basecoeffs);
168
169  n_Delete(&bb,basecoeffs);
170  return bim;
171}
172
173bigintmat * bimSub(bigintmat * a, bigintmat * b)
174{
175  if (a->cols() != b->cols()) return NULL;
176  if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
177
178  const int mn = si_min(a->rows(),b->rows());
179  const int ma = si_max(a->rows(),b->rows());
180
181  const coeffs basecoeffs = a->basecoeffs();
182
183  int i;
184
185  if (a->cols() == 1)
186  {
187    bigintmat * bim = new bigintmat(ma, 1, basecoeffs);
188
189    for (i=0; i<mn; i++)
190      bim->rawset(i, n_Sub((*a)[i], (*b)[i], basecoeffs), basecoeffs);
191
192    if (ma > mn)
193    {
194      if (ma == a->rows())
195      {
196        for(i=mn; i<ma; i++)
197          bim->set(i, (*a)[i], basecoeffs);
198      }
199      else
200        for(i=mn; i<ma; i++)
201        {
202          number n = b->get(i);
203          n = n_Neg(n, basecoeffs);
204          bim->rawset(i, n, basecoeffs);
205        }
206    }
207    return bim;
208  }
209
210  if (mn != ma) return NULL;
211
212  bigintmat * bim = new bigintmat(mn, a->cols(), basecoeffs);
213
214  for (i=0; i<mn*a->cols(); i++)
215    bim->rawset(i, n_Sub((*a)[i], (*b)[i], basecoeffs), basecoeffs);
216
217  return bim;
218}
219
220bigintmat * bimSub(bigintmat * a, int b)
221{
222
223  const int mn = a->rows()*a->cols();
224
225  const coeffs basecoeffs = a->basecoeffs();
226  number bb=n_Init(b,basecoeffs);
227
228  int i;
229
230  bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
231
232  for (i=0; i<mn; i++)
233    bim->rawset(i, n_Sub((*a)[i], bb, basecoeffs), basecoeffs);
234
235  n_Delete(&bb,basecoeffs);
236  return bim;
237}
238
239bigintmat * bimMult(bigintmat * a, bigintmat * b)
240{
241  const int ca = a->cols();
242  const int cb = b->cols();
243
244  const int ra = a->rows();
245  const int rb = b->rows();
246
247  if (ca != rb)
248  {
249#ifndef NDEBUG
250    Werror("wrong bigintmat sizes at multiplication a * b: acols: %d != brows: %d\n", ca, rb);
251#endif
252    return NULL;
253  }
254
255  assume (ca == rb);
256
257  if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
258
259  const coeffs basecoeffs = a->basecoeffs();
260
261  int i, j, k;
262
263  number sum;
264
265  bigintmat * bim = new bigintmat(ra, cb, basecoeffs);
266
267  for (i=1; i<=ra; i++)
268    for (j=1; j<=cb; j++)
269    {
270      sum = n_Init(0, basecoeffs);
271
272      for (k=1; k<=ca; k++)
273      {
274        number prod = n_Mult( BIMATELEM(*a, i, k), BIMATELEM(*b, k, j), basecoeffs);
275
276        number sum2 = n_Add(sum, prod, basecoeffs); // no inplace add :(
277
278        n_Delete(&sum, basecoeffs); n_Delete(&prod, basecoeffs);
279
280        sum = sum2;
281      }
282      bim->rawset(i, j, sum, basecoeffs);
283    }
284  return bim;
285}
286
287bigintmat * bimMult(bigintmat * a, int b)
288{
289
290  const int mn = a->rows()*a->cols();
291
292  const coeffs basecoeffs = a->basecoeffs();
293  number bb=n_Init(b,basecoeffs);
294
295  int i;
296
297  bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
298
299  for (i=0; i<mn; i++)
300    bim->rawset(i, n_Mult((*a)[i], bb, basecoeffs), basecoeffs);
301
302  n_Delete(&bb,basecoeffs);
303  return bim;
304}
305
306bigintmat * bimMult(bigintmat * a, number b, const coeffs cf)
307{
308  if (cf!=a->basecoeffs()) return NULL;
309
310  const int mn = a->rows()*a->cols();
311
312  const coeffs basecoeffs = a->basecoeffs();
313
314  int i;
315
316  bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
317
318  for (i=0; i<mn; i++)
319    bim->rawset(i, n_Mult((*a)[i], b, basecoeffs), basecoeffs);
320
321  return bim;
322}
323
324// ----------------------------------------------------------------- //
325// Korrekt?
326
327intvec * bim2iv(bigintmat * b)
328{
329  intvec * iv = new intvec(b->rows(), b->cols(), 0);
330  for (int i=0; i<(b->rows())*(b->cols()); i++)
331    (*iv)[i] = n_Int((*b)[i], b->basecoeffs()); // Geht das so?
332  return iv;
333}
334
335bigintmat * iv2bim(intvec * b, const coeffs C)
336{
337  const int l = (b->rows())*(b->cols());
338  bigintmat * bim = new bigintmat(b->rows(), b->cols(), C);
339
340  for (int i=0; i < l; i++)
341    bim->rawset(i, n_Init((*b)[i], C), C);
342
343  return bim;
344}
345
346// ----------------------------------------------------------------- //
347
348int bigintmat::compare(const bigintmat* op) const
349{
350  assume (basecoeffs() == op->basecoeffs() );
351
352#ifndef NDEBUG
353  if (basecoeffs() != op->basecoeffs() )
354    WerrorS("wrong bigintmat comparison: different basecoeffs!\n");
355#endif
356
357  if ((col!=1) ||(op->cols()!=1))
358  {
359    if((col!=op->cols())
360       || (row!=op->rows()))
361      return -2;
362  }
363
364  int i;
365  for (i=0; i<si_min(row*col,op->rows()*op->cols()); i++)
366  {
367    if ( n_Greater(v[i], (*op)[i], basecoeffs()) )
368      return 1;
369    else if (! n_Equal(v[i], (*op)[i], basecoeffs()))
370      return -1;
371  }
372
373  for (; i<row; i++)
374  {
375    if ( n_GreaterZero(v[i], basecoeffs()) )
376      return 1;
377    else if (! n_IsZero(v[i], basecoeffs()) )
378      return -1;
379  }
380  for (; i<op->rows(); i++)
381  {
382    if ( n_GreaterZero((*op)[i], basecoeffs()) )
383      return -1;
384    else if (! n_IsZero((*op)[i], basecoeffs()) )
385      return 1;
386  }
387  return 0;
388}
389
390
391bigintmat * bimCopy(const bigintmat * b)
392{
393  if (b == NULL)
394    return NULL;
395
396  return new bigintmat(b);
397}
398
399char* bigintmat::String()
400{
401  assume (rows() > 0);
402  assume (cols() > 0);
403
404  StringSetS("");
405  const int l = cols()*rows();
406
407  n_Write(v[0], basecoeffs());
408  for (int i = 1; i < l; i++)
409  {
410    StringAppendS(","); n_Write(v[i], basecoeffs());
411  }
412  /* if (i != col*row-1)
413  {
414  StringAppendS(",");
415  if ((i+1)%col == 0)
416  StringAppendS("\n");
417  }   */
418  return StringAppendS("");
419}
420
421int intArrSum(int * a, int length)
422{
423  int sum = 0;
424  for (int i=0; i<length; i++)
425    sum += a[i];
426  return sum;
427}
428
429int findLongest(int * a, int length)
430{
431  int l = 0;
432  int index;
433  for (int i=0; i<length; i++)
434  {
435    if (a[i] > l)
436    {
437      l = a[i];
438      index = i;
439    }
440  }
441  return index;
442}
443
444int getShorter (int * a, int l, int j, int cols, int rows)
445{
446  int sndlong = 0;
447  int min;
448  for (int i=0; i<rows; i++)
449  {
450    int index = cols*i+j;
451    if ((a[index] > sndlong) && (a[index] < l))
452    {
453      min = floor(log10(cols))+floor(log10(rows))+5;
454      if ((a[index] < min) && (min < l))
455        sndlong = min;
456      else
457        sndlong = a[index];
458    }
459  }
460  if (sndlong == 0)
461  {
462    min = floor(log10(cols))+floor(log10(rows))+5;
463    if (min < l)
464      sndlong = min;
465    else
466      sndlong = 1;
467  }
468  return sndlong;
469}
470
471
472int * bigintmat::getwid(int maxwid)
473{
474  int const c = /*2**/(col-1)+1;
475  if (col + c > maxwid-1) return NULL;
476  int * wv = (int*)omAlloc(sizeof(int)*col*row);
477  int * cwv = (int*)omAlloc(sizeof(int)*col);
478  for (int j=0; j<col; j++)
479  {
480    cwv[j] = 0;
481    for (int i=0; i<row; i++)
482    {
483      StringSetS("");
484      n_Write(v[col*i+j], basecoeffs());
485      char * tmp = StringAppendS("");
486      const int _nl = strlen(tmp);
487      wv[col*i+j] = _nl;
488      if (_nl > cwv[j])
489        cwv[j]=_nl;
490    }
491  }
492
493  // Groesse verkleinern, bis < maxwid
494  while (intArrSum(cwv, col)+c > maxwid)
495  {
496    int j = findLongest(cwv, col);
497    cwv[j] = getShorter(wv, cwv[j], j, col, row);
498  }
499  omFree(wv);
500  return cwv;
501}
502
503void bigintmat::pprint(int maxwid)
504{
505  if ((col==0) || (row==0))
506    PrintS("");
507  else
508  {
509    int * colwid = getwid(maxwid);
510    if (colwid == NULL)
511    {
512      WerrorS("not enough space to print bigintmat");
513      return;
514    }
515    char * ps;
516    int slength = 0;
517    for (int j=0; j<col; j++)
518      slength += colwid[j]*row;
519    slength += col*row+row;
520    ps = (char*) omAlloc0(sizeof(char)*(slength));
521    int pos = 0;
522    for (int i=0; i<col*row; i++)
523    {
524      StringSetS("");
525      n_Write(v[i], basecoeffs());
526      char * temp = StringAppendS("");
527      char * ts = omStrDup(temp);
528      const int _nl = strlen(ts);
529      int cj = i%col;
530      if (_nl > colwid[cj])
531      {
532        StringSetS("");
533        int ci = floor(i/col);
534        StringAppend("[%d,%d]", ci+1, cj+1);
535        char *tmp = StringAppendS("");
536        char * ph = omStrDup(tmp);
537        int phl = strlen(ph);
538        if (phl > colwid[cj])
539        {
540          for (int j=0; j<colwid[cj]-1; j++)
541            ps[pos+j] = ' ';
542          ps[pos+colwid[cj]-1] = '*';
543        }
544        else
545        {
546          for (int j=0; j<colwid[cj]-phl; j++)
547            ps[pos+j] = ' ';
548          for (int j=0; j<phl; j++)
549            ps[pos+colwid[cj]-phl+j] = ph[j];
550        }
551        omFree(ph);
552      }
553      else  // Mit Leerzeichen auffÃŒllen und zahl reinschreiben
554      {
555        for (int j=0; j<colwid[cj]-_nl; j++)
556          ps[pos+j] = ' ';
557        for (int j=0; j<_nl; j++)
558          ps[pos+colwid[cj]-_nl+j] = ts[j];
559      }
560      // ", " und (evtl) "\n" einfÃŒgen
561      if ((i+1)%col == 0)
562      {
563        if (i != col*row-1)
564        {
565          ps[pos+colwid[cj]] = ',';
566          ps[pos+colwid[cj]+1] = '\n';
567          pos += colwid[cj]+2;
568        }
569      }
570      else
571      {
572        ps[pos+colwid[cj]] = ',';
573        pos += colwid[cj]+1;
574      }
575
576      omFree(ts);  // Hier ts zerstÃ¶ren
577    }
578    PrintS(ps);
579   // omFree(ps);
580  }
581}
582
583// Ungetestet
584static void bimRowContent(bigintmat *bimat, int rowpos, int colpos)
585{
586  const coeffs basecoeffs = bimat->basecoeffs();
587
588  number tgcd, m;
589  int i=bimat->cols();
590
591  loop
592  {
593    tgcd = n_Copy(BIMATELEM(*bimat,rowpos,i--), basecoeffs);
594    if (! n_IsZero(tgcd, basecoeffs)) break;
595    if (i<colpos) return;
596  }
597  if ((! n_GreaterZero(tgcd, basecoeffs)) && (! n_IsZero(tgcd, basecoeffs))) tgcd = n_Neg(tgcd, basecoeffs);
598  if ( n_IsOne(tgcd,basecoeffs)) return;
599  loop
600  {
601    m = n_Copy(BIMATELEM(*bimat,rowpos,i--), basecoeffs);
602    if (! n_IsZero(m,basecoeffs))
603    {
604      number tp1 = n_Gcd(tgcd, m, basecoeffs);
605      n_Delete(&tgcd, basecoeffs);
606      tgcd = tp1;
607    }
608    if ( n_IsOne(tgcd,basecoeffs)) return;
609    if (i<colpos) break;
610  }
611  for (i=bimat->cols();i>=colpos;i--)
612  {
613    number tp2 = n_Div(BIMATELEM(*bimat,rowpos,i), tgcd,basecoeffs);
614    n_Delete(&BIMATELEM(*bimat,rowpos,i), basecoeffs);
615    BIMATELEM(*bimat,rowpos,i) = tp2;
616  }
617  n_Delete(&tgcd, basecoeffs);
618  n_Delete(&m, basecoeffs);
619}
620
621static void bimReduce(bigintmat *bimat, int rpiv, int colpos,
623{
624  const coeffs basecoeffs = bimat->basecoeffs();
625
626  number tgcd, ce, m1, m2;
627  int j, i;
628  number piv = BIMATELEM(*bimat,rpiv,colpos);
629
631  {
632    ce = n_Copy(BIMATELEM(*bimat,j,colpos),basecoeffs);
633    if (! n_IsZero(ce, basecoeffs))
634    {
635      n_Delete(&BIMATELEM(*bimat,j,colpos), basecoeffs);
636      BIMATELEM(*bimat,j,colpos) = n_Init(0, basecoeffs);
637      m1 = n_Copy(piv,basecoeffs);
638      m2 = n_Copy(ce,basecoeffs);
639      tgcd = n_Gcd(m1, m2, basecoeffs);
640      if (! n_IsOne(tgcd,basecoeffs))
641      {
642        number tp1 = n_Div(m1, tgcd,basecoeffs);
643        number tp2 = n_Div(m2, tgcd,basecoeffs);
644        n_Delete(&m1, basecoeffs);
645        n_Delete(&m2, basecoeffs);
646        m1 = tp1;
647        m2 = tp2;
648      }
649      for (i=bimat->cols();i>colpos;i--)
650      {
651        n_Delete(&BIMATELEM(*bimat,j,i), basecoeffs);
652        number tp1 = n_Mult(BIMATELEM(*bimat,j,i), m1,basecoeffs);
653        number tp2 = n_Mult(BIMATELEM(*bimat,rpiv,i), m2,basecoeffs);
654        BIMATELEM(*bimat,j,i) = n_Sub(tp1, tp2,basecoeffs);
655        n_Delete(&tp1, basecoeffs);
656        n_Delete(&tp2, basecoeffs);
657      }
658      bimRowContent(bimat, j, colpos+1);
659      n_Delete(&m1, basecoeffs);
660      n_Delete(&m2, basecoeffs);
661    }
662    n_Delete(&ce, basecoeffs);
663  }
664}
665
666
667
668
669
Note: See TracBrowser for help on using the repository browser.