Ramón Casero Cañas on Thu, 16 Jan 2003 09:37:16 +0100 (MET)


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

thank you for the inversion


Thank you all for the help. At last my function takes a C++ matrix, 
inverts it mudulo2, and overwrites the input matrix with the result. The 
code is below, just in case you want to have a look at this PARI's 
`Hello world'. Of course you're going to say `it's so easy, and 
everything that's hard was already done'. Well, it's been very 
interesting getting to know PARI/GP, and I hope using it further.

I have also tried to do it the u_FpM_inv() way, but I'm still struggling 
trying to get that it compiles. I have some directions from some of you 
about this problem, so I'll try to overcome it.

Best,

Ramón.

--------------------------------------------------

bool inverse_mat_mod2(bool **_a, unsigned _size) {
   return inverse_mat_mod2(_a, _size, true);
}

bool inverse_mat_mod2(bool **_a, const unsigned _size, const bool _CHECK) {
   GEN a;           // temp matrix by columns
   unsigned i, j;   // counter
   unsigned sizepari = _size + 1;
   GEN z;           // auxiliary element of the matrix
   GEN v;           // auxiliary element of the matrix
   long paristacksize = 500000;
   long total_size = sizepari * sizepari + _size * _size * 9; // thanks, 
Karim

   paristacksize = sizeof(long) * total_size * 10; // thanks again
   pari_init(paristacksize, 2); // init PARI stack

   /* allocate memory for the temp matrix
    * note: remember that PARI matrices of mxn need (m+1)x(n+1), and
    * that elements go a[1], a[2], ..., while a[0] is for type
    * information */
   std::cout << "inverse_mat_mod2: size = " << _size << std::endl;
   a = cgetg(sizepari, t_MAT); // `a' has _size columns
   for (i = 1; i <= _size; ++i) { // allocate `a' columns of _size elements
     a[i] = lgetg(sizepari, t_COL);
   }
   for (i = 1; i <= _size; ++i) { // allocate space for elements
     v = (GEN)a[i];
     for (j = 1; j <= _size; ++j) {
       v[j] = lgetg(3, t_INTMOD);
       z = (GEN)v[j];
       z[1] = lstoi(2); // modulo 2
       if (_a[j-1][i-1]) z[2] = un; else z[2] = zero;//lstoi(0);
     }
   }
   std::cout << "inverse_mat_mod2: assigment done" << std::endl;

   /* invert matrix */
   GEN ainv;
   ainv = ginv(a);

   /* check whether the inverse is valid */
   if (_CHECK) {
     GEN a_ainv;
     a_ainv = gmul(a, ainv);
     for (i = 1; i <= _size; ++i) {
       v = (GEN)a_ainv[i];
       for (j = 1; j <= _size; ++j) {
         z = (GEN)v[j];
         if ((i == j) != itos((GEN)z[2])) return false;
       }
     }
   }

   /* for no _CHECK or valid inversion, overwrite input array with
      inverse matrix */
   for (i = 1; i <= _size; ++i) {
     v = (GEN)ainv[i];
     for (j = 1; j <= _size; ++j) {
       z = (GEN)v[j];
       _a[j-1][i-1] = (bool)itos((GEN)z[2]);
     }
   }

   return true;
}

-- 
+++ GMX - Mail, Messaging & more  http://www.gmx.net +++
NEU: Mit GMX ins Internet. Rund um die Uhr für 1 ct/ Min. surfen!