1 /* gauss.c ; solve small linear systems over GF(2)
2 *
3 * Algorithm: Gaussian elimination as described in ...
4
5 Copyright 2007, 2008, 2009, 2010 Pierrick Gaudry, Francois Morain, Emmanuel Thom\'e, Paul Zimmermann
6
7 This file is part of CADO-NFS.
8
9 CADO-NFS is free software; you can redistribute it and/or modify it
10 under the terms of the GNU Lesser General Public License as published
11 by the Free Software Foundation; either version 2.1 of the License, or
12 (at your option) any later version.
13
14 CADO-NFS is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more
18 details.
19
20 You should have received a copy of the GNU Lesser General Public
21 License along with CADO-NFS; see the file COPYING. If not, write to
22 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
23 Boston, MA 02110-1301, USA.
24 */
25
26
27 /*
28 Compile with one of these:
29
30 gcc -I../ -I../utils -I../build/x86_64 -L../build/x86_64/utils gauss.c -lgmp -lutils -lm
31
32 gcc -O4 -DNDEBUG -funroll-loops -I../ -I../utils -I../build/x86_64 -L../build/x86_64/utils -DMULTI_ROW=3 -DNB_TEST=10 gauss.c -lgmp -lutils -lm
33
34
35 Note: for small sizes (100-200), MULTI_ROW=2 seems to be fine,
36 for larger sizes (also for small sizes on alpha), try MULTI_ROW=3
37
38 With NO_MAIN option, a .o file is generated, exporting only the
39 kernel function. In that case, the code does not use GMP at all.
40
41 The exported function is the following:
42
43 int kernel(mp_limb_t* mat, mp_limb_t** ker, int nrows, int ncols,
44 int limbs_per_row, int limbs_per_col);
45
46 which returns the dimension of the kernel and fill-in the pointer ker.
47 If just the dimension of kernel is wanted, set ker=NULL.
48
49 */
50
51 /*===========================================================================*/
52 /* includes, defines, and declarations... */
53 /* Warning: there are some (static) global variables */
54 /*===========================================================================*/
55
56 #include "cado.h" // IWYU pragma: keep
57 #include <stdio.h>
58 #include <stdlib.h>
59 #include <string.h>
60 #include <gmp.h> /* only used for setting a random matrix */
61 #include "bblas_gauss.h"
62 #include "macros.h" // ASSERT
63
64 #ifndef MULTI_ROW
65 #define MULTI_ROW 3
66 #endif
67
68 /* If the flag below is set, the interface to kernel() changes slightly, as the
69 * mp_limb_t ** ker argument is expected to hold space for pointers BUT
70 * these pointers need not be initialized to anything. The number of pointers
71 * eventually returned, and pointing to data, is exactly the dimension of
72 * the kernel.
73 *
74 * For tiny matrices, this feature does not make sense, as a flat storage
75 * for the kernel vectors is enough.
76 */
77 #define SAVE_KERNEL_MEMORY 0
78
79 #ifndef __GMP_H__
80 typedef unsigned long int mp_limb_t;
81 #endif
82
83 /* some global variables storing the data */
84 static int NROWS = 0;
85 static int NCOLS = 0;
86 static int LIMBS_PER_ROW = 0;
87 static mp_limb_t* matrix = NULL;
88 static mp_limb_t** ptr_rows = NULL;
89
90 /* functions declaration */
91
92 static void check_soundness();
93
94 #if 0 /* ununsed currently */
95 static void addPartialRows(int row, int pivot, mp_limb_t mask,
96 int j_cur, mp_limb_t **ptr);
97 #endif
98 static inline void addRows(int row, int pivot, mp_limb_t mask,
99 /*int j_cur,*/ mp_limb_t **ptr);
100 static inline void add2Rows(int row, int row2, int pivot, mp_limb_t mask,
101 /*int j_cur,*/ mp_limb_t **ptr_current);
102 static inline void add3Rows(int row, int row2, int row3, int pivot,
103 mp_limb_t mask, /* int j_current, */
104 mp_limb_t **ptr_current);
105 static inline int getPivot(mp_limb_t **ptr, mp_limb_t mask);
106
107 /*===========================================================================*/
108 /* functions definition */
109 /*===========================================================================*/
110
check_soundness()111 static void check_soundness(){
112 #if 0 // ndef NDEBUG
113 mp_limb_t x = 2342326;
114 int i;
115 ASSERT (LIMBS_PER_ROW*ULONG_BITS >= NCOLS);
116 /* touch everywhere in the matrix to look for SEGV's */
117 for (i = 0; i < NROWS; ++i)
118 x += matrix[LIMBS_PER_ROW*i];
119 if (x == matrix[LIMBS_PER_ROW*NROWS-1])
120 ASSERT(0); /* very unlikely to happen, but forces no optimizations */
121 #endif
122 }
123
simple_addrows(mp_limb_t * p,int d,int s,int stride)124 void simple_addrows(mp_limb_t * p, int d, int s, int stride)
125 {
126 if (!p) return;
127 mp_limb_t * src = p + s * stride;
128 mp_limb_t * dst = p + d * stride;
129 int j;
130 for(j = 0 ; j < stride ; j++) {
131 dst[j] ^= src[j];
132 }
133 }
134
135 /*
136 * Kernel computation.
137 *
138 * Put the vectors of the kernel in the (preallocated) table of cols ker.
139 * Return the dimension of the kernel.
140 * If just the dimension is wanted, then pass ker = NULL, and in that case,
141 * limbs_per_cols is meaningless.
142 *
143 * Variables description:
144 * ptr_current : table containing for each row the address of
145 * the limb which we are currently dealing with.
146 * These are incremented after ULONG_BITS
147 * pivot elimination.
148 * col_current : index of the current column
149 * j_current : the index of the limb containing the current column
150 * mask1 : mask with only one 1 set at the current column
151 * mask2 : mask with 1's set at column higher than the current
152 * set_pivot : table. at index j, contains the index of the row that
153 * used for pivoting column j. Initialized with -1.
154 * set_used : table. at index i, contains 1 if the i-th row has
155 * been used for pivoting, at some point, 0 otherwise.
156 * rank : the rank of the matrix
157 * i, j : used in for(;;) loops
158 */
159
160 struct gaussian_elimination_data {
161 int rank;
162 int nrows;
163 int ncols;
164 mp_limb_t* mat;
165 mp_limb_t* lmat;
166 int limbs_per_row;
167 int limbs_per_col;
168 int *set_pivot;
169 int *set_used;
170 };
171
172 #define ADDBIT_SLOW(r, l, i, j) \
173 r[i * l + j / ULONG_BITS] ^= 1UL << (j % ULONG_BITS)
174
gaussian_elimination(struct gaussian_elimination_data * G)175 void gaussian_elimination(struct gaussian_elimination_data * G)
176 {
177 int i;
178
179 /* store the data in the global variables */
180 matrix = G->mat;
181 NROWS = G->nrows;
182 NCOLS = G->ncols;
183 LIMBS_PER_ROW = G->limbs_per_row;
184 int limbs_per_col = G->limbs_per_col;
185 ptr_rows = (mp_limb_t **)malloc((NROWS+1)*sizeof(mp_limb_t *));
186 for (i = 0; i < NROWS+1; ++i)
187 ptr_rows[i] = matrix + LIMBS_PER_ROW*i;
188
189
190 /* if G->lmat == NULL, then don't bug the user if he just gave 0 for the
191 * otherwise unused limbs_per_col value */
192 ASSERT (G->lmat == NULL || limbs_per_col*ULONG_BITS >= NROWS);
193
194 G->set_pivot = (int *)malloc(NCOLS*sizeof(int));
195 G->set_used = (int *)malloc(NROWS*sizeof(int));
196
197 for (i = 0; i < NCOLS; ++i)
198 G->set_pivot[i] = -1;
199 for (i = 0; i < NROWS; ++i)
200 G->set_used[i] = 0;
201
202 mp_limb_t **ptr_current;
203 int j_current, col_current;
204 mp_limb_t mask1, mask2;
205 #ifdef VERBOSE
206 double st0 = seconds();
207 double st;
208 #endif
209
210 // shut up.
211 #ifdef VERBOSE
212 fprintf (stderr, "Using MULTI_ROW=%d\n", MULTI_ROW);
213 #endif
214
215 check_soundness();
216
217 /* initialize ptr_current, set_used and set_pivot */
218 ptr_current = (mp_limb_t **)malloc((NROWS)*sizeof(mp_limb_t *));
219 for (i = 0; i < NROWS; ++i)
220 ptr_current[i] = ptr_rows[i];
221
222 /* Main Loop: for each column, find the pivot and eliminate */
223 j_current = 0;
224 mask1 = 1UL;
225 mask2 = (-(1UL)) ^ (1UL);
226 col_current = 0;
227 while (col_current < NCOLS) {
228 int pivot;
229
230 /* Get the pivot */
231 pivot = getPivot(ptr_current, mask1);
232 if (pivot != -1) {
233 #ifndef NDEBUG
234 for (i = 0; i < pivot; ++i)
235 ASSERT (!((*ptr_current[i]) & mask1 ));
236 ASSERT (((*ptr_current[pivot]) & mask1 ));
237 ASSERT (G->set_used[pivot] == 0);
238 #endif
239 G->set_pivot[col_current] = pivot;
240 G->set_used[pivot] = 1;
241
242 /* Eliminate with the pivot */
243 #if MULTI_ROW == 1
244 for (i = pivot + 1; i < NROWS; ++i) {
245 /* is there a 1 in position (i, col_current) ? */
246 if ( (*ptr_current[i]) & mask1 ) {
247 addRows(i, pivot, mask1, /*j_current,*/ ptr_current);
248 simple_addrows(G->lmat, i, pivot, limbs_per_col);
249 }
250 }
251 #elif MULTI_ROW == 2
252 /* try with two rows at a time */
253 {
254 int i1;
255 i = pivot + 1;
256 i1 = -1;
257 while (i < NROWS) {
258 if ( (*ptr_current[i]) & mask1 ) {
259 if (i1 >= 0) {
260 add2Rows(i1, i, pivot, mask1, /*j_current,*/ ptr_current);
261 simple_addrows(G->lmat, i, pivot, limbs_per_col);
262 simple_addrows(G->lmat, i1, pivot, limbs_per_col);
263 i1 = -1;
264 } else {
265 i1 = i;
266 }
267 } /* end if */
268 ++i;
269 } /* end while */
270 if (i1 >= 0) {
271 addRows(i1, pivot, mask1, /*j_current,*/ ptr_current);
272 simple_addrows(G->lmat, i1, pivot, limbs_per_col);
273 }
274 } /* end block */
275 #elif MULTI_ROW == 3
276 /* try with three rows at a time */
277 {
278 int i1, i2;
279 i = pivot + 1;
280 i1 = -1; i2 = -1;
281 while (i < NROWS) {
282 if ( (*ptr_current[i]) & mask1 ) {
283 if (i1 >= 0) {
284 if (i2 >= 0) {
285 add3Rows(i1, i2, i, pivot, mask1, /*j_current,*/ ptr_current);
286 simple_addrows(G->lmat, i, pivot, limbs_per_col);
287 simple_addrows(G->lmat, i1, pivot, limbs_per_col);
288 simple_addrows(G->lmat, i2, pivot, limbs_per_col);
289 i1 = -1; i2 = -1;
290 }
291 else {
292 i2 = i;
293 }
294 } else {
295 i1 = i;
296 }
297 } /* end if */
298 ++i;
299 } /* end while */
300 if (i1 >= 0) {
301 if (i2 >= 0) {
302 add2Rows(i1, i2, pivot, mask1, /*j_current,*/ ptr_current);
303 simple_addrows(G->lmat, i1, pivot, limbs_per_col);
304 simple_addrows(G->lmat, i2, pivot, limbs_per_col);
305 } else {
306 addRows(i1, pivot, mask1, /*j_current,*/ ptr_current);
307 simple_addrows(G->lmat, i1, pivot, limbs_per_col);
308 }
309 }
310 } /* end block */
311 #else
312 #error MULTI_ROW must be 1, 2, or 3.
313 #endif
314
315 /* Purge the pivot line */
316 addRows(pivot, pivot, mask1, /*j_current,*/ ptr_current);
317 #if 0
318 #ifdef VERBOSE
319 if (G->lmat) {
320 printf("L:=\n");
321 printMatrix(G->lmat, NROWS, NROWS, limbs_per_col);
322 printf(";\n");
323 }
324 #endif
325 #endif
326 }
327
328 /* increment */
329 ++col_current;
330 mask1 <<= 1;
331 mask2 <<= 1;
332 if (!mask1) { /* we have done ULONG_BITS operations */
333 mask1 = 1UL;
334 mask2 = (-(1UL)) ^ (1UL);
335 j_current++;
336 ASSERT (j_current <= LIMBS_PER_ROW);
337 for (i = 0; i < NROWS; ++i)
338 ptr_current[i]++;
339 }
340
341 /* some verbosity... */
342 if ((col_current % 128) == 0)
343 {
344 #ifdef VERBOSE
345 st = (seconds () - st0); /* time in seconds */
346 fprintf (stderr, "done %d pivots in %1.0fs (est. %1.0fs)\n",
347 col_current,
348 st, (double) NCOLS * st / (double) col_current);
349 #endif
350 }
351
352 } /* end while */
353 free(ptr_current);
354 }
355
356
kernel(mp_limb_t * mat,mp_limb_t ** ker,int nrows,int ncols,int limbs_per_row,int limbs_per_col)357 int kernel(mp_limb_t* mat, mp_limb_t** ker, int nrows, int ncols,
358 int limbs_per_row, int limbs_per_col)
359 {
360 int i,j;
361
362 struct gaussian_elimination_data G[1] = {{
363 .mat = mat,
364 .nrows = nrows, .ncols = ncols,
365 .limbs_per_row = limbs_per_row,
366 .limbs_per_col = limbs_per_col,
367 .rank = 0, }};
368 gaussian_elimination(G);
369
370 #ifdef VERBOSE
371 printf("Pivots:=[");
372 for (i = 0; i < NCOLS-1; ++i)
373 printf("%d,", G->set_pivot[i]);
374 printf("%d];\n", G->set_pivot[NCOLS-1]);
375
376 printf("usedrows:= [");
377 for (i = 0; i < NROWS-1; ++i)
378 printf("%d,", G->set_used[i]);
379 printf("%d];\n", G->set_used[NROWS-1]);
380 #endif
381
382
383 /* if ker == NULL, then don't bug the user if he just gave 0 for the
384 * otherwise unused limbs_per_col value */
385 ASSERT (ker == NULL || limbs_per_col*ULONG_BITS >= NROWS);
386 /* Explore the set of unused rows, get the G->rank */
387 G->rank = NROWS;
388 for (i = 0; i < NROWS; ++i)
389 if (!G->set_used[i]) {
390 G->rank--;
391 if (ker != NULL) {
392 /* each unused row gives one element of the kernel */
393 #if SAVE_KERNEL_MEMORY
394 ker[0] = (mp_limb_t *) malloc(limbs_per_col * sizeof(mp_limb_t));
395 #endif
396
397 /* ker[0] contains the space for the new element of the kernel */
398 mp_limb_t *ptr = ptr_rows[i];
399
400 for (j = 0; j < limbs_per_col; ++j)
401 ker[0][j] = 0;
402 ker[0][i / ULONG_BITS] = 1UL << (i % ULONG_BITS);
403
404 mp_limb_t mask1 = 1UL;
405 j = 0;
406 while (j < NCOLS) {
407 if ((*ptr) & mask1) { /* have to recover the pivoting */
408 int pivot;
409 pivot = G->set_pivot[j];
410 (ker[0])[pivot / ULONG_BITS] |=
411 (1UL << (pivot % ULONG_BITS));
412 } /* end if */
413 /* increment */
414 ++j;
415 mask1 <<= 1;
416 if (!mask1) { /* we have done ULONG_BITS operations */
417 mask1 = 1UL;
418 ptr++;
419 }
420 } /* end while */
421 ker++;
422 } /* end if ker != NULL */
423 } /* end if / for, loop over all unused */
424
425 free(ptr_rows);
426 free(G->set_pivot);
427 free(G->set_used);
428
429 return NROWS - G->rank;
430 } /* end function kernel */
431
432 /* puts in lmat a full rank matrix of size nrows*nrows, and return a rank
433 * value r, such that the r first (least significant bits) rows of the
434 * matrix lmat * mat form a basis of the row span of the matrix mat,
435 * while the last (nrows-r) (most significant bits) rows are zero.
436 *
437 * lmat is expected to be already allocated, and to have a stride equal
438 * to limbs_per_col.
439 *
440 * elim_table is an extra field, which should only be used when ncols is
441 * small. If not null, one expects there a malloc()'ed area (see size
442 * below), where each block of 16*limbs_per_col limbs contains
443 * precomputed row combinations. Row combination of index i, in the k-th
444 * such block, is such that the multiplication of this row by the
445 * original matrix mat yields a row whose bits of indices k to k+4 are
446 * exactly set to i. This can be used to do a projection.
447 */
spanned_basis(mp_limb_t * lmat,mp_limb_t * mat,int nrows,int ncols,int limbs_per_row,int limbs_per_col,mp_limb_t * elim_table)448 int spanned_basis(mp_limb_t * lmat, mp_limb_t * mat, int nrows, int ncols,
449 int limbs_per_row, int limbs_per_col, mp_limb_t * elim_table)
450 {
451 int i;
452
453 #if 0
454 #ifdef VERBOSE
455 printf("M:=\n");
456 printMatrix(mat, nrows, ncols, limbs_per_row);
457 printf(";\n");
458 #endif
459 #endif
460 mp_limb_t * lmat2 = malloc(nrows * limbs_per_col * sizeof(mp_limb_t));
461 memset(lmat2, 0, nrows * limbs_per_col * sizeof(mp_limb_t));
462 for (i = 0; i < nrows; ++i)
463 ADDBIT_SLOW(lmat2, limbs_per_col, i, i);
464
465 struct gaussian_elimination_data G[1] = {{
466 .mat = mat, .lmat = lmat2,
467 .nrows = nrows, .ncols = ncols,
468 .limbs_per_row = limbs_per_row,
469 .limbs_per_col = limbs_per_col,
470 .rank = 0, }};
471 gaussian_elimination(G);
472
473
474 G->rank = 0;
475 for (i = 0; i < NROWS; ++i) {
476 if (G->set_used[i]) {
477 G->rank++;
478 }
479 }
480
481 int r = 0;
482 int h = G->rank;
483
484 for (i = 0; i < NROWS; ++i) {
485 if (G->set_used[i]) {
486 memcpy(lmat + r * limbs_per_col, lmat2 + i * limbs_per_col, limbs_per_col * sizeof(mp_limb_t));
487 r++;
488 } else {
489 memcpy(lmat + h * limbs_per_col, lmat2 + i * limbs_per_col, limbs_per_col * sizeof(mp_limb_t));
490 h++;
491 }
492 }
493
494 if (elim_table) {
495 // clearing is done progressively below
496 // memset(elim_table, 0, iceildiv(NCOLS, 4) * 16 * limbs_per_col);
497 mp_limb_t * zcol = malloc(limbs_per_col * sizeof(mp_limb_t));
498 memset(zcol, 0, limbs_per_col * sizeof(mp_limb_t));
499 mp_limb_t * e = elim_table;
500 for (i = 0; i < NCOLS ; i+=4, e+=16 * limbs_per_col) {
501 memset(e, 0, 16 * limbs_per_col);
502 mp_limb_t * killer[4];
503 for(int di = 0;di<4 ; di++) {
504 if (i+di >= NCOLS || G->set_pivot[i+di]<0)
505 killer[di] = zcol;
506 else
507 killer[di] = lmat2 + G->set_pivot[i+di] * limbs_per_col;
508 }
509 /* Do a gray code walk to fill in the table. */
510 mp_limb_t * ex, * ey;
511 int d = 0;
512 memset(e, 0, limbs_per_col); ex=e;
513 #define DO_BIT(z) do { \
514 d^=1<<z; ey=e+d*limbs_per_col; \
515 for(int j = 0 ; j < limbs_per_col ; j++) { \
516 ey[j] = ex[j] ^ killer[z][j]; \
517 } \
518 ex=ey; \
519 } while (0)
520 /* zero */ DO_BIT(0); DO_BIT(1); DO_BIT(0);
521 DO_BIT(2); DO_BIT(0); DO_BIT(1); DO_BIT(0);
522 DO_BIT(3); DO_BIT(0); DO_BIT(1); DO_BIT(0);
523 DO_BIT(2); DO_BIT(0); DO_BIT(1); DO_BIT(0);
524 }
525 #undef DO_BIT
526 free(zcol);
527 }
528
529 free(lmat2);
530 free(ptr_rows);
531 free(G->set_pivot);
532 free(G->set_used);
533
534 return r;
535 }
536
537
538 /* add the pivot row to the given row, expect for the given column
539 * mask gives the bit we have to keep
540 * j_current is the number of the limb which contains the column
541 */
542
addRows(int row,int pivot,mp_limb_t mask,mp_limb_t ** ptr_current)543 static inline void addRows(int row, int pivot, mp_limb_t mask,
544 /*int j_current,*/ mp_limb_t **ptr_current) {
545 int i;
546 mp_limb_t *ptr1, *ptr2;
547
548 ptr1 = ptr_rows[row];
549 ptr2 = ptr_rows[pivot];
550
551 for (i = 0; i < LIMBS_PER_ROW; ++i)
552 ptr1[i] ^= ptr2[i];
553 *ptr_current[row] |= mask;
554 } /* end function addRows */
555
556
557 /* add the pivot row to the two given rows, expect for the given column
558 * mask gives the bit we have to keep
559 * j_current is the number of the limb which contains the column
560 */
561
add2Rows(int row,int row2,int pivot,mp_limb_t mask,mp_limb_t ** ptr_current)562 static inline void add2Rows(int row, int row2, int pivot, mp_limb_t mask,
563 /*int j_current,*/ mp_limb_t **ptr_current)
564 {
565 // int i;
566 mp_limb_t *ptr1, *ptr2, *ptr3, *ptr_lim;
567
568 ptr1 = ptr_rows[row];
569 ptr2 = ptr_rows[row2];
570 ptr3 = ptr_rows[pivot];
571 ptr_lim = ptr_rows[pivot + 1];
572 while (ptr3 != ptr_lim) {
573 *ptr1++ ^= *ptr3;
574 *ptr2++ ^= *ptr3++;
575 }
576
577 *ptr_current[row] |= mask;
578 *ptr_current[row2] |= mask;
579 } /* end function add2Rows */
580
581
add3Rows(int row,int row2,int row3,int pivot,mp_limb_t mask,mp_limb_t ** ptr_current)582 static inline void add3Rows(int row, int row2, int row3, int pivot,
583 mp_limb_t mask, /* int j_current, */
584 mp_limb_t **ptr_current) {
585 // int i;
586 mp_limb_t *ptr1, *ptr2, *ptr3, *ptr_piv, *ptr_lim;
587
588 ptr1 = ptr_rows[row];
589 ptr2 = ptr_rows[row2];
590 ptr3 = ptr_rows[row3];
591 ptr_piv = ptr_rows[pivot];
592 ptr_lim = ptr_rows[pivot + 1];
593 while (ptr_piv != ptr_lim) {
594 mp_limb_t aux = *ptr_piv++;
595 *ptr1++ ^= aux;
596 *ptr2++ ^= aux;
597 *ptr3++ ^= aux;
598 }
599
600 *ptr_current[row] |= mask;
601 *ptr_current[row2] |= mask;
602 *ptr_current[row3] |= mask;
603 } /* end function add3Rows */
604
605 #if 0 /* unused currently */
606 /* add the pivot row to the given row, starting with the given column+1
607 Note: the mask and ptr_current allow to do this quickly. */
608 /*
609 TODO: do several elimination simultaneously to reduce the cost of the loop
610 */
611
612 static void addPartialRows(int row, int pivot, mp_limb_t mask, int j_current,
613 mp_limb_t **ptr_current) {
614 // int i;
615 mp_limb_t *ptr1, *ptr2, *ptrlim;
616
617 *ptr_current[row] ^= (*ptr_current[pivot] & mask);
618
619 ptr1 = ptr_current[row] + 1;
620 ptr2 = ptr_current[pivot] + 1;
621 ptrlim = ptr_rows[row+1];
622 while (ptr1 < ptrlim)
623 *ptr1++ ^= *ptr2++;
624 } /* end function addPartialRows */
625 #endif
626
627 /* Return the index of the first row with a one in the desired column. */
628 /* If the column is empty, return -1 */
629
getPivot(mp_limb_t ** ptr,mp_limb_t mask)630 static inline int getPivot(mp_limb_t **ptr, mp_limb_t mask) {
631 int i;
632 for (i = 0; i < NROWS; ++i)
633 if (*ptr[i]&mask)
634 return i;
635
636 return -1;
637 }
638
639