1 #include "vectorsF2.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4 
5 #define WL vectorsF2_WL
6 
7 #define MC 0x80000000UL    /* permet de diagonaliser la matrice dans Diag() */
8 
9 unsigned long MMC[WL] =
10    { MC, MC >> 1, MC >> 2, MC >> 3, MC >> 4, MC >> 5, MC >> 6, MC >> 7,
11    MC >> 8, MC >> 9, MC >> 10,
12    MC >> 11, MC >> 12, MC >> 13, MC >> 14, MC >> 15, MC >> 16, MC >> 17,
13    MC >> 18, MC >> 19, MC >> 20,
14    MC >> 21, MC >> 22, MC >> 23, MC >> 24, MC >> 25, MC >> 26, MC >> 27,
15    MC >> 28, MC >> 29, MC >> 30, MC >> 31
16 };
17 
InverseMatrix(Matrix * InvM,Matrix * M)18 lebool InverseMatrix (Matrix * InvM, Matrix * M)
19 {
20 
21    Matrix Temp;
22    int j, rang;
23    if (M->nblignes != M->l) {
24       printf ("Matrix M is not square!\n");
25       exit (1);
26    }
27    AllocMat (&Temp, M->nblignes, M->l, 2);
28    for (j = 0; j < M->l; j++)
29       CopyBV (&(Temp.lignes[j][0]), &(M->lignes[j][0]));
30    for (j = 0; j < M->l; j++) {
31       BVCanonic (&(Temp.lignes[j][1]), j);
32    }
33    /* DispMat(&Temp,2,M->l,M->nblignes,0); */
34    rang = CompleteElimination (&Temp, M->nblignes, M->l, 2);
35    /* DispMat(&Temp,2,M->l,M->nblignes,0); */
36    /* printf("rang=%d",rang); */
37    for (j = 0; j < M->l; j++)
38       CopyBV (&(InvM->lignes[j][0]), &(Temp.lignes[j][1]));
39    return (rang == M->l);
40    FreeMat (&Temp);
41 }
42 
43 
44 /* ********************************************************************** */
45 /* lebool Diag( Matrix m, int kg,				          */
46 /*               int t, int l, int *gr )                                  */
47 /* Evalue si la matrice de travail m sur kg lignes est de plein rang en   */
48 /* la diagonalisant.  On procede sur t BitVect en considerant les l       */
49 /* premiers bits de chacun.  La fonction retourne TRUE si la matrice m    */
50 /* est de plein rang t*l et *gr est inchange. La fonction retourne        */
51 /* FALSE sinon et *gr prend pour valeur le numero du BitVect ou il y a    */
52 /* eu echec moins un ( = dimension pour laquelle on a resolution l ).     */
53 /* ********************************************************************** */
Diag(Matrix * m,int kg,int t,int l,int * gr)54 lebool Diag (Matrix * m, int kg, int t, int l, int *gr)
55 {
56    int i, j, cl, rang;
57 
58    rang = 0;
59 
60    /* On diagonalise la matrice des entrees (i,j) sur l bits */
61    /* avec 0 <= i < kg et 0 <= j < t .  */
62 
63    for (j = 0; j < t; j++) {
64       cl = 1;
65       while (cl <= l) {
66          /* On cherche dans la j-eme colonne, commencant a la */
67          /* rang-eme ligne, la premiere entree dont le bit le plus */
68          /* significatif est non nul.  Bref, on cherche un pivot.  */
69          i = rang;
70 
71          while ((i < kg)
72             && (m->lignes[i][j].vect[(cl - 1) / WL] < MMC[(cl - 1) % WL]))
73             i++;
74          if (i < kg) {            /* pivot trouve ... */
75             ExchangeVect (m, rang, i);
76             for (i = rang + 1; i < kg; i++) {
77                if (m->lignes[i][j].vect[(cl - 1) / WL] & MMC[(cl - 1) % WL])
78                   XorVect (m, i, rang, j, m->t);
79             }
80             rang++;
81          } else {                 /* pas de pivot trouve ... ==> pas de plein
82                                      rang ... */
83 
84             *gr = j;              /* no de groupe ou il y a echec moins un */
85             return FALSE;         /* c'est j car on indexe a partir de 0 .  */
86          }
87          cl++;
88       }
89    }
90    return TRUE;                   /* on a trouve tous les pivots ==> plein
91                                      rang ! */
92 }
93 
94 
CompleteElimination(Matrix * m,int nblignes,int l,int t)95 int CompleteElimination (Matrix * m, int nblignes, int l, int t)
96 {
97    int i, j, cl, rang;
98 
99    rang = 0;
100 
101    j = 0;
102    while (j < t) {
103       cl = 0;
104       while (cl < l) {
105          /* On cherche dans la j-eme colonne, commencant a la */
106          /* rang-eme ligne, la premiere entree dont le bit le plus */
107          /* significatif est non nul.  Bref, on cherche un pivot.  */
108          i = rang;
109          while ((i < nblignes)
110             && ((((m->lignes)[i])[j]).vect[(cl) / WL] < MMC[(cl) % WL]))
111             i++;
112          if (i < nblignes) {      /* pivot trouve ... */
113             ExchangeVect (m, rang, i);
114             for (i = 0; i < nblignes; i++)
115                if (i != rang)
116                   if ((((m->lignes)[i])[j]).vect[(cl) / WL] & MMC[(cl) % WL])
117                      XorVect (m, i, rang, 0, m->t);
118 
119             rang++;
120             if (rang == nblignes)
121                return rang;
122          } else
123             return rang;
124          cl++;
125       }
126       j++;
127    }
128    return rang;
129 }
130 
131 
GaussianElimination(Matrix * m,int nblignes,int l,int t)132 int GaussianElimination (Matrix * m, int nblignes, int l, int t)
133 {
134    int i, j, cl, rang;
135 
136    rang = 0;
137 
138    j = 0;
139    while (j < t) {
140       cl = 0;
141       while (cl < l) {
142          /* On cherche dans la j-eme colonne, commencant a la */
143          /* rang-eme ligne, la premiere entree dont le bit le plus */
144          /* significatif est non nul.  Bref, on cherche un pivot.  */
145          i = rang;
146          while ((i < nblignes)
147             && ((((m->lignes)[i])[j]).vect[(cl) / WL] < MMC[(cl) % WL]))
148             i++;
149          if (i < nblignes) {      /* pivot trouve ... */
150             ExchangeVect (m, rang, i);
151             for (i = rang + 1; i < nblignes; i++)
152                if ((((m->lignes)[i])[j]).vect[(cl) / WL] & MMC[(cl) % WL])
153                   XorVect (m, i, rang, 0, m->t);
154 
155             rang++;
156             if (rang == nblignes)
157                return rang;
158          }
159          cl++;
160       }
161       j++;
162    }
163    return rang;
164 }
165 
166 
SpecialGaussianElimination(Matrix * m,int nblignes,int l,int t,int * indices)167 int SpecialGaussianElimination (Matrix * m, int nblignes, int l, int t,
168    int *indices)
169 {
170    int i, j, cl, rang;
171 
172    rang = 0;
173 
174    j = 0;
175    while (j < t) {
176       cl = 0;
177       while (cl < l) {
178          /* On cherche dans la j-eme colonne, commencant a la */
179          /* rang-eme ligne, la premiere entree dont le bit le plus */
180          /* significatif est non nul.  Bref, on cherche un pivot.  */
181          i = rang;
182          while ((i < nblignes)
183             && ((((m->lignes)[i])[indices[j]]).vect[(cl) / WL] <
184                MMC[(cl) % WL]))
185             i++;
186          if (i < nblignes) {      /* pivot trouve ... */
187             ExchangeVect (m, rang, i);
188             for (i = rang + 1; i < nblignes; i++)
189                if ((((m->lignes)[i])[indices[j]]).vect[(cl) /
190                      WL] & MMC[(cl) % WL])
191                   XorVect (m, i, rang, 0, m->t);
192 
193             rang++;
194             if (rang == nblignes)
195                return rang;
196          }
197          cl++;
198       }
199       j++;
200    }
201    return rang;
202 }
203 
204 
MultMatrixByBV(BitVect * A,Matrix * M,BitVect * B)205 void MultMatrixByBV (BitVect * A, Matrix * M, BitVect * B)
206 {
207    int i, j, res;
208    if (M->l < B->n * WL) {
209       printf ("Error in MultMatrixByBV(): sizes do not match\n");
210       exit (1);
211    }
212    if (A->n * WL < M->nblignes) {
213       printf ("Error in MultMatrixByBV(): sizes do not match\n");
214       exit (1);
215    }
216    if (M->t != 1) {
217       printf ("Error in MultMatrixByBV(): Not implemented for M->t > 1\n");
218       exit (1);
219    }
220    PutBVToZero (A);
221    for (i = 0; i < M->nblignes; i++) {
222       res = 0;
223       for (j = 0; j < M->l; j++)
224          res += ValBitBV (&(M->lignes[i][0]), j) * ValBitBV (B, j);
225       res %= 2;
226       PutBitBV (A, i, res);
227    }
228 }
229 
230 
231 /* ********************************************************************** */
232 /* int GetBitBV (BitVect A, int noBit)			  */
233 /* Fonction qui permet de prendre la valeur du noBit-ieme bit             */
234 /* (l'indexation commence au bit 0)					  */
235 /* ********************************************************************** */
ValBitBV(BitVect * A,int noBit)236 int ValBitBV (BitVect * A, int noBit)
237 {
238    int k;
239    unsigned long mask;
240    k = noBit / WL;
241    mask = 0x80000000UL >> (noBit - k * WL);
242    if (A->vect[k] & mask)
243       return 1;
244    else
245       return 0;
246 }
247 
248 
249 /* ********************************************************************** */
250 /* void SetBitBV (BitVect A, int noBit, int valBit)		  */
251 /* Fonction qui permet de rendre la valeur du noBit-ieme bit `a valBit     */
252 /* (l'indexation commence au bit 0) (valBit=1 ou valBit=0)		  */
253 /* ********************************************************************** */
PutBitBV(BitVect * A,int noBit,int valBit)254 void PutBitBV (BitVect * A, int noBit, int valBit)
255 {
256    int k;
257    unsigned long mask;
258 
259    k = noBit / WL;
260    if (valBit == 1) {
261       mask = 0x80000000UL >> (noBit - k * WL);
262       A->vect[k] |= mask;
263    } else {
264       mask = 0xffffffffUL ^ (0x80000000UL >> (noBit - k * WL));
265       A->vect[k] &= mask;
266    }
267 }
268 
269 
270 /* ********************************************************************** */
271 /* SetBVToZero( BitVect A)        					  */
272 /* Initialise le vecteur de bit a zero       				  */
273 /* ********************************************************************** */
PutBVToZero(BitVect * A)274 void PutBVToZero (BitVect * A)
275 {
276    int i;
277    for (i = 0; i < A->n; i++)
278       A->vect[i] = 0UL;
279 }
280 
281 
282 /* ********************************************************************** */
283 /* void CopyBV(BitVect A, BitVect B)       				  */
284 /* Copie le contenu de B dans A (A=B)       				  */
285 /* ********************************************************************** */
CopyBV(BitVect * A,BitVect * B)286 void CopyBV (BitVect * A, BitVect * B)
287 {
288    int i;
289 
290    if (A->n != B->n) {
291       printf
292     ("Error in CopyBV(): vectors of different dimensions! (%d and %d bits)\n",
293          A->n * WL, B->n * WL);
294       exit (1);
295    }
296 
297    if (B->n == 0) {
298       printf ("Nothing to copy!\n");
299       exit (1);
300    }
301    for (i = 0; i < B->n; i++)
302       A->vect[i] = B->vect[i];
303 }
CopyBVPart(BitVect * A,BitVect * B,int l)304 void CopyBVPart (BitVect * A, BitVect * B, int l)
305 {
306 
307    int i, n;
308 
309    n = (l - 1) / WL + 1;
310 
311    if (A->n < n) {
312       printf ("Error in CopyBVPart() : The vector A is not large enough!\n");
313       exit (1);
314    }
315    if (B->n == 0) {
316       printf ("Nothing to copy!\n");
317       exit (1);
318    }
319 
320    for (i = 0; i < n; i++)
321       A->vect[i] = B->vect[i];
322 
323    if (l % WL) {
324       BitVect m;
325       AllocBV (&m, A->n * WL);
326       Mask (&m, l);
327       ANDBVSelf (A, &m);
328       FreeBV (&m);
329    }
330 }
331 
332 
333 /* ********************************************************************** */
334 /* void EgalBV(BitVect A, BitVect B)       				  */
335 /* Compare le contenu de B avec celui de A.  Retourne TRUE si les deux    */
336 /* contiennent la meme information.       				  */
337 /* ********************************************************************** */
CompareBV(BitVect * A,BitVect * B)338 lebool CompareBV (BitVect * A, BitVect * B)
339 {
340 
341    int i;
342 
343    if (A->n != B->n) {
344       printf ("Error in EgalBV(): Vectors of different sizes\n");
345       exit (1);
346    }
347 
348    for (i = 0; i < A->n; i++)
349       if (A->vect[i] != B->vect[i])
350          return FALSE;
351    return TRUE;
352 }
353 
354 
BVisZero(BitVect * A)355 lebool BVisZero (BitVect * A)
356 {
357    int j = 0;
358    while (j < A->n)
359       if (A->vect[j++] != 0UL)
360          return FALSE;
361    return TRUE;
362 }
363 
364 
365 /* ********************************************************************** */
366 /* void XORBV(BitVect A, BitVect B, BitVect C)      			  */
367 /* Cette fonction effectue A = B ^ C       				  */
368 /* ********************************************************************** */
XORBV(BitVect * A,BitVect * B,BitVect * C)369 void XORBV (BitVect * A, BitVect * B, BitVect * C)
370 {
371 
372    int i;
373 
374    if ((A->n != B->n) || (B->n != C->n)) {
375       printf ("Error in XORBV(): Vectors of different sizes\n");
376       exit (1);
377    }
378 
379    for (i = 0; i < B->n; i++)
380       A->vect[i] = B->vect[i] ^ C->vect[i];
381 }
382 
383 
384 /* ********************************************************************** */
385 /* void XOR2BV(BitVect A, BitVect B, BitVect C, BitVect D)    		  */
386 /* Cette fonction effectue A = B ^ C ^ D      				  */
387 /* ********************************************************************** */
XOR2BV(BitVect * A,BitVect * B,BitVect * C,BitVect * D)388 void XOR2BV (BitVect * A, BitVect * B, BitVect * C, BitVect * D)
389 {
390 
391    int i;
392 
393    if ((A->n != B->n) || (B->n != C->n) || (C->n != D->n)) {
394       printf ("Error in XOR2BV(): Vectors of different sizes\n");
395       exit (1);
396    }
397 
398    for (i = 0; i < B->n; i++)
399       A->vect[i] = B->vect[i] ^ C->vect[i] ^ D->vect[i];
400 }
401 
402 
403 /* ********************************************************************** */
404 /* void ANDBV(BitVect A, BitVect B, BitVect C)      			  */
405 /* Cette fonction effectue A = B & C       				  */
406 /* ********************************************************************** */
ANDBV(BitVect * A,BitVect * B,BitVect * C)407 void ANDBV (BitVect * A, BitVect * B, BitVect * C)
408 {
409 
410    int i;
411 
412    if ((A->n != B->n) || (B->n != C->n)) {
413       printf ("Error in ANDBV(): Vectors of different sizes\n");
414       exit (1);
415    }
416 
417    for (i = 0; i < B->n; i++)
418       A->vect[i] = B->vect[i] & C->vect[i];
419 }
420 
421 
ANDBVMask(BitVect * A,BitVect * B,int t)422 void ANDBVMask (BitVect * A, BitVect * B, int t)
423 {
424    int n, m, j;
425    if (A->n != B->n) {
426       printf ("Error in ANDBVMask(): Vectors of different sizes\n");
427       exit (1);
428    }
429 
430    if (t > B->n * WL)
431       CopyBV (A, B);
432    else if (t == 0)
433       PutBVToZero (A);
434    else {
435       n = t / WL;
436       m = t - n * WL;
437       for (j = 0; j < n; j++) {
438          A->vect[j] = B->vect[j];
439 
440       }
441       if (m != 0) {
442          A->vect[j] = B->vect[j] & (0xffffffffUL << (WL - m));
443          j++;
444       }
445       /* printf("n=%d j=%d %d m=%d ",n,j,A->n,m); */
446       for (; j < A->n; j++) {
447          A->vect[j] = 0UL;
448       }
449    }
450 }
451 
452 
ANDBVInvMask(BitVect * A,BitVect * B,int t)453 void ANDBVInvMask (BitVect * A, BitVect * B, int t)
454 {
455    int n, m, j;
456    if (A->n != B->n) {
457       printf ("Error in ANDBV(): Vectors of different sizes\n");
458       exit (1);
459    }
460    if (t > B->n * WL)
461       PutBVToZero (A);
462    else if (t == 0)
463       CopyBV (A, B);
464    else {
465       n = t / WL;
466       m = t - n * WL;
467       for (j = 0; j < n; j++)
468          A->vect[j] = 0UL;
469       if (m == 0)
470          A->vect[j] = B->vect[j];
471       else {
472          A->vect[j] = B->vect[j] & (0xffffffffUL >> m);
473 
474       }
475       j++;
476       for (; j < A->n; j++)
477          A->vect[j] = B->vect[j];
478    }
479 }
480 
481 
482 
483 
484 /* ********************************************************************** */
485 /* void ANDBVSelf(BitVect A, BitVect B)       				  */
486 /* Cette fonction effectue A &=B     					  */
487 /* ********************************************************************** */
ANDBVSelf(BitVect * A,BitVect * B)488 void ANDBVSelf (BitVect * A, BitVect * B)
489 {
490    int i;
491 
492    if ((A->n != B->n)) {
493       printf ("Error in ANDBVSelf(): Vectors of different sizes\n");
494       exit (1);
495    }
496 
497    for (i = 0; i < B->n; i++)
498       A->vect[i] &= B->vect[i];
499 }
500 
501 
502 /* ********************************************************************** */
503 /* void XORBVSelf(BitVect A, BitVect B)     				  */
504 /* Cette fonction effectue A ^=B    					  */
505 /* ********************************************************************** */
XORBVSelf(BitVect * A,BitVect * B)506 void XORBVSelf (BitVect * A, BitVect * B)
507 {
508    int i;
509 
510    if ((A->n != B->n)) {
511       printf ("Error in XORBVSelf(): Vectors of different sizes\n");
512       exit (1);
513    }
514 
515    for (i = 0; i < B->n; i++)
516       A->vect[i] ^= B->vect[i];
517 }
518 
519 
520 /* ********************************************************************** */
521 /* void BVLShift( BitVect R, BitVect A, int n )                      */
522 /* Effectue : R = A << n ;                                                */
523 /* ********************************************************************** */
BVLShift(BitVect * R,BitVect * A,int n)524 void BVLShift (BitVect * R, BitVect * A, int n)
525 {
526    int i;
527    int WLmn;
528    unsigned long temp;
529 
530    if ((R->n != A->n)) {
531       printf ("Error in BVLShift(): Vectors of different sizes\n");
532       exit (1);
533    }
534 
535    for (i = 0; i < A->n; i++)
536       R->vect[i] = A->vect[i];
537    while (n >= 32) {
538       for (i = 1; i < A->n; i++)
539          R->vect[i - 1] = R->vect[i];
540       R->vect[A->n - 1] = 0UL;
541       n -= 32;
542    }
543    if (n > 0) {
544       WLmn = WL - n;
545       R->vect[0] <<= n;
546       for (i = 1; i < A->n; i++) {
547          temp = R->vect[i] >> WLmn;
548          R->vect[i - 1] |= temp;
549          R->vect[i] <<= n;
550       }
551    }
552 
553 }
554 
555 
556 /* ********************************************************************** */
557 /* void BVRShift( BitVect R, BitVect A, int n )                      */
558 /* Effectue : R = A >> n ;                                                */
559 /* ********************************************************************** */
BVRShift(BitVect * R,BitVect * A,int n)560 void BVRShift (BitVect * R, BitVect * A, int n)
561 {
562    int i;
563    int WLmn;
564    unsigned long temp;
565    if ((R->n != A->n)) {
566       printf ("Error in BVRShift(): Vectors of different sizes\n");
567       exit (1);
568    }
569 
570    for (i = 0; i < A->n; i++)
571       R->vect[i] = A->vect[i];
572    while (n >= 32) {
573       for (i = A->n; i > 1; i--)
574          R->vect[i - 1] = R->vect[i - 2];
575       R->vect[0] = 0UL;
576       n -= 32;
577    }
578    if (n > 0) {
579       WLmn = WL - n;
580       R->vect[A->n - 1] >>= n;
581       for (i = A->n - 2; i >= 0; i--) {
582          temp = R->vect[i] << WLmn;
583          R->vect[i + 1] |= temp;
584          R->vect[i] >>= n;
585       }
586    }
587 }
588 
589 
590 /* ********************************************************************** */
591 /* void BVLShiftSelf( BitVect R, int n )                             */
592 /* Effectue : R <<= n ;                                                   */
593 /* ********************************************************************** */
BVLShiftSelf(BitVect * R,int n)594 void BVLShiftSelf (BitVect * R, int n)
595 {
596    int i;
597    int WLmn;
598    unsigned long temp;
599 
600    while (n >= 32) {
601       for (i = 1; i < R->n; i++)
602          R->vect[i - 1] = R->vect[i];
603       R->vect[R->n - 1] = 0UL;
604       n -= 32;
605    }
606    if (n > 0) {
607       WLmn = WL - n;
608       R->vect[0] <<= n;
609       for (i = 1; i < R->n; i++) {
610          temp = R->vect[i] >> WLmn;
611          R->vect[i - 1] |= temp;
612          R->vect[i] <<= n;
613       }
614    }
615 }
616 
617 
618 /* ********************************************************************** */
619 /* void BVLS1Self( BitVect R )                                            */
620 /* Effectue : R <<= 1 ;                                                   */
621 /* Version specialisee de BVLShiftSelf pour utilisation frequente.        */
622 /* ********************************************************************** */
BVLS1Self(BitVect * R)623 void BVLS1Self (BitVect * R)
624 {
625    int i;
626    R->vect[0] <<= 1;
627    for (i = 1; i < R->n; i++) {
628       if (R->vect[i] & MC)
629          R->vect[i - 1] |= 0x1UL;
630       R->vect[i] <<= 1;
631    }
632 }
633 
634 
635 /* ********************************************************************** */
636 /* void BVRShiftSelf( BitVect R, int n )                             */
637 /* Effectue : R >>= n ;                                                   */
638 /* ********************************************************************** */
BVRShiftSelf(BitVect * R,int n)639 void BVRShiftSelf (BitVect * R, int n)
640 {
641    int i;
642    int WLmn;
643    unsigned long temp;
644 
645    while (n >= 32) {
646       for (i = R->n - 1; i > 0; i--)
647          R->vect[i] = R->vect[i - 1];
648       R->vect[0] = 0UL;
649       n -= 32;
650    }
651    if (n > 0) {
652       WLmn = WL - n;
653       R->vect[R->n - 1] >>= n;
654       for (i = R->n - 2; i >= 0; i--) {
655          temp = R->vect[i] << WLmn;
656          R->vect[i + 1] |= temp;
657          R->vect[i] >>= n;
658       }
659    }
660 }
661 
662 
663 /* ********************************************************************** */
664 /* void invertBV(Bitvect A)                                               */
665 /* fait A ~=A                                                             */
666 /* ********************************************************************** */
InverseBV(BitVect * A)667 void InverseBV (BitVect * A)
668 {
669    int i;
670    for (i = 0; i < A->n; i++)
671       A->vect[i] = ~A->vect[i];
672 }
673 
674 
675 /* ********************************************************************** */
676 /* lebool CheckCD( BitVect ds1, BitVect ds2 )                            */
677 /* Verifie si les ensembles ds1 et ds2 ont des bits communs.              */
678 /* ********************************************************************** */
VerifBitsCommuns(BitVect * ds1,BitVect * ds2)679 lebool VerifBitsCommuns (BitVect * ds1, BitVect * ds2)
680 {
681    int i;
682    unsigned long temp = 0UL;
683    if ((ds1->n != ds2->n)) {
684       printf ("Error in VerifBitsCommuns(): Vectors of different sizes\n");
685       exit (1);
686    }
687    for (i = 0; i < ds1->n; i++)
688       temp |= (ds1->vect[i] & ds2->vect[i]);
689    if (temp)
690       return TRUE;
691    else
692       return FALSE;
693 }
694 
695 
696 /* -------------------------------------------- */
697 /* Fonctions pour la manipulation des matrices. */
698 /* -------------------------------------------- */
BVCanonic(BitVect * A,int l)699 void BVCanonic (BitVect * A, int l)
700 {
701    int n;
702    PutBVToZero (A);
703    n = l / WL;
704    if (n > A->n) {
705       printf
706          ("Error in  BVCanonic(): vector A is not long enough to store  BVCanonic[%d].\n",
707          l);
708       exit (1);
709    }
710    A->vect[n] = 0x80000000UL >> (l - n * WL);
711 }
712 
Mask(BitVect * A,int l)713 void Mask (BitVect * A, int l)
714 {
715 
716    InvMask (A, l);
717    InverseBV (A);
718 }
719 
InvMask(BitVect * A,int l)720 void InvMask (BitVect * A, int l)
721 {
722    AllOnes (A);
723    BVRShiftSelf (A, l);
724 }
725 
AllOnes(BitVect * A)726 void AllOnes (BitVect * A)
727 {
728    int i;
729    for (i = 0; i < A->n; i++)
730       A->vect[i] = 0xffffffffUL;
731 }
732 
AllocBV(BitVect * A,int l)733 void AllocBV (BitVect * A, int l)
734 {
735    int n;
736    n = (l - 1) / WL + 1;
737    A->vect = (unsigned long *) calloc ((size_t) n, sizeof (unsigned long));
738    A->n = n;
739 }
740 
FreeBV(BitVect * A)741 void FreeBV (BitVect * A)
742 {
743    if (A->vect != NULL) {
744       free (A->vect);
745    }
746    A->n = 0;
747 }
748 
AllocMat(Matrix * m,int nblines,int l,int t)749 void AllocMat (Matrix * m, int nblines, int l, int t)
750 {
751    int i, j;
752    m->lignes = (BitVect **) calloc ((size_t) nblines, sizeof (BitVect *));
753    for (i = 0; i < nblines; i++) {
754       if (!(m->lignes[i] = (BitVect *) calloc ((size_t) t, sizeof (BitVect)))) {
755          printf ("\n*** Memoire insuffisante pour AllocMat() ! nl=%d***\n",
756             nblines);
757          exit (1);
758       }
759       for (j = 0; j < t; j++)
760          AllocBV (&(m->lignes[i][j]), l);
761 
762    }
763    m->nblignes = nblines;
764    m->t = t;
765    m->l = l;
766 
767 }
768 
769 
770 /* ********************************************************************** */
771 /* void FreeSpace( Matrix m, int nl )                               */
772 /* Libere l'espace des nl vecteurs de la matrice m.                       */
773 /* ********************************************************************** */
FreeMat(Matrix * m)774 void FreeMat (Matrix * m)
775 {
776    int i, j;
777    for (i = 0; i < m->nblignes; i++) {
778       for (j = 0; j < m->t; j++)
779          FreeBV (&(m->lignes[i][j]));
780       free (m->lignes[i]);
781    }
782    free (m->lignes);
783    m->nblignes = 0;
784    m->l = 0;
785    m->t = 0;
786 
787 }
788 
789 
790 /* ********************************************************************** */
791 /* void CopyMat( Matrix m, Matrix ms, int nl, int t )                   */
792 /* Cette procedure sert a copier les t premiers BitVect des nl premieres  */
793 /* lignes de la matrice ms dans la matrice de travail m.                  */
794 /* ********************************************************************** */
CopyMat(Matrix * m,Matrix * ms,int nl,int t)795 void CopyMat (Matrix * m, Matrix * ms, int nl, int t)
796 {
797    int i, j;
798    if (m == NULL) {
799       AllocMat (m, ms->nblignes, ms->l, ms->t);
800    } else if ((ms->nblignes < nl) || (ms->t < t)) {
801       printf ("Error in CopyMat(): source matrix too small %d\n",
802          ms->nblignes / ms->t);
803       exit (1);
804    } else if ((m->nblignes < nl) || (m->t < t)) {
805       printf ("Error in CopyMat(): destination matrix too small\n");
806       exit (1);
807    }
808    for (i = 0; i < nl; i++)
809       for (j = 0; j < t; j++) {
810          CopyBV (&(m->lignes[i][j]), &(ms->lignes[i][j]));
811       }
812 }
813 
814 
815 /* ********************************************************************** */
816 /* void CopyNTupleMat( Matrix m, Matrix ms, int nl,     */
817 /*                     int *colonnes, int t )                   */
818 /* Cette procedure sert a copier les t-1 BitVect indiqu'es par le vecteur  */
819 /* *colonnes plus la colonne 0 des nl premieres lignes de la matrice ms   */
820 /* dans la matrice de travail m.       */
821 /* ********************************************************************** */
CopyNTupleMat(Matrix * m,Matrix * ms,int nl,int * colonnes,int t)822 void CopyNTupleMat (Matrix * m, Matrix * ms, int nl, int *colonnes, int t)
823 {
824 
825    int i, j, k, n;
826 
827    if (m == NULL)
828       AllocMat (m, ms->nblignes, ms->l, t);
829    else {
830       if ((ms->nblignes != m->nblignes) || (ms->l != m->l))
831          printf ("Error in CopieNTupleMat(): matrices of different sizes\n");
832    }
833    n = (ms->l - 1) / WL;
834    for (i = 0; i < nl; i++) {
835       for (k = 0; k <= n; k++)
836          (m->lignes[i])[0].vect[k] = (ms->lignes[i])[0].vect[k];
837       for (j = 1; j < t; j++)
838          for (k = 0; k <= n; k++)
839             (m->lignes[i])[j].vect[k] =
840                (ms->lignes[i])[colonnes[j - 1]].vect[k];
841 
842    }
843 }
844 
845 
846 /* ********************************************************************** */
847 /* void SwapVect( Matrix m, int i, int j )                     */
848 /* Pour interchanger les lignes i et j de la matrice de travail m.        */
849 /* ********************************************************************** */
ExchangeVect(Matrix * m,int i,int j)850 void ExchangeVect (Matrix * m, int i, int j)
851 {
852    BitVect *p;
853    if (i != j) {
854       p = m->lignes[i];
855       m->lignes[i] = m->lignes[j];
856       m->lignes[j] = p;
857    }
858 }
859 
860 
TransposeMatrices(Matrix * T,Matrix * M,int mmax,int smax,int L)861 void TransposeMatrices (Matrix * T, Matrix * M, int mmax, int smax, int L)
862 {
863 
864    int s, l, m;
865 
866    for (s = 0; s < smax; s++)
867       for (l = 0; l < L; l++) {
868          PutBVToZero (&T->lignes[l][s]);
869          for (m = 0; m < mmax; m++) {
870             /* printf("m=%d l=%d s=%d\n",m,l,s);fflush(stdout); */
871             if (M->lignes[m][s].vect[0] & (0x80000000UL >> l)) {
872                T->lignes[l][s].vect[0] |= (0x80000000UL >> m);
873             }
874          }
875       }
876 }
877 
878 
879 /* ********************************************************************** */
880 /* void XorVect( Matrix m,                                               */
881 /*               int r, int s, int min, int max )     */
882 /* Effectue un Xor entre la s-eme et la r-eme ligne de la matrice de      */
883 /* travail m pour les colonnes (BitVect) min a max-1 seulement.           */
884 /* Le resultat est mis dans la r-eme ligne. ( m[r] ^= m[s] )              */
885 /* ********************************************************************** */
XorVect(Matrix * m,int r,int s,int min,int max)886 void XorVect (Matrix * m, int r, int s, int min, int max)
887 {
888    int j;
889    for (j = min; j < max; j++)
890       XORBVSelf (&(m->lignes[r][j]), &(m->lignes[s][j]));
891 }
892 
893 
894 /* ********************************************************************** */
895 /* void displaymat(Matrix m, int t, int l, int kg)        */
896 /* Affiche la matrice m sur kg lignes par t x l colonnes  		  */
897 /* ********************************************************************** */
DispMat(Matrix * m,int t,int l,int kg,lebool mathematica)898 void DispMat (Matrix * m, int t, int l, int kg, lebool mathematica)
899 {
900    int i, j;
901 
902    i = kg;
903 
904    printf ("\n");
905    if (mathematica)
906       printf ("{");
907    for (i = 0; i < kg; i++) {
908       if (!mathematica)
909          printf ("[");
910       for (j = 0; j < t; j++) {
911          DispBitVect (&(m->lignes[i][j]), l, mathematica);
912       }
913       if (mathematica) {
914          if (i != kg - 1)
915             printf (",\n");
916          else
917             printf ("}\n");
918       } else
919          printf ("]\n");
920    }
921    printf ("\n\n");
922 
923 }
924 
925 
926 /* ********************************************************************** */
927 /* void displaybitvect(BitVect A, int l)    			  */
928 /* Affiche le BitVect A sur l bits seulement    			  */
929 /* ********************************************************************** */
DispBitVect(BitVect * A,int l,int mathematica)930 void DispBitVect (BitVect * A, int l, int mathematica)
931 {
932    int j;
933    unsigned Un;
934    Un = 1UL;
935    j = 0;
936    if (mathematica) {
937       printf ("{");
938       while (j < l - 1) {
939          printf ("%ld,",
940             (A->vect[j / WL] >> (((WL * A->n) - j - 1) % WL)) & Un);
941          j++;
942       }
943       printf ("%ld}", (A->vect[j / WL] >> (((WL * A->n) - j - 1) % WL)) & Un);
944    } else
945       while (j < l) {
946          printf ("%ld",
947             (A->vect[j / WL] >> (((WL * A->n) - j - 1) % WL)) & Un);
948          j++;
949       }
950 }
951 
952 
953 /* ********************************************************************** */
954 /* MultMatrixByMatrix ( Matrix *A, Matrix *B, Matrix *C)        	  */
955 /* Fait la multiplication matricielle A = B x C    			  */
956 /* ********************************************************************** */
957 
MultMatrixByMatrix(Matrix * A,Matrix * B,Matrix * C)958 void MultMatrixByMatrix (Matrix * A, Matrix * B, Matrix * C)
959 {
960    int i, j;
961    if (B->l != C->nblignes) {
962       printf ("Tailles de matrices non-compatibles, kaput.\n");
963       exit (1);
964    }
965 
966    if (A->nblignes != B->nblignes || A->l != C->l) {
967       printf ("Matrice preallouee de mauvaise taille.\n");
968       exit (1);
969    }
970 
971    for (i = 0; i < A->nblignes; i++)
972       PutBVToZero (A->lignes[i]);
973 
974    for (i = 0; i < B->nblignes; i++)
975       for (j = 0; j < B->l; j++) {
976          if (ValBitBV (B->lignes[i], j))
977             XORBVSelf (A->lignes[i], C->lignes[j]);
978       }
979 }
980 
981 
982 /* ********************************************************************** */
983 /* MatrixTwoPow ( Matrix *A, Matrix *B, unsigned int e)          	  */
984 /* Fait l'exponentiation : A = B^(2^e)          			  */
985 /* ********************************************************************** */
986 
MatrixTwoPow(Matrix * A,Matrix * B,unsigned int e)987 void MatrixTwoPow (Matrix * A, Matrix * B, unsigned int e)
988 {
989    unsigned int i;
990    Matrix tempMatrix;
991    Matrix *AA = &tempMatrix;
992 
993    if (B->nblignes != B->l) {
994       printf ("Matrice non carree.\n");
995       exit (1);
996    }
997    if (A->nblignes != B->nblignes || A->l != B->l) {
998       printf ("Matrice preallouee de mauvaise taille.\n");
999       exit (1);
1000    }
1001 
1002    AllocMat (AA, B->nblignes, B->l, 1);
1003 
1004    if (e == 0) {
1005       CopyMat (A, B, B->nblignes, 1);
1006       return;
1007    }
1008    /* A = B^(2^1) */
1009    MultMatrixByMatrix (A, B, B);
1010 
1011    for (i = 1; i < e - 1; i += 2) {
1012      /* AA = A * A */
1013       MultMatrixByMatrix (AA, A, A);
1014 
1015       /* A = AA * AA */
1016       MultMatrixByMatrix (A, AA, AA);
1017    }
1018 
1019    if (i == e - 1) {
1020      /* AA = A * A */
1021       MultMatrixByMatrix (AA, A, A);
1022       CopyMat (A, AA, AA->nblignes, 1);
1023    }
1024 
1025    FreeMat (AA);
1026 }
1027 
1028 /* ********************************************************************** */
1029 /* MatrixPow ( Matrix *A, Matrix *B, unsigned int e)            	  */
1030 /* Fait l'exponentiation : A = B^e               			  */
1031 /* ********************************************************************** */
1032 
1033 #ifdef USE_LONGLONG
MatrixPow(Matrix * A,Matrix * B,longlong e)1034 void MatrixPow (Matrix * A, Matrix * B, longlong e)
1035 #else
1036 void MatrixPow (Matrix * A, Matrix * B, long e)
1037 #endif
1038 {
1039    int i;
1040    Matrix C;
1041    Matrix D;
1042 
1043    if (B->nblignes != B->l) {
1044       printf ("Matrice non carree.\n");
1045       exit (1);
1046    }
1047 
1048    if (A->nblignes != B->nblignes || A->l != B->l) {
1049       printf ("Matrice preallouee de mauvaise taille.\n");
1050       exit (1);
1051    }
1052 
1053    AllocMat (&C, B->nblignes, B->l, 1);
1054 
1055    if (e < 0) {
1056       InverseMatrix (&C, B);
1057       MatrixPow (A, &C, -e);
1058       FreeMat (&C);
1059       return;
1060    }
1061    AllocMat (&D, B->nblignes, B->l, 1);
1062 
1063    /* A = I */
1064    for (i = 0; i < A->nblignes; i++)
1065       BVCanonic (A->lignes[i], i);
1066 
1067    /* C = B^1 */
1068    CopyMat (&C, B, B->nblignes, 1);
1069 
1070    while (e) {
1071       if (e & 1) {
1072          CopyMat (&D, A, B->nblignes, 1);
1073          MultMatrixByMatrix (A, &D, &C);
1074       }
1075 
1076       e >>= 1;
1077 
1078       if (e) {
1079          CopyMat (&D, &C, B->nblignes, 1);
1080          MultMatrixByMatrix (&C, &D, &D);
1081       }
1082    }
1083 
1084    FreeMat (&C);
1085    FreeMat (&D);
1086 }
1087 
1088 
1089 /* FIN vectorsF2.c */
1090