My Project
Loading...
Searching...
No Matches
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
 
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?)
 
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
 
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.
 
void bimMult (bigintmat *a, bigintmat *b, bigintmat *c)
 Multipliziert Matrix a und b und speichert Ergebnis in c.
 
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.
 
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.
 
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:4111
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:551

◆ 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:654
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:459
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:542

◆ bimChangeCoeff()

bigintmat * bimChangeCoeff ( bigintmat a,
coeffs  cnew 
)

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

Definition at line 1800 of file bigintmat.cc.

1801{
1802 coeffs cold = a->basecoeffs();
1803 bigintmat *b = new bigintmat(a->rows(), a->cols(), cnew);
1804 // Erzeugt Karte von alten coeffs nach neuen
1806 number t1;
1807 number t2;
1808 // apply map to all entries.
1809 for (int i=1; i<=a->rows(); i++)
1810 {
1811 for (int j=1; j<=a->cols(); j++)
1812 {
1813 t1 = a->get(i, j);
1814 t2 = f(t1, cold, cnew);
1815 b->set(i, j, t2);
1816 n_Delete(&t1, cold);
1817 n_Delete(&t2, cnew);
1818 }
1819 }
1820 return b;
1821}
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:704
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition coeffs.h:80
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 2044 of file bigintmat.cc.

2045{
2046 coeffs Z = A->basecoeffs(),
2047 Q = nInitChar(n_Q, 0);
2048 number den = n_Init(1, Z);
2049 nMapFunc f = n_SetMap(Q, Z);
2050
2051 for(int i=1; i<= A->rows(); i++)
2052 {
2053 for(int j=1; j<= A->cols(); j++)
2054 {
2055 number ad = n_Mult(den, A->view(i, j), Z);
2056 number re = n_IntMod(ad, N, Z);
2057 n_Delete(&ad, Z);
2058 number q = n_Farey(re, N, Z);
2059 n_Delete(&re, Z);
2060 if (!q)
2061 {
2062 n_Delete(&ad, Z);
2063 n_Delete(&den, Z);
2064 return NULL;
2065 }
2066
2067 number d = n_GetDenom(q, Q),
2068 n = n_GetNumerator(q, Q);
2069
2070 n_Delete(&q, Q);
2071 n_Delete(&ad, Z);
2072 number dz = f(d, Q, Z),
2073 nz = f(n, Q, Z);
2074 n_Delete(&d, Q);
2075 n_Delete(&n, Q);
2076
2077 if (!n_IsOne(dz, Z))
2078 {
2079 L->skalmult(dz, Z);
2080 n_InpMult(den, dz, Z);
2081#if 0
2082 PrintS("den increasing to ");
2083 n_Print(den, Z);
2084 PrintLn();
2085#endif
2086 }
2087 n_Delete(&dz, Z);
2088 L->rawset(i, j, nz);
2089 }
2090 }
2091
2092 nKillChar(Q);
2093 PrintS("bimFarey worked\n");
2094#if 0
2095 L->Print();
2096 PrintS("\n * 1/");
2097 n_Print(den, Z);
2098 PrintLn();
2099#endif
2100 return den;
2101}
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:439
bool skalmult(number b, coeffs c)
Multipliziert zur Matrix den Skalar b hinzu.
Definition bigintmat.cc:934
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:640
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:607
@ 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:673
static FORCE_INLINE number n_Farey(number a, number b, const coeffs r)
Definition coeffs.h:771
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition numbers.cc:419
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:632
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:645
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:612
void nKillChar(coeffs r)
undo all initialisations
Definition numbers.cc:574
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition coeffs.h:472
void PrintS(const char *s)
Definition reporter.cc:284
void PrintLn()
Definition reporter.cc:310
#define A
Definition sirandom.c:24
#define Q
Definition sirandom.c:26

◆ 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:650
fq_nmod_poly_t prod
Definition facHensel.cc:100
#define assume(x)
Definition mod2.h:387
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 1928 of file bigintmat.cc.

1929{
1930 if (!nCoeffs_are_equal(a->basecoeffs(), b->basecoeffs()))
1931 {
1932 WerrorS("Error in bimMult. Coeffs do not agree!");
1933 return;
1934 }
1935 if ((a->rows() != c->rows()) || (b->cols() != c->cols()) || (a->cols() != b->rows()))
1936 {
1937 WerrorS("Error in bimMult. Dimensions do not agree!");
1938 return;
1939 }
1940 bigintmat *tmp = bimMult(a, b);
1941 c->copy(tmp);
1942
1943 delete tmp;
1944}
bool nCoeffs_are_equal(coeffs r, coeffs s)
bigintmat * bimMult(bigintmat *a, bigintmat *b)
Definition bigintmat.cc:255
bool copy(bigintmat *b)
Kopiert Einträge von b auf Bigintmat.
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:4091

