My Project
Macros | Functions
bigintmat.cc File Reference
#include "misc/auxiliary.h"
#include "coeffs/bigintmat.h"
#include "misc/intvec.h"
#include "coeffs/rmodulon.h"
#include <cmath>

Go to the source code of this file.

Macros

#define swap(_i, _j)
 
#define MIN(a, b)   (a < b ? a : b)
 

Functions

static coeffs numbercoeffs (number n, coeffs c)
 create Z/nA of type n_Zn More...
 
bool operator== (const bigintmat &lhr, const bigintmat &rhr)
 
bool operator!= (const bigintmat &lhr, const bigintmat &rhr)
 
bigintmatbimAdd (bigintmat *a, bigintmat *b)
 Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ? @Note: NULL as a result means an error (non-compatible matrices?) More...
 
bigintmatbimAdd (bigintmat *a, long b)
 
bigintmatbimSub (bigintmat *a, bigintmat *b)
 
bigintmatbimSub (bigintmat *a, long b)
 
bigintmatbimMult (bigintmat *a, bigintmat *b)
 
bigintmatbimMult (bigintmat *a, long b)
 
bigintmatbimMult (bigintmat *a, number b, const coeffs cf)
 
intvecbim2iv (bigintmat *b)
 
bigintmativ2bim (intvec *b, const coeffs C)
 
bigintmatbimCopy (const bigintmat *b)
 same as copy constructor - apart from it being able to accept NULL as input More...
 
static int intArrSum (int *a, int length)
 
static int findLongest (int *a, int length)
 
static int getShorter (int *a, int l, int j, int cols, int rows)
 
bigintmatbimChangeCoeff (bigintmat *a, coeffs cnew)
 Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen. More...
 
void bimMult (bigintmat *a, bigintmat *b, bigintmat *c)
 Multipliziert Matrix a und b und speichert Ergebnis in c. More...
 
static void reduce_mod_howell (bigintmat *A, bigintmat *b, bigintmat *eps, bigintmat *x)
 
static bigintmatprependIdentity (bigintmat *A)
 
static number bimFarey (bigintmat *A, number N, bigintmat *L)
 
static number solveAx_dixon (bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern)
 
static number solveAx_howell (bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
 
number solveAx (bigintmat *A, bigintmat *b, bigintmat *x)
 solve Ax=b*d. x needs to be pre-allocated to the same number of columns as b. the minimal denominator d is returned. Currently available for Z, Q and Z/nZ (and possibly for all fields: d=1 there) Beware that the internal functions can find the kernel as well - but the interface is lacking. More...
 
void diagonalForm (bigintmat *A, bigintmat **S, bigintmat **T)
 
int kernbase (bigintmat *a, bigintmat *c, number p, coeffs q)
 a basis for the nullspace of a mod p: only used internally in Round2. Don't use it. More...
 
bool nCoeffs_are_equal (coeffs r, coeffs s)
 

Macro Definition Documentation

◆ MIN

#define MIN (   a,
  b 
)    (a < b ? a : b)

◆ swap

#define swap (   _i,
  _j 
)
Value:
int __i = (_i), __j=(_j); \
number c = v[__i]; \
v[__i] = v[__j]; \
v[__j] = c \
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39

Function Documentation

◆ bim2iv()

intvec* bim2iv ( bigintmat b)

Definition at line 341 of file bigintmat.cc.

342 {
343  intvec * iv = new intvec(b->rows(), b->cols(), 0);
344  for (int i=0; i<(b->rows())*(b->cols()); i++)
345  (*iv)[i] = n_Int((*b)[i], b->basecoeffs()); // Geht das so?
346  return iv;
347 }
int i
Definition: cfEzgcd.cc:132
CanonicalForm b
Definition: cfModGcd.cc:4103
Definition: intvec.h:23
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition: coeffs.h:547

◆ bimAdd() [1/2]

bigintmat* bimAdd ( bigintmat a,
bigintmat b 
)

Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ? @Note: NULL as a result means an error (non-compatible matrices?)

Definition at line 182 of file bigintmat.cc.

183 {
184  if (a->cols() != b->cols()) return NULL;
185  if (a->rows() != b->rows()) return NULL;
186  if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
187 
188  const coeffs basecoeffs = a->basecoeffs();
189 
190  int i;
191 
192  bigintmat * bim = new bigintmat(a->rows(), a->cols(), basecoeffs);
193 
194  for (i=a->rows()*a->cols()-1;i>=0; i--)
195  bim->rawset(i, n_Add((*a)[i], (*b)[i], basecoeffs), basecoeffs);
196 
197  return bim;
198 }
Matrices of numbers.
Definition: bigintmat.h:51
int cols() const
Definition: bigintmat.h:144
int rows() const
Definition: bigintmat.h:145
void rawset(int i, number n, const coeffs C=NULL)
replace an entry with the given number n (only delete old). NOTE: starts at [0]. Should be named set_...
Definition: bigintmat.h:196
coeffs basecoeffs() const
Definition: bigintmat.h:146
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
Definition: coeffs.h:650
The main handler for Singular numbers which are suitable for Singular polynomials.
#define NULL
Definition: omList.c:12

◆ bimAdd() [2/2]

bigintmat* bimAdd ( bigintmat a,
long  b 
)

Definition at line 199 of file bigintmat.cc.

200 {
201 
202  const int mn = si_min(a->rows(),a->cols());
203 
204  const coeffs basecoeffs = a->basecoeffs();
205  number bb=n_Init(b,basecoeffs);
206 
207  int i;
208 
209  bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
210 
211  for (i=1; i<=mn; i++)
212  BIMATELEM(*bim,i,i)=n_Add(BIMATELEM(*a,i,i), bb, basecoeffs);
213 
214  n_Delete(&bb,basecoeffs);
215  return bim;
216 }
static int si_min(const int a, const int b)
Definition: auxiliary.h:125
#define BIMATELEM(M, I, J)
Definition: bigintmat.h:133
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:455
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:538

