source: git/libpolys/coeffs/bigintmat.cc @ 75f10d

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