source: git/libpolys/coeffs/bigintmat.cc @ eea50a

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