◆ bimChangeCoeff()

bigintmat* bimChangeCoeff ( bigintmat a,
coeffs  cnew 
)

Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.

Definition at line 1804 of file bigintmat.cc.

1805 {
1806  coeffs cold = a->basecoeffs();
1807  bigintmat *b = new bigintmat(a->rows(), a->cols(), cnew);
1808  // Erzeugt Karte von alten coeffs nach neuen
1809  nMapFunc f = n_SetMap(cold, cnew);
1810  number t1;
1811  number t2;
1812  // apply map to all entries.
1813  for (int i=1; i<=a->rows(); i++)
1814  {
1815  for (int j=1; j<=a->cols(); j++)
1816  {
1817  t1 = a->get(i, j);
1818  t2 = f(t1, cold, cnew);
1819  b->set(i, j, t2);
1820  n_Delete(&t1, cold);
1821  n_Delete(&t2, cnew);
1822  }
1823  }
1824  return b;
1825 }
FILE * f
Definition: checklibs.c:9
number get(int i, int j) const
get a copy of an entry. NOTE: starts at [1,1]
Definition: bigintmat.cc:119
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition: coeffs.h:700
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
int j
Definition: facHensel.cc:110

◆ bimCopy()

bigintmat* bimCopy ( const bigintmat b)

same as copy constructor - apart from it being able to accept NULL as input

Definition at line 405 of file bigintmat.cc.

406 {
407  if (b == NULL)
408  return NULL;
409 
410  return new bigintmat(b);
411 }

◆ bimFarey()

static number bimFarey ( bigintmat A,
number  N,
bigintmat L 
)
static

Definition at line 2048 of file bigintmat.cc.

2049 {
2050  coeffs Z = A->basecoeffs(),
2051  Q = nInitChar(n_Q, 0);
2052  number den = n_Init(1, Z);
2053  nMapFunc f = n_SetMap(Q, Z);
2054 
2055  for(int i=1; i<= A->rows(); i++)
2056  {
2057  for(int j=1; j<= A->cols(); j++)
2058  {
2059  number ad = n_Mult(den, A->view(i, j), Z);
2060  number re = n_IntMod(ad, N, Z);
2061  n_Delete(&ad, Z);
2062  number q = n_Farey(re, N, Z);
2063  n_Delete(&re, Z);
2064  if (!q)
2065  {
2066  n_Delete(&ad, Z);
2067  n_Delete(&den, Z);
2068  return NULL;
2069  }
2070 
2071  number d = n_GetDenom(q, Q),
2072  n = n_GetNumerator(q, Q);
2073 
2074  n_Delete(&q, Q);
2075  n_Delete(&ad, Z);
2076  number dz = f(d, Q, Z),
2077  nz = f(n, Q, Z);
2078  n_Delete(&d, Q);
2079  n_Delete(&n, Q);
2080 
2081  if (!n_IsOne(dz, Z))
2082  {
2083  L->skalmult(dz, Z);
2084  n_InpMult(den, dz, Z);
2085 #if 0
2086  PrintS("den increasing to ");
2087  n_Print(den, Z);
2088  PrintLn();
2089 #endif
2090  }
2091  n_Delete(&dz, Z);
2092  L->rawset(i, j, nz);
2093  }
2094  }
2095 
2096  nKillChar(Q);
2097  PrintS("bimFarey worked\n");
2098 #if 0
2099  L->Print();
2100  PrintS("\n * 1/");
2101  n_Print(den, Z);
2102  PrintLn();
2103 #endif
2104  return den;
2105 }
CanonicalForm den(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
void Print()
IO: simply prints the matrix to the current output (screen?)
Definition: bigintmat.cc:443
bool skalmult(number b, coeffs c)
Multipliziert zur Matrix den Skalar b hinzu.
Definition: bigintmat.cc:938
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
Definition: coeffs.h:636
static FORCE_INLINE number n_GetDenom(number &n, const coeffs r)
return the denominator of n (if elements of r are by nature not fractional, result is 1)
Definition: coeffs.h:603
@ n_Q
rational (GMP) numbers
Definition: coeffs.h:30
void n_Print(number &a, const coeffs r)
print a number (BEWARE of string buffers!) mostly for debugging
Definition: numbers.cc:646
static FORCE_INLINE number n_Farey(number a, number b, const coeffs r)
Definition: coeffs.h:767
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:392
static FORCE_INLINE number n_IntMod(number a, number b, const coeffs r)
for r a field, return n_Init(0,r) always: n_Div(a,b,r)*b+n_IntMod(a,b,r)==a n_IntMod(a,...
Definition: coeffs.h:628
static FORCE_INLINE void n_InpMult(number &a, number b, const coeffs r)
multiplication of 'a' and 'b'; replacement of 'a' by the product a*b
Definition: coeffs.h:641
static FORCE_INLINE number n_GetNumerator(number &n, const coeffs r)
return the numerator of n (if elements of r are by nature not fractional, result is n)
Definition: coeffs.h:608
void nKillChar(coeffs r)
undo all initialisations
Definition: numbers.cc:547
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:468
STATIC_VAR jList * Q
Definition: janet.cc:30
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
#define A
Definition: sirandom.c:24

◆ bimMult() [1/4]

bigintmat* bimMult ( bigintmat a,
bigintmat b 
)

Definition at line 255 of file bigintmat.cc.

256 {
257  const int ca = a->cols();
258  const int cb = b->cols();
259 
260  const int ra = a->rows();
261  const int rb = b->rows();
262 
263  if (ca != rb)
264  {
265 #ifndef SING_NDEBUG
266  Werror("wrong bigintmat sizes at multiplication a * b: acols: %d != brows: %d\n", ca, rb);
267 #endif
268  return NULL;
269  }
270 
271  assume (ca == rb);
272 
273  if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
274 
275  const coeffs basecoeffs = a->basecoeffs();
276 
277  int i, j, k;
278 
279  number sum;
280 
281  bigintmat * bim = new bigintmat(ra, cb, basecoeffs);
282 
283  for (i=1; i<=ra; i++)
284  for (j=1; j<=cb; j++)
285  {
286  sum = n_Init(0, basecoeffs);
287 
288  for (k=1; k<=ca; k++)
289  {
290  number prod = n_Mult( BIMATELEM(*a, i, k), BIMATELEM(*b, k, j), basecoeffs);
291 
292  n_InpAdd(sum, prod, basecoeffs);
293 
294  n_Delete(&prod, basecoeffs);
295  }
296  bim->rawset(i, j, sum, basecoeffs);
297  }
298  return bim;
299 }
int k
Definition: cfEzgcd.cc:99
static FORCE_INLINE void n_InpAdd(number &a, number b, const coeffs r)
addition of 'a' and 'b'; replacement of 'a' by the sum a+b
Definition: coeffs.h:646
fq_nmod_poly_t prod
Definition: facHensel.cc:100
#define assume(x)
Definition: mod2.h:389
void Werror(const char *fmt,...)
Definition: reporter.cc:189

◆ bimMult() [2/4]

void bimMult ( bigintmat a,
bigintmat b,
bigintmat c 
)

Multipliziert Matrix a und b und speichert Ergebnis in c.

Definition at line 1932 of file bigintmat.cc.

1933 {
1934  if (!nCoeffs_are_equal(a->basecoeffs(), b->basecoeffs()))
1935  {
1936  WerrorS("Error in bimMult. Coeffs do not agree!");
1937  return;
1938  }
1939  if ((a->rows() != c->rows()) || (b->cols() != c->cols()) || (a->cols() != b->rows()))
1940  {
1941  WerrorS("Error in bimMult. Dimensions do not agree!");
1942  return;
1943  }
1944  bigintmat *tmp = bimMult(a, b);
1945  c->copy(tmp);
1946 
1947  delete tmp;
1948 }
bool nCoeffs_are_equal(coeffs r, coeffs s)
Definition: bigintmat.cc:2645
bigintmat * bimMult(bigintmat *a, bigintmat *b)
Definition: bigintmat.cc:255
bool copy(bigintmat *b)
Kopiert Einträge von b auf Bigintmat.
Definition: bigintmat.cc:1259
void WerrorS(const char *s)
Definition: feFopen.cc:24

◆ bimMult() [3/4]

bigintmat* bimMult ( bigintmat a,
long  b 
)

Definition at line 301 of file bigintmat.cc.

302 {
303 
304  const int mn = a->rows()*a->cols();
305 
306  const coeffs basecoeffs = a->basecoeffs();
307  number bb=n_Init(b,basecoeffs);
308 
309  int i;
310 
311  bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
312 
313  for (i=0; i<mn; i++)
314  bim->rawset(i, n_Mult((*a)[i], bb, basecoeffs), basecoeffs);
315 
316  n_Delete(&bb,basecoeffs);
317  return bim;
318 }

◆ bimMult() [4/4]

bigintmat* bimMult ( bigintmat a,
number  b,
const coeffs  cf 
)

Definition at line 320 of file bigintmat.cc.

321 {
322  if (cf!=a->basecoeffs()) return NULL;
323 
324  const int mn = a->rows()*a->cols();
325 
326  const coeffs basecoeffs = a->basecoeffs();
327 
328  int i;
329 
330  bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
331 
332  for (i=0; i<mn; i++)
333  bim->rawset(i, n_Mult((*a)[i], b, basecoeffs), basecoeffs);
334 
335  return bim;
336 }
CanonicalForm cf
Definition: cfModGcd.cc:4083

◆ bimSub() [1/2]

bigintmat* bimSub ( bigintmat a,
bigintmat b 
)

Definition at line 218 of file bigintmat.cc.

219 {
220  if (a->cols() != b->cols()) return NULL;
221  if (a->rows() != b->rows()) return NULL;
222  if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
223 
224  const coeffs basecoeffs = a->basecoeffs();
225 
226  int i;
227 
228  bigintmat * bim = new bigintmat(a->rows(), a->cols(), basecoeffs);
229 
230  for (i=a->rows()*a->cols()-1;i>=0; i--)
231  bim->rawset(i, n_Sub((*a)[i], (*b)[i], basecoeffs), basecoeffs);
232 
233  return bim;
234 }
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of 'a' and 'b', i.e., a-b
Definition: coeffs.h:655

◆ bimSub() [2/2]

bigintmat* bimSub ( bigintmat a,
long  b 
)

Definition at line 236 of file bigintmat.cc.

237 {
238  const int mn = si_min(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=1; i<=mn; i++)
248  BIMATELEM(*bim,i,i)=n_Sub(BIMATELEM(*a,i,i), bb, basecoeffs);
249 
250  n_Delete(&bb,basecoeffs);
251  return bim;
252 }

◆ diagonalForm()

void diagonalForm ( bigintmat A,
bigintmat **  S,
bigintmat **  T 
)

Definition at line 2475 of file bigintmat.cc.

2476 {
2477  bigintmat * t, *s, *a=A;
2478  coeffs R = a->basecoeffs();
2479  if (T)
2480  {
2481  *T = new bigintmat(a->cols(), a->cols(), R),
2482  (*T)->one();
2483  t = new bigintmat(*T);
2484  }
2485  else
2486  {
2487  t = *T;
2488  }
2489 
2490  if (S)
2491  {
2492  *S = new bigintmat(a->rows(), a->rows(), R);
2493  (*S)->one();
2494  s = new bigintmat(*S);
2495  }
2496  else
2497  {
2498  s = *S;
2499  }
2500 
2501  int flip=0;
2502  do
2503  {
2504  bigintmat * x, *X;
2505  if (flip)
2506  {
2507  x = s;
2508  X = *S;
2509  }
2510  else
2511  {
2512  x = t;
2513  X = *T;
2514  }
2515 
2516  if (x)
2517  {
2518  x->one();
2519  bigintmat * r = new bigintmat(a->rows()+a->cols(), a->cols(), R);
2520  bigintmat * rw = new bigintmat(1, a->cols(), R);
2521  for(int i=0; i<a->cols(); i++)
2522  {
2523  x->getrow(i+1, rw);
2524  r->setrow(i+1, rw);
2525  }
2526  for (int i=0; i<a->rows(); i++)
2527  {
2528  a->getrow(i+1, rw);
2529  r->setrow(i+a->cols()+1, rw);
2530  }
2531  r->hnf();
2532  for(int i=0; i<a->cols(); i++)
2533  {
2534  r->getrow(i+1, rw);
2535  x->setrow(i+1, rw);
2536  }
2537  for(int i=0; i<a->rows(); i++)
2538  {
2539  r->getrow(i+a->cols()+1, rw);
2540  a->setrow(i+1, rw);
2541  }
2542  delete rw;
2543  delete r;
2544 
2545 #if 0
2546  Print("X: %ld\n", X);
2547  X->Print();
2548  Print("\nx: %ld\n", x);
2549  x->Print();
2550 #endif
2551  bimMult(X, x, X);
2552 #if 0
2553  Print("\n2:X: %ld %ld %ld\n", X, *S, *T);
2554  X->Print();
2555  Print("\n2:x: %ld\n", x);
2556  x->Print();
2557  PrintLn();
2558 #endif
2559  }
2560  else
2561  {
2562  a->hnf();
2563  }
2564 
2565  int diag = 1;
2566  for(int i=a->rows(); diag && i>0; i--)
2567  {
2568  for(int j=a->cols(); j>0; j--)
2569  {
2570  if ((a->rows()-i)!=(a->cols()-j) && !n_IsZero(a->view(i, j), R))
2571  {
2572  diag = 0;
2573  break;
2574  }
2575  }
2576  }
2577 #if 0
2578  PrintS("Diag ? %d\n", diag);
2579  a->Print();
2580  PrintLn();
2581 #endif
2582  if (diag) break;
2583 
2584  a = a->transpose(); // leaks - I need to write inpTranspose
2585  flip = 1-flip;
2586  } while (1);
2587  if (flip)
2588  a = a->transpose();
2589 
2590  if (S) *S = (*S)->transpose();
2591  if (s) delete s;
2592  if (t) delete t;
2593  A->copy(a);
2594 }
Variable x
Definition: cfModGcd.cc:4082
void hnf()
transforms INPLACE to HNF
Definition: bigintmat.cc:1660
bigintmat * transpose()
Definition: bigintmat.cc:37
void setrow(int i, bigintmat *m)
Setzt i-te Zeile gleich übergebenem Vektor (Matrix) m.
Definition: bigintmat.cc:860
number view(int i, int j) const
view an entry an entry. NOTE: starts at [1,1]
Definition: bigintmat.cc:127
void one()
Macht Matrix (Falls quadratisch) zu Einheitsmatrix.
Definition: bigintmat.cc:1325
void getrow(int i, bigintmat *a)
Schreibt i-te Zeile in Vektor (Matrix) a.
Definition: bigintmat.cc:791
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition: coeffs.h:464
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
std::pair< ideal, ring > flip(const ideal I, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal, const gfan::ZVector adjustedInteriorPoint, const gfan::ZVector adjustedFacetNormal)
Definition: flip.cc:17
STATIC_VAR jList * T
Definition: janet.cc:30
#define R
Definition: sirandom.c:27

◆ findLongest()

static int findLongest ( int *  a,
int  length 
)
static

Definition at line 537 of file bigintmat.cc.

538 {
539  int l = 0;
540  int index;
541  for (int i=0; i<length; i++)
542  {
543  if (a[i] > l)
544  {
545  l = a[i];
546  index = i;
547  }
548  }
549  return index;
550 }
int l
Definition: cfEzgcd.cc:100
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:257
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592

◆ getShorter()

static int getShorter ( int *  a,
int  l,
int  j,
int  cols,
int  rows 
)
static

Definition at line 552 of file bigintmat.cc.

553 {
554  int sndlong = 0;
555  int min;
556  for (int i=0; i<rows; i++)
557  {
558  int index = cols*i+j;
559  if ((a[index] > sndlong) && (a[index] < l))
560  {
561  min = floor(log10((double)cols))+floor(log10((double)rows))+5;
562  if ((a[index] < min) && (min < l))
563  sndlong = min;
564  else
565  sndlong = a[index];
566  }
567  }
568  if (sndlong == 0)
569  {
570  min = floor(log10((double)cols))+floor(log10((double)rows))+5;
571  if (min < l)
572  sndlong = min;
573  else
574  sndlong = 1;
575  }
576  return sndlong;
577 }
static int min(int a, int b)
Definition: fast_mult.cc:268
const signed long floor(const ampf< Precision > &x)
Definition: amp.h:773
const ampf< Precision > log10(const ampf< Precision > &x)
Definition: amp.h:1022

◆ intArrSum()

static int intArrSum ( int *  a,
int  length 
)
static

Definition at line 529 of file bigintmat.cc.

530 {
531  int sum = 0;
532  for (int i=0; i<length; i++)
533  sum += a[i];
534  return sum;
535 }

◆ iv2bim()

bigintmat* iv2bim ( intvec b,
const coeffs  C 
)

Definition at line 349 of file bigintmat.cc.

350 {
351  const int l = (b->rows())*(b->cols());
352  bigintmat * bim = new bigintmat(b->rows(), b->cols(), C);
353 
354  for (int i=0; i < l; i++)
355  bim->rawset(i, n_Init((*b)[i], C), C);
356 
357  return bim;
358 }

◆ kernbase()

int kernbase ( bigintmat a,
bigintmat c,
number  p,
coeffs  q 
)

a basis for the nullspace of a mod p: only used internally in Round2. Don't use it.

Definition at line 2600 of file bigintmat.cc.

2601 {
2602 #if 0
2603  PrintS("Kernel of ");
2604  a->Print();
2605  PrintS(" modulo ");
2606  n_Print(p, q);
2607  PrintLn();
2608 #endif
2609 
2610  coeffs coe = numbercoeffs(p, q);
2611  bigintmat *m = bimChangeCoeff(a, coe), *U, *V;
2612  diagonalForm(m, &U, &V);
2613 #if 0
2614  PrintS("\ndiag form: ");
2615  m->Print();
2616  PrintS("\nU:\n");
2617  U->Print();
2618  PrintS("\nV:\n");
2619  V->Print();
2620  PrintLn();
2621 #endif
2622 
2623  int rg = 0;
2624 #undef MIN
2625 #define MIN(a,b) (a < b ? a : b)
2626  for(rg=0; rg<MIN(m->rows(), m->cols()) && !n_IsZero(m->view(m->rows()-rg,m->cols()-rg), coe); rg++);
2627 
2628  bigintmat * k = new bigintmat(m->cols(), m->rows(), coe);
2629  for(int i=0; i<rg; i++)
2630  {
2631  number A = n_Ann(m->view(m->rows()-i, m->cols()-i), coe);
2632  k->set(m->cols()-i, i+1, A);
2633  n_Delete(&A, coe);
2634  }
2635  for(int i=rg; i<m->cols(); i++)
2636  {
2637  k->set(m->cols()-i, i+1-rg, n_Init(1, coe));
2638  }
2639  bimMult(V, k, k);
2640  c->copy(bimChangeCoeff(k, q));
2641  return c->cols();
2642 }
#define MIN(a, b)
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.
Definition: bigintmat.cc:1804
void diagonalForm(bigintmat *A, bigintmat **S, bigintmat **T)
Definition: bigintmat.cc:2475
static coeffs numbercoeffs(number n, coeffs c)
create Z/nA of type n_Zn
Definition: bigintmat.cc:21
int m
Definition: cfEzgcd.cc:128
int p
Definition: cfModGcd.cc:4078
static FORCE_INLINE number n_Ann(number a, const coeffs r)
if r is a ring with zero divisors, return an annihilator!=0 of b otherwise return NULL
Definition: coeffs.h:679

◆ nCoeffs_are_equal()

bool nCoeffs_are_equal ( coeffs  r,
coeffs  s 
)

Definition at line 2645 of file bigintmat.cc.

2646 {
2647  if ((r == NULL) || (s == NULL))
2648  return false;
2649  if (r == s)
2650  return true;
2651  if ((getCoeffType(r)==n_Z) && (getCoeffType(s)==n_Z))
2652  return true;
2653  if ((getCoeffType(r)==n_Zp) && (getCoeffType(s)==n_Zp))
2654  {
2655  if (r->ch == s->ch)
2656  return true;
2657  else
2658  return false;
2659  }
2660  // n_Zn stimmt wahrscheinlich noch nicht
2661  if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2662  {
2663  if (r->ch == s->ch)
2664  return true;
2665  else
2666  return false;
2667  }
2668  if ((getCoeffType(r)==n_Q) && (getCoeffType(s)==n_Q))
2669  return true;
2670  // FALL n_Zn FEHLT NOCH!
2671  //if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2672  return false;
2673 }
@ n_Zn
only used if HAVE_RINGS is defined
Definition: coeffs.h:44
@ n_Zp
\F{p < 2^31}
Definition: coeffs.h:29
@ n_Z
only used if HAVE_RINGS is defined
Definition: coeffs.h:43
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:421

◆ numbercoeffs()

static coeffs numbercoeffs ( number  n,
coeffs  c 
)
static

create Z/nA of type n_Zn

Definition at line 21 of file bigintmat.cc.

22 {
23  mpz_t p;
24  number2mpz(n, c, p);
25  ZnmInfo *pp = new ZnmInfo;
26  pp->base = p;
27  pp->exp = 1;
28  coeffs nc = nInitChar(n_Zn, (void*)pp);
29  mpz_clear(p);
30  delete pp;
31  return nc;
32 }
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
static FORCE_INLINE void number2mpz(number n, coeffs c, mpz_t m)
Definition: coeffs.h:987

◆ operator!=()

bool operator!= ( const bigintmat lhr,
const bigintmat rhr 
)

Definition at line 176 of file bigintmat.cc.

177 {
178  return !(lhr==rhr);
179 }

◆ operator==()

bool operator== ( const bigintmat lhr,
const bigintmat rhr 
)

Definition at line 159 of file bigintmat.cc.

160 {
161  if (&lhr == &rhr) { return true; }
162  if (lhr.cols() != rhr.cols()) { return false; }
163  if (lhr.rows() != rhr.rows()) { return false; }
164  if (lhr.basecoeffs() != rhr.basecoeffs()) { return false; }
165 
166  const int l = (lhr.rows())*(lhr.cols());
167 
168  for (int i=0; i < l; i++)
169  {
170  if (!n_Equal(lhr[i], rhr[i], lhr.basecoeffs())) { return false; }
171  }
172 
173  return true;
174 }
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff 'a' and 'b' represent the same number; they may have different representations.
Definition: coeffs.h:460

◆ prependIdentity()

static bigintmat* prependIdentity ( bigintmat A)
static

Definition at line 2036 of file bigintmat.cc.

2037 {
2038  coeffs R = A->basecoeffs();
2039  bigintmat *m = new bigintmat(A->rows()+A->cols(), A->cols(), R);
2040  m->copySubmatInto(A, 1, 1, A->rows(), A->cols(), A->cols()+1, 1);
2041  number one = n_Init(1, R);
2042  for(int i=1; i<= A->cols(); i++)
2043  m->set(i,i,one);
2044  n_Delete(&one, R);
2045  return m;
2046 }

◆ reduce_mod_howell()

static void reduce_mod_howell ( bigintmat A,
bigintmat b,
bigintmat eps,
bigintmat x 
)
static

Definition at line 1950 of file bigintmat.cc.

1951 {
1952  //write b = Ax + eps where eps is "small" in the sense of bounded by the
1953  //pivot entries in H. H does not need to be Howell (or HNF) but need
1954  //to be triagonal in the same direction.
1955  //b can have multiple columns.
1956 #if 0
1957  PrintS("reduce_mod_howell: A:\n");
1958  A->Print();
1959  PrintS("\nb:\n");
1960  b->Print();
1961 #endif
1962 
1963  coeffs R = A->basecoeffs();
1964  assume(x->basecoeffs() == R);
1965  assume(b->basecoeffs() == R);
1966  assume(eps->basecoeffs() == R);
1967  if (!A->cols())
1968  {
1969  x->zero();
1970  eps->copy(b);
1971 
1972 #if 0
1973  PrintS("\nx:\n");
1974  x->Print();
1975  PrintS("\neps:\n");
1976  eps->Print();
1977  PrintS("\n****************************************\n");
1978 #endif
1979  return;
1980  }
1981 
1982  bigintmat * B = new bigintmat(b->rows(), 1, R);
1983  for(int i=1; i<= b->cols(); i++)
1984  {
1985  int A_col = A->cols();
1986  b->getcol(i, B);
1987  for(int j = B->rows(); j>0; j--)
1988  {
1989  number Ai = A->view(A->rows() - B->rows() + j, A_col);
1990  if (n_IsZero(Ai, R) &&
1991  n_IsZero(B->view(j, 1), R))
1992  {
1993  continue; //all is fine: 0*x = 0
1994  }
1995  else if (n_IsZero(B->view(j, 1), R))
1996  {
1997  x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
1998  A_col--;
1999  }
2000  else if (n_IsZero(Ai, R))
2001  {
2002  A_col--;
2003  }
2004  else
2005  {
2006  // "solve" ax=b, possibly enlarging d
2007  number Bj = B->view(j, 1);
2008  number q = n_Div(Bj, Ai, R);
2009  x->rawset(x->rows() - B->rows() + j, i, q);
2010  for(int k=j; k>B->rows() - A->rows(); k--)
2011  {
2012  //B[k] = B[k] - x[k]A[k][j]
2013  number s = n_Mult(q, A->view(A->rows() - B->rows() + k, A_col), R);
2014  B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
2015  n_Delete(&s, R);
2016  }
2017  A_col--;
2018  }
2019  if (!A_col)
2020  {
2021  break;
2022  }
2023  }
2024  eps->setcol(i, B);
2025  }
2026  delete B;
2027 #if 0
2028  PrintS("\nx:\n");
2029  x->Print();
2030  PrintS("\neps:\n");
2031  eps->Print();
2032  PrintS("\n****************************************\n");
2033 #endif
2034 }
void setcol(int j, bigintmat *m)
Setzt j-te Spalte gleich übergebenem Vektor (Matrix) m.
Definition: bigintmat.cc:826
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:615
b *CanonicalForm B
Definition: facBivar.cc:52

◆ solveAx()

number solveAx ( bigintmat A,
bigintmat b,
bigintmat x 
)

solve Ax=b*d. x needs to be pre-allocated to the same number of columns as b. the minimal denominator d is returned. Currently available for Z, Q and Z/nZ (and possibly for all fields: d=1 there) Beware that the internal functions can find the kernel as well - but the interface is lacking.

Definition at line 2430 of file bigintmat.cc.

2431 {
2432 #if 0
2433  PrintS("Solve Ax=b for A=\n");
2434  A->Print();
2435  PrintS("\nb = \n");
2436  b->Print();
2437  PrintS("\nx = \n");
2438  x->Print();
2439  PrintLn();
2440 #endif
2441 
2442  coeffs R = A->basecoeffs();
2443  assume (R == b->basecoeffs());
2444  assume (R == x->basecoeffs());
2445  assume ((x->cols() == b->cols()) && (x->rows() == A->cols()) && (A->rows() == b->rows()));
2446 
2447  switch (getCoeffType(R))
2448  {
2449  #ifdef HAVE_RINGS
2450  case n_Z:
2451  return solveAx_dixon(A, b, x, NULL);
2452  case n_Zn:
2453  case n_Znm:
2454  case n_Z2m:
2455  return solveAx_howell(A, b, x, NULL);
2456  #endif
2457  case n_Zp:
2458  case n_Q:
2459  case n_GF:
2460  case n_algExt:
2461  case n_transExt:
2462  WarnS("have field, should use Gauss or better");
2463  break;
2464  default:
2465  if (R->cfXExtGcd && R->cfAnn)
2466  { //assume it's Euclidean
2467  return solveAx_howell(A, b, x, NULL);
2468  }
2469  WerrorS("have no solve algorithm");
2470  break;
2471  }
2472  return NULL;
2473 }
static number solveAx_dixon(bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern)
Definition: bigintmat.cc:2108
static number solveAx_howell(bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
Definition: bigintmat.cc:2298
@ n_GF
\GF{p^n < 2^16}
Definition: coeffs.h:32
@ n_Znm
only used if HAVE_RINGS is defined
Definition: coeffs.h:45
@ n_algExt
used for all algebraic extensions, i.e., the top-most extension in an extension tower is algebraic
Definition: coeffs.h:35
@ n_Z2m
only used if HAVE_RINGS is defined
Definition: coeffs.h:46
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition: coeffs.h:38
#define WarnS
Definition: emacs.cc:78

◆ solveAx_dixon()

static number solveAx_dixon ( bigintmat A,
bigintmat B,
bigintmat x,
bigintmat kern 
)
static

Definition at line 2108 of file bigintmat.cc.

2108  {
2109  coeffs R = A->basecoeffs();
2110 
2111  assume(getCoeffType(R) == n_Z);
2112 
2113  number p = n_Init(536870909, R); // PreviousPrime(2^29); not clever
2114  coeffs Rp = numbercoeffs(p, R); // R/pR
2115  bigintmat *Ap = bimChangeCoeff(A, Rp),
2116  *m = prependIdentity(Ap),
2117  *Tp, *Hp;
2118  delete Ap;
2119 
2120  m->howell();
2121  Hp = new bigintmat(A->rows(), A->cols(), Rp);
2122  Hp->copySubmatInto(m, A->cols()+1, 1, A->rows(), A->cols(), 1, 1);
2123  Tp = new bigintmat(A->cols(), A->cols(), Rp);
2124  Tp->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
2125 
2126  int i, j;
2127 
2128  for(i=1; i<= A->cols(); i++)
2129  {
2130  for(j=m->rows(); j>A->cols(); j--)
2131  {
2132  if (!n_IsZero(m->view(j, i), Rp)) break;
2133  }
2134  if (j>A->cols()) break;
2135  }
2136 // Print("Found nullity (kern dim) of %d\n", i-1);
2137  bigintmat * kp = new bigintmat(A->cols(), i-1, Rp);
2138  kp->copySubmatInto(Tp, 1, 1, A->cols(), i-1, 1, 1);
2139  kp->howell();
2140 
2141  delete m;
2142 
2143  //Hp is the mod-p howell form
2144  //Tp the transformation, mod p
2145  //kp a basis for the kernel, in howell form, mod p
2146 
2147  bigintmat * eps_p = new bigintmat(B->rows(), B->cols(), Rp),
2148  * x_p = new bigintmat(A->cols(), B->cols(), Rp),
2149  * fps_p = new bigintmat(kp->cols(), B->cols(), Rp);
2150 
2151  //initial solution
2152 
2153  number zero = n_Init(0, R);
2154  x->skalmult(zero, R);
2155  n_Delete(&zero, R);
2156 
2157  bigintmat * b = new bigintmat(B);
2158  number pp = n_Init(1, R);
2159  i = 1;
2160  do
2161  {
2162  bigintmat * b_p = bimChangeCoeff(b, Rp), * s;
2163  bigintmat * t1, *t2;
2164  reduce_mod_howell(Hp, b_p, eps_p, x_p);
2165  delete b_p;
2166  if (!eps_p->isZero())
2167  {
2168  PrintS("no solution, since no modular solution\n");
2169 
2170  delete eps_p;
2171  delete x_p;
2172  delete Hp;
2173  delete kp;
2174  delete Tp;
2175  delete b;
2176  n_Delete(&pp, R);
2177  n_Delete(&p, R);
2178  nKillChar(Rp);
2179 
2180  return NULL;
2181  }
2182  t1 = bimMult(Tp, x_p);
2183  delete x_p;
2184  x_p = t1;
2185  reduce_mod_howell(kp, x_p, x_p, fps_p); //we're not all interested in fps_p
2186  s = bimChangeCoeff(x_p, R);
2187  t1 = bimMult(A, s);
2188  t2 = bimSub(b, t1);
2189  t2->skaldiv(p);
2190  delete b;
2191  delete t1;
2192  b = t2;
2193  s->skalmult(pp, R);
2194  t1 = bimAdd(x, s);
2195  delete s;
2196  x->swapMatrix(t1);
2197  delete t1;
2198 
2199  if(kern && i==1)
2200  {
2201  bigintmat * ker = bimChangeCoeff(kp, R);
2202  t1 = bimMult(A, ker);
2203  t1->skaldiv(p);
2204  t1->skalmult(n_Init(-1, R), R);
2205  b->appendCol(t1);
2206  delete t1;
2207  x->appendCol(ker);
2208  delete ker;
2209  x_p->extendCols(kp->cols());
2210  eps_p->extendCols(kp->cols());
2211  fps_p->extendCols(kp->cols());
2212  }
2213 
2214  n_InpMult(pp, p, R);
2215 
2216  if (b->isZero())
2217  {
2218  //exact solution found, stop
2219  delete eps_p;
2220  delete fps_p;
2221  delete x_p;
2222  delete Hp;
2223  delete kp;
2224  delete Tp;
2225  delete b;
2226  n_Delete(&pp, R);
2227  n_Delete(&p, R);
2228  nKillChar(Rp);
2229 
2230  return n_Init(1, R);
2231  }
2232  else
2233  {
2234  bigintmat *y = new bigintmat(x->rows(), x->cols(), R);
2235  number d = bimFarey(x, pp, y);
2236  if (d)
2237  {
2238  bigintmat *c = bimMult(A, y);
2239  bigintmat *bd = new bigintmat(B);
2240  bd->skalmult(d, R);
2241  if (kern)
2242  {
2243  bd->extendCols(kp->cols());
2244  }
2245  if (*c == *bd)
2246  {
2247  x->swapMatrix(y);
2248  delete y;
2249  delete c;
2250  if (kern)
2251  {
2252  y = new bigintmat(x->rows(), B->cols(), R);
2253  c = new bigintmat(x->rows(), kp->cols(), R);
2254  x->splitcol(y, c);
2255  x->swapMatrix(y);
2256  delete y;
2257  kern->swapMatrix(c);
2258  delete c;
2259  }
2260 
2261  delete bd;
2262 
2263  delete eps_p;
2264  delete fps_p;
2265  delete x_p;
2266  delete Hp;
2267  delete kp;
2268  delete Tp;
2269  delete b;
2270  n_Delete(&pp, R);
2271  n_Delete(&p, R);
2272  nKillChar(Rp);
2273 
2274  return d;
2275  }
2276  delete c;
2277  delete bd;
2278  n_Delete(&d, R);
2279  }
2280  delete y;
2281  }
2282  i++;
2283  } while (1);
2284  delete eps_p;
2285  delete fps_p;
2286  delete x_p;
2287  delete Hp;
2288  delete kp;
2289  delete Tp;
2290  n_Delete(&pp, R);
2291  n_Delete(&p, R);
2292  nKillChar(Rp);
2293  return NULL;
2294 }
static void reduce_mod_howell(bigintmat *A, bigintmat *b, bigintmat *eps, bigintmat *x)
Definition: bigintmat.cc:1950
static bigintmat * prependIdentity(bigintmat *A)
Definition: bigintmat.cc:2036
bigintmat * bimSub(bigintmat *a, bigintmat *b)
Definition: bigintmat.cc:218
static number bimFarey(bigintmat *A, number N, bigintmat *L)
Definition: bigintmat.cc:2048
bigintmat * bimAdd(bigintmat *a, bigintmat *b)
Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ? @Note: NULL as a result means an error (non-compati...
Definition: bigintmat.cc:182
CF_NO_INLINE bool isZero() const
void swapMatrix(bigintmat *a)
Definition: bigintmat.cc:1566
int isZero()
Definition: bigintmat.cc:1363
void extendCols(int i)
append i zero-columns to the matrix
Definition: bigintmat.cc:1076
void skaldiv(number b)
Macht Ganzzahldivision aller Matrixeinträge mit b.
Definition: bigintmat.cc:1861
void copySubmatInto(bigintmat *, int sr, int sc, int nr, int nc, int tr, int tc)
copy the submatrix of b, staring at (a,b) having n rows, m cols into the given matrix at pos....
Definition: bigintmat.cc:1287
void howell()
dito, but Howell form (only different for zero-divsors)
Definition: bigintmat.cc:1585
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53

◆ solveAx_howell()

static number solveAx_howell ( bigintmat A,
bigintmat b,
bigintmat x,
bigintmat kern 
)
static

Definition at line 2298 of file bigintmat.cc.

2299 {
2300  // try to solve Ax=b, more precisely, find
2301  // number d
2302  // bigintmat x
2303  // sth. Ax=db
2304  // where d is small-ish (divides the determinant of A if this makes sense)
2305  // return 0 if there is no solution.
2306  //
2307  // if kern is non-NULL, return a basis for the kernel
2308 
2309  //Algo: we do row-howell (triangular matrix). The idea is
2310  // Ax = b <=> AT T^-1x = b
2311  // y := T^-1 x, solve AT y = b
2312  // and return Ty.
2313  //Howell does not compute the trafo, hence we need to cheat:
2314  //B := (I_n | A^t)^t, then the top part of the Howell form of
2315  //B will give a useful trafo
2316  //Then we can find x by back-substitution and lcm/gcd to find the denominator
2317  //The defining property of Howell makes this work.
2318 
2319  coeffs R = A->basecoeffs();
2321  m->howell(); // since m contains the identity, we'll have A->cols()
2322  // many cols.
2323  number den = n_Init(1, R);
2324 
2325  bigintmat * B = new bigintmat(A->rows(), 1, R);
2326  for(int i=1; i<= b->cols(); i++)
2327  {
2328  int A_col = A->cols();
2329  b->getcol(i, B);
2330  B->skalmult(den, R);
2331  for(int j = B->rows(); j>0; j--)
2332  {
2333  number Ai = m->view(m->rows()-B->rows() + j, A_col);
2334  if (n_IsZero(Ai, R) &&
2335  n_IsZero(B->view(j, 1), R))
2336  {
2337  continue; //all is fine: 0*x = 0
2338  }
2339  else if (n_IsZero(B->view(j, 1), R))
2340  {
2341  x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
2342  A_col--;
2343  }
2344  else if (n_IsZero(Ai, R))
2345  {
2346  delete m;
2347  delete B;
2348  n_Delete(&den, R);
2349  return 0;
2350  }
2351  else
2352  {
2353  // solve ax=db, possibly enlarging d
2354  // so x = db/a
2355  number Bj = B->view(j, 1);
2356  number g = n_Gcd(Bj, Ai, R);
2357  number xi;
2358  if (n_Equal(Ai, g, R))
2359  { //good: den stable!
2360  xi = n_Div(Bj, Ai, R);
2361  }
2362  else
2363  { //den <- den * (a/g), so old sol. needs to be adjusted
2364  number inc_d = n_Div(Ai, g, R);
2365  n_InpMult(den, inc_d, R);
2366  x->skalmult(inc_d, R);
2367  B->skalmult(inc_d, R);
2368  xi = n_Div(Bj, g, R);
2369  n_Delete(&inc_d, R);
2370  } //now for the back-substitution:
2371  x->rawset(x->rows() - B->rows() + j, i, xi);
2372  for(int k=j; k>0; k--)
2373  {
2374  //B[k] = B[k] - x[k]A[k][j]
2375  number s = n_Mult(xi, m->view(m->rows()-B->rows() + k, A_col), R);
2376  B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
2377  n_Delete(&s, R);
2378  }
2379  n_Delete(&g, R);
2380  A_col--;
2381  }
2382  if (!A_col)
2383  {
2384  if (B->isZero()) break;
2385  else
2386  {
2387  delete m;
2388  delete B;
2389  n_Delete(&den, R);
2390  return 0;
2391  }
2392  }
2393  }
2394  }
2395  delete B;
2396  bigintmat *T = new bigintmat(A->cols(), A->cols(), R);
2397  T->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
2398  if (kern)
2399  {
2400  int i, j;
2401  for(i=1; i<= A->cols(); i++)
2402  {
2403  for(j=m->rows(); j>A->cols(); j--)
2404  {
2405  if (!n_IsZero(m->view(j, i), R)) break;
2406  }
2407  if (j>A->cols()) break;
2408  }
2409  Print("Found nullity (kern dim) of %d\n", i-1);
2410  bigintmat * ker = new bigintmat(A->rows(), i-1, R);
2411  ker->copySubmatInto(T, 1, 1, A->rows(), i-1, 1, 1);
2412  kern->swapMatrix(ker);
2413  delete ker;
2414  }
2415  delete m;
2416  bigintmat * y = bimMult(T, x);
2417  x->swapMatrix(y);
2418  delete y;
2419  x->simplifyContentDen(&den);
2420 #if 0
2421  PrintS("sol = 1/");
2422  n_Print(den, R);
2423  PrintS(" *\n");
2424  x->Print();
2425  PrintLn();
2426 #endif
2427  return den;
2428 }
g
Definition: cfModGcd.cc:4090
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ,...
Definition: coeffs.h:664