Introduction of the lifting of multiplication matrices (over the rationals)

src/msolve/duplicate.c

@@ -1,3 +1,22 @@
+/* This file is part of msolve.
+ *
+ * msolve is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * msolve is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with msolve.  If not, see <https://www.gnu.org/licenses/>
+ *
+ * Authors:
+ * Jérémy Berthomieu
+ * Christian Eder
+ * Mohab Safey El Din */
 
 static inline void duplicate_linear_data(int nthreads, int nvars, int nlins,
                                          nvars_t **blinvars, uint32_t **blineqs,

src/msolve/hilbert.c

@@ -2836,6 +2836,51 @@ static inline sp_matfglm_t * build_matrixn_from_bs(int32_t *lmb, long dquot,
   return matrix;
 }
 
+
+static inline void compute_modular_matrix(sp_matfglm_t *matrix,
+                                          mpq_matfglm_t mpz_mat,
+                                          uint32_t prime){
+  uint32_t len_xn = matrix->nrows; 
+  uint32_t dquot = matrix->ncols;
+  matrix->charac = prime;
+  uint64_t len1 = dquot * matrix->nrows;
+  int32_t len2 = dquot - matrix->nrows;
+
+
+  for(uint32_t i = 0; i < mpz_mat->nrows; i++){
+    uint64_t lc = mpz_fdiv_ui(mpz_mat->denoms[i], prime);
+    lc = mod_p_inverse_32(lc, prime);
+    uint32_t nc = i*mpz_mat->ncols;
+    for(uint32_t j = 0 ; j < mpz_mat->ncols; j++){
+      uint32_t mod = mpz_fdiv_ui(mpz_mat->dense_mat[2*(nc+j)], prime);
+      matrix->dense_mat[nc+j] =( ((uint64_t)mod) * lc )% prime;
+    }
+  }
+  for(int32_t i = 0; i < (dquot-len_xn); i++){
+    matrix->triv_idx[i] = mpz_mat->triv_idx[i];
+  }
+  for(int32_t i = 0; i < len2; i++){
+    matrix->triv_pos[i] = mpz_mat->triv_pos[i];
+  }
+  for(uint32_t i = 0; i < len_xn; i++){
+    matrix->dense_idx[i] = mpz_mat->dense_idx[i];
+  }
+  for(uint32_t i = 0; i < len_xn; i++){
+    matrix->dst[i] = mpz_mat->dst[i];
+  }
+#ifdef DEBUGLIFTMAT
+  fprintf(stderr, "\nModular matrix (prime = %u)\n", prime);
+  for(int i = 0; i < matrix->nrows; i++){
+    int nc = i*matrix->ncols;
+    for(int j = 0; j < matrix->ncols; j++){
+      fprintf(stderr, "%u, ", matrix->dense_mat[nc+j]);
+    }
+    fprintf(stderr, "\n");
+  }
+#endif
+}
+
+
 /**
 
    lmb is the monomial basis (of the quotient ring) given by ascending order.
@@ -2863,32 +2908,29 @@ static inline void build_matrixn_from_bs_trace_application(sp_matfglm_t *matrix,
 							   const int nv,
                                                            const long fc){
 
-  long len_xn = matrix->nrows; //get_div_xn(bexp_lm, bs->lml, nv, div_xn);
-
+  uint32_t len_xn = matrix->nrows; 
   matrix->charac = fc;
-  /* matrix->ncols = dquot; */
-  /* matrix->nrows = len_xn; */
-  long len1 = dquot * matrix->nrows;
-  long len2 = dquot - matrix->nrows;
+  uint64_t len1 = dquot * matrix->nrows;
+  int32_t len2 = dquot - matrix->nrows;
 
-  for(long i = 0; i < len1; i++){
+  for(uint64_t i = 0; i < len1; i++){
     matrix->dense_mat[i] = 0;
   }
-  for(long i = 0; i < len2; i++){
+  for(int32_t i = 0; i < (dquot-len_xn); i++){
     matrix->triv_idx[i] = 0;
   }
-  for(long i = 0; i < len2; i++){
+  for(int32_t i = 0; i < len2; i++){
     matrix->triv_pos[i] = 0;
   }
-  for(long i = 0; i < len_xn; i++){
+  for(uint32_t i = 0; i < len_xn; i++){
     matrix->dense_idx[i] = 0;
   }
-  for(long i = 0; i < len_xn; i++){
+  for(uint32_t i = 0; i < len_xn; i++){
     matrix->dst[i] = 0;
   }
 
   long pos = 0, k = 0;
-  for(long i = 0; i < bs->lml; i++){
+  for(uint32_t i = 0; i < bs->lml; i++){
     long len = bs->hm[bs->lmps[i]][LENGTH];
     if(i==div_xn[k]){
       len_gb_xn[k]=len;
@@ -2905,7 +2947,7 @@ static inline void build_matrixn_from_bs_trace_application(sp_matfglm_t *matrix,
   long l_dens = 0;
   long nrows = 0;
   long count = 0;
-  for(long i = 0; i < dquot; i++){
+  for(uint32_t i = 0; i < dquot; i++){
 
     long pos = -1;
     int32_t *exp = lmb + (i * nv);
@@ -2972,8 +3014,8 @@ static inline void build_matrixn_from_bs_trace_application(sp_matfglm_t *matrix,
     }
   }
   //Ici on support que les entres de matrix->dst sont initialisees a 0
-  for(long i = 0; i < matrix->nrows; i++){
-    for(long j = matrix->ncols - 1; j >= 0; j--){
+  for(uint32_t i = 0; i < matrix->nrows; i++){
+    for(int32_t j = matrix->ncols - 1; j >= 0; j--){
       if(matrix->dense_mat[i*matrix->ncols + j] == 0){
         matrix->dst[i]++;
       }

src/msolve/iofiles.c

@@ -670,7 +670,7 @@ static int get_coefficient_ff_and_term_from_line(char *line, int32_t nterms,
   char *prev_pos = NULL;
   size_t term_size = 50000;
   char *term  = (char *)malloc(term_size * sizeof(char));
-  long int cf_tmp  = 0; /** temp for coefficient value, possibly coeff is negative. */
+  int64_t cf_tmp  = 0; /** temp for coefficient value, possibly coeff is negative. */
 
   prev_pos = line;
   get_term(line, &prev_pos, &term, &term_size);
@@ -690,14 +690,14 @@ static int get_coefficient_ff_and_term_from_line(char *line, int32_t nterms,
       }
     }
     while (iv_tmp < 0) {
-      iv_tmp  +=  (int32_t)field_char; //MS change int -> long int
+      iv_tmp  +=  field_char; //MS change int -> long int
     }
     gens->cfs[pos]  = (int32_t)iv_tmp;
     store_exponent(term, gens, pos*gens->nvars);
     for(int j = 1; j < nterms; j++){
       get_term(line, &prev_pos, &term, &term_size);
       if (term != NULL) {
-        cf_tmp  = (int)strtol(term, NULL, 10);
+        cf_tmp  = (int64_t)strtol(term, NULL, 10);
 
         if (cf_tmp == 0) {
           if (term[0] == '-') {
@@ -707,9 +707,9 @@ static int get_coefficient_ff_and_term_from_line(char *line, int32_t nterms,
           }
         }
         while (cf_tmp < 0) {
-          cf_tmp  += (int)field_char;
+          cf_tmp  += field_char;
         }
-        gens->cfs[pos+j] = cf_tmp;
+        gens->cfs[pos+j] = (int32_t)(cf_tmp % field_char);
         store_exponent(term, gens, (pos+j)*gens->nvars);
       }
       //      store_exponent(term, basis, ht);
@@ -718,7 +718,6 @@ static int get_coefficient_ff_and_term_from_line(char *line, int32_t nterms,
     return 0;
   }
   free(term);
-
   return 1;
 }
 
@@ -930,7 +929,7 @@ static inline void get_data_from_file(char *fn, int32_t *nr_vars,
   *nr_gens = get_ngenerators(fn);
 
   const int max_line_size  = 1073741824;
-  char *line  = (char *)calloc((nelts_t)max_line_size, sizeof(char));
+  char *line  = (char *)malloc((nelts_t)max_line_size * sizeof(char));
 
   FILE *fh  = fopen(fn,"r");
 
@@ -967,8 +966,8 @@ static inline void get_data_from_file(char *fn, int32_t *nr_vars,
 }
 
 static inline void display_gens_ff(FILE *fh, data_gens_ff_t *gens){
-  long pos = 0;
-  int c;
+  int64_t pos = 0;
+  int32_t c;
   for(long i = 0; i < gens->ngens; i++){
     for(long j = 0; j < gens->lens[i]-1; j++){
       c = gens->cfs[pos+j];
@@ -983,7 +982,7 @@ static inline void display_gens_ff(FILE *fh, data_gens_ff_t *gens){
     }
     c = gens->cfs[pos+gens->lens[i]-1];
     if(c!= 1){
-      fprintf(fh, "%d", c);
+      fprintf(fh, "%d*", c);
       display_monomial(fh, gens, pos+gens->lens[i]-1, &gens->exps);
     }
     else{

src/msolve/libmsolve.c

@@ -33,3 +33,4 @@
 #include "../usolve/libusolve.h"
 #include "../neogb/libneogb.h"
 #include "msolve.c"
+