◆ 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:659

◆ 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 2471 of file bigintmat.cc.

2472{
2473 bigintmat * t, *s, *a=A;
2474 coeffs R = a->basecoeffs();
2475 if (T)
2476 {
2477 *T = new bigintmat(a->cols(), a->cols(), R),
2478 (*T)->one();
2479 t = new bigintmat(*T);
2480 }
2481 else
2482 {
2483 t = *T;
2484 }
2485
2486 if (S)
2487 {
2488 *S = new bigintmat(a->rows(), a->rows(), R);
2489 (*S)->one();
2490 s = new bigintmat(*S);
2491 }
2492 else
2493 {
2494 s = *S;
2495 }
2496
2497 int flip=0;
2498 do
2499 {
2500 bigintmat * x, *X;
2501 if (flip)
2502 {
2503 x = s;
2504 X = *S;
2505 }
2506 else
2507 {
2508 x = t;
2509 X = *T;
2510 }
2511
2512 if (x)
2513 {
2514 x->one();
2515 bigintmat * r = new bigintmat(a->rows()+a->cols(), a->cols(), R);
2516 bigintmat * rw = new bigintmat(1, a->cols(), R);
2517 for(int i=0; i<a->cols(); i++)
2518 {
2519 x->getrow(i+1, rw);
2520 r->setrow(i+1, rw);
2521 }
2522 for (int i=0; i<a->rows(); i++)
2523 {
2524 a->getrow(i+1, rw);
2525 r->setrow(i+a->cols()+1, rw);
2526 }
2527 r->hnf();
2528 for(int i=0; i<a->cols(); i++)
2529 {
2530 r->getrow(i+1, rw);
2531 x->setrow(i+1, rw);
2532 }
2533 for(int i=0; i<a->rows(); i++)
2534 {
2535 r->getrow(i+a->cols()+1, rw);
2536 a->setrow(i+1, rw);
2537 }
2538 delete rw;
2539 delete r;
2540
2541#if 0
2542 Print("X: %ld\n", X);
2543 X->Print();
2544 Print("\nx: %ld\n", x);
2545 x->Print();
2546#endif
2547 bimMult(X, x, X);
2548#if 0
2549 Print("\n2:X: %ld %ld %ld\n", X, *S, *T);
2550 X->Print();
2551 Print("\n2:x: %ld\n", x);
2552 x->Print();
2553 PrintLn();
2554#endif
2555 }
2556 else
2557 {
2558 a->hnf();
2559 }
2560
2561 int diag = 1;
2562 for(int i=a->rows(); diag && i>0; i--)
2563 {
2564 for(int j=a->cols(); j>0; j--)
2565 {
2566 if ((a->rows()-i)!=(a->cols()-j) && !n_IsZero(a->view(i, j), R))
2567 {
2568 diag = 0;
2569 break;
2570 }
2571 }
2572 }
2573#if 0
2574 PrintS("Diag ? %d\n", diag);
2575 a->Print();
2576 PrintLn();
2577#endif
2578 if (diag) break;
2579
2580 a = a->transpose(); // leaks - I need to write inpTranspose
2581 flip = 1-flip;
2582 } while (1);
2583 if (flip)
2584 a = a->transpose();
2585
2586 if (S) *S = (*S)->transpose();
2587 if (s) delete s;
2588 if (t) delete t;
2589 A->copy(a);
2590}
Variable x
Definition cfModGcd.cc:4090
void hnf()
transforms INPLACE to HNF
bigintmat * transpose()
Definition bigintmat.cc:37
void setrow(int i, bigintmat *m)
Setzt i-te Zeile gleich übergebenem Vektor (Matrix) m.
Definition bigintmat.cc:856
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.
void getrow(int i, bigintmat *a)
Schreibt i-te Zeile in Vektor (Matrix) a.
Definition bigintmat.cc:787
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition coeffs.h:468
#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 533 of file bigintmat.cc.

534{
535 int l = 0;
536 int index;
537 for (int i=0; i<length; i++)
538 {
539 if (a[i] > l)
540 {
541 l = a[i];
542 index = i;
543 }
544 }
545 return index;
546}
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)

◆ getShorter()

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

Definition at line 548 of file bigintmat.cc.

549{
550 int sndlong = 0;
551 int min;
552 for (int i=0; i<rows; i++)
553 {
554 int index = cols*i+j;
555 if ((a[index] > sndlong) && (a[index] < l))
556 {
557 min = floor(log10((double)cols))+floor(log10((double)rows))+5;
558 if ((a[index] < min) && (min < l))
559 sndlong = min;
560 else
561 sndlong = a[index];
562 }
563 }
564 if (sndlong == 0)
565 {
566 min = floor(log10((double)cols))+floor(log10((double)rows))+5;
567 if (min < l)
568 sndlong = min;
569 else
570 sndlong = 1;
571 }
572 return sndlong;
573}
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 525 of file bigintmat.cc.

526{
527 int sum = 0;
528 for (int i=0; i<length; i++)
529 sum += a[i];
530 return sum;
531}

◆ 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 2596 of file bigintmat.cc.

2597{
2598#if 0
2599 PrintS("Kernel of ");
2600 a->Print();
2601 PrintS(" modulo ");
2602 n_Print(p, q);
2603 PrintLn();
2604#endif
2605
2606 coeffs coe = numbercoeffs(p, q);
2607 bigintmat *m = bimChangeCoeff(a, coe), *U, *V;
2608 diagonalForm(m, &U, &V);
2609#if 0
2610 PrintS("\ndiag form: ");
2611 m->Print();
2612 PrintS("\nU:\n");
2613 U->Print();
2614 PrintS("\nV:\n");
2615 V->Print();
2616 PrintLn();
2617#endif
2618
2619 int rg = 0;
2620#undef MIN
2621#define MIN(a,b) (a < b ? a : b)
2622 for(rg=0; rg<MIN(m->rows(), m->cols()) && !n_IsZero(m->view(m->rows()-rg,m->cols()-rg), coe); rg++);
2623
2624 bigintmat * k = new bigintmat(m->cols(), m->rows(), coe);
2625 for(int i=0; i<rg; i++)
2626 {
2627 number A = n_Ann(m->view(m->rows()-i, m->cols()-i), coe);
2628 k->set(m->cols()-i, i+1, A);
2629 n_Delete(&A, coe);
2630 }
2631 for(int i=rg; i<m->cols(); i++)
2632 {
2633 k->set(m->cols()-i, i+1-rg, n_Init(1, coe));
2634 }
2635 bimMult(V, k, k);
2636 c->copy(bimChangeCoeff(k, q));
2637 return c->cols();
2638}
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.
#define MIN(a, b)
void diagonalForm(bigintmat *A, bigintmat **S, bigintmat **T)
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:4086
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:683

◆ nCoeffs_are_equal()

bool nCoeffs_are_equal ( coeffs  r,
coeffs  s 
)

Definition at line 2641 of file bigintmat.cc.

2642{
2643 if ((r == NULL) || (s == NULL))
2644 return false;
2645 if (r == s)
2646 return true;
2647 if ((getCoeffType(r)==n_Z) && (getCoeffType(s)==n_Z))
2648 return true;
2649 if ((getCoeffType(r)==n_Zp) && (getCoeffType(s)==n_Zp))
2650 {
2651 if (r->ch == s->ch)
2652 return true;
2653 else
2654 return false;
2655 }
2656 // n_Zn stimmt wahrscheinlich noch nicht
2657 if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2658 {
2659 if (r->ch == s->ch)
2660 return true;
2661 else
2662 return false;
2663 }
2664 if ((getCoeffType(r)==n_Q) && (getCoeffType(s)==n_Q))
2665 return true;
2666 // FALL n_Zn FEHLT NOCH!
2667 //if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2668 return false;
2669}
@ 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:429

◆ 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:991

◆ 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:464

◆ prependIdentity()

static bigintmat * prependIdentity ( bigintmat A)
static

Definition at line 2032 of file bigintmat.cc.

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

◆ reduce_mod_howell()

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

Definition at line 1946 of file bigintmat.cc.

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

2427{
2428#if 0
2429 PrintS("Solve Ax=b for A=\n");
2430 A->Print();
2431 PrintS("\nb = \n");
2432 b->Print();
2433 PrintS("\nx = \n");
2434 x->Print();
2435 PrintLn();
2436#endif
2437
2438 coeffs R = A->basecoeffs();
2439 assume (R == b->basecoeffs());
2440 assume (R == x->basecoeffs());
2441 assume ((x->cols() == b->cols()) && (x->rows() == A->cols()) && (A->rows() == b->rows()));
2442
2443 switch (getCoeffType(R))
2444 {
2445 #ifdef HAVE_RINGS
2446 case n_Z:
2447 return solveAx_dixon(A, b, x, NULL);
2448 case n_Zn:
2449 case n_Znm:
2450 case n_Z2m:
2451 return solveAx_howell(A, b, x, NULL);
2452 #endif
2453 case n_Zp:
2454 case n_Q:
2455 case n_GF:
2456 case n_algExt:
2457 case n_transExt:
2458 WarnS("have field, should use Gauss or better");
2459 break;
2460 default:
2461 if (R->cfXExtGcd && R->cfAnn)
2462 { //assume it's Euclidean
2463 return solveAx_howell(A, b, x, NULL);
2464 }
2465 WerrorS("have no solve algorithm");
2466 break;
2467 }
2468 return NULL;
2469}
static number solveAx_dixon(bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern)
static number solveAx_howell(bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
@ 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 2104 of file bigintmat.cc.

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

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