source: git/libpolys/coeffs/bigintmat.cc @ 6caad65

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