Skip to content

Commit

Permalink
Merge pull request #114 from mohabsafey/bug-fixes
Browse files Browse the repository at this point in the history 
debug (adds missing test when lifting the GB)
  • Loading branch information
@mohabsafey
mohabsafey authored Dec 30, 2023
2 parents 3c8ec30 + d22f8dc commit 0dd6ec8
Show file tree
Hide file tree
Showing 4 changed files with 109 additions and 52 deletions.
6 changes: 4 additions & 2 deletions src/msolve/duplicate.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -353,14 +353,16 @@ static inline void duplicate_data_mthread_gbtrace(int nthreads,
int32_t **leadmons_current,
trace_t **btrace){

const long len = num_gb[0] * (st->nvars-st->nev);

const len_t len = num_gb[0] * (st->nvars);

for(int i = 0; i < nthreads; i++){
leadmons_current[i] = (int32_t *)calloc(len, sizeof(int32_t));
}
/* leadmons_ori[0] has already been allocated*/
for(int i = 1; i < nthreads; i++){
leadmons_ori[i] = (int32_t *)calloc(len, sizeof(int32_t));
for(long j = 0; j < len; j++){
for(len_t j = 0; j < len; j++){
leadmons_ori[i][j] = leadmons_ori[0][j];
}
}
Expand Down
149 changes: 102 additions & 47 deletions src/msolve/lifting-gb.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,13 @@

#include<flint/fmpz.h>

#define NBCHECK 2
#ifndef MIN
#define MIN(x, y) ((x) > (y) ? (y) : (x))
#endif
#ifndef MAX
#define MAX(a,b) (((a)>(b))?(a):(b))
#endif

typedef struct{
uint32_t len; /* length of the encoded polynomial */
Expand Down Expand Up @@ -78,8 +84,8 @@ typedef struct{
mpz_t *den; /* lifted denominator */
mpz_t gden; /* guessed denominator */
mpz_t tmp;
int32_t start; /* indicates smallest index of poly whose coef has been lifted but not checked*/
int32_t end; /* indicates largest index of poly whose coef has been lifted but not checked*/
int32_t start; /* smallest index of poly whose witness coef has been lifted but not checked*/
int32_t end; /* largest index of poly whose witness coef has been lifted but not checked*/
int *check1; /* tells whether lifted data are ok with one more prime */
int *check2; /* tells whether lifted data are ok with two more primes */
int32_t S;
Expand All @@ -98,19 +104,19 @@ static inline void data_lift_init(data_lift_t dl,
dl->S = 0;
int32_t i;

dl->steps = calloc(nsteps, sizeof(int32_t));
dl->steps = (int32_t *)calloc(nsteps, sizeof(int32_t));
for(i = 0; i < nsteps; i++){
dl->steps[i] = steps[i];
}
dl->cstep = 0;
dl->lend = npol;
dl->crt_mult = 0;
dl->crt = malloc(sizeof(mpz_t) * dl->npol);
dl->crt = (mpz_t *)malloc(sizeof(mpz_t) * dl->npol);
for(int32_t i = 0; i < dl->npol; i++){
mpz_init(dl->crt[i]);
}
dl->recon = 0;
dl->coef = calloc(npol, sizeof(mpz_t) );
dl->coef =(int32_t *) calloc(npol, sizeof(int32_t) );

dl->num = malloc(sizeof(mpz_t) * npol);
for(i = 0; i < npol; i++){
Expand All @@ -123,7 +129,7 @@ static inline void data_lift_init(data_lift_t dl,
mpz_init_set_ui(dl->gden, 1);
mpz_init(dl->tmp);
dl->start = 0;
dl->end = 0;
dl->end = npol;
dl->check1 = calloc(npol, sizeof(int));
dl->check2 = calloc(npol, sizeof(int));

Expand Down Expand Up @@ -161,26 +167,26 @@ static inline void gb_modpoly_init(gb_modpoly_t modgbs,
int32_t *lm, int32_t *basis){
modgbs->alloc = alloc;
modgbs->nprimes = 0;
modgbs->primes = calloc(alloc, sizeof(uint64_t));
modgbs->cf_64 = calloc(alloc, sizeof(uint64_t));
modgbs->primes = (uint64_t *)calloc(alloc, sizeof(uint64_t));
modgbs->cf_64 = (uint64_t *)calloc(alloc, sizeof(uint64_t));
modgbs->ld = ld;
modgbs->nv = nv;
modgbs->modpolys = malloc(sizeof(modpolys_struct) * ld);
modgbs->modpolys = (modpolys_t *)malloc(sizeof(modpolys_t) * ld);

modgbs->mb = basis;
modgbs->ldm = calloc(nv*ld, sizeof(int32_t));
modgbs->ldm = (int32_t *)calloc(nv*ld, sizeof(int32_t));
for(int32_t i = 0; i < ld; i++){
for(int j = 0; j < nv; j++){
modgbs->ldm[i*nv+j] = lm[i*nv+j];
}
}
for(uint32_t i = 0; i < ld; i++){
modgbs->modpolys[i]->len = lens[i];
modgbs->modpolys[i]->cf_32 = malloc(sizeof(uint32_t **)*lens[i]);
modgbs->modpolys[i]->cf_zz = malloc(sizeof(mpz_t)*lens[i]);
modgbs->modpolys[i]->cf_qq = malloc(sizeof(mpz_t)*2*lens[i]);
modgbs->modpolys[i]->cf_32 = (uint32_t **)malloc(sizeof(uint32_t *)*lens[i]);
modgbs->modpolys[i]->cf_zz = (mpz_t *)malloc(sizeof(mpz_t)*lens[i]);
modgbs->modpolys[i]->cf_qq = (mpz_t *)malloc(sizeof(mpz_t)*2*lens[i]);
for(uint32_t j = 0; j < lens[i]; j++){
modgbs->modpolys[i]->cf_32[j] = calloc(sizeof(uint32_t), alloc);
modgbs->modpolys[i]->cf_32[j] = calloc(alloc, sizeof(uint32_t));
mpz_init(modgbs->modpolys[i]->cf_zz[j]);
}
for(uint32_t j = 0; j < 2 * lens[i]; j++){
Expand Down Expand Up @@ -538,7 +544,7 @@ static inline int modpgbs_set(gb_modpoly_t modgbs,
}
modgbs->modpolys[i]->cf_32[bc][modgbs->nprimes] = c;
bc--;
}
}
}

modgbs->nprimes++;
Expand Down Expand Up @@ -834,12 +840,11 @@ static inline void incremental_dlift_crt_full(gb_modpoly_t modgbs, data_lift_t d

/* all primes are assumed to be good primes */
mpz_mul_ui(prod_p, mod_p, (uint32_t)newprime);
for(int32_t k = dl->lstart; k < modgbs->ld; k++){
uint64_t c = modgbs->modpolys[k]->cf_32[coef[k]][modgbs->nprimes - 1 ];

mpz_CRT_ui(dl->crt[k], dl->crt[k], mod_p,
c, newprime, prod_p, dl->tmp, 1);
for(int32_t k = 0; k < dl->end; k++){
uint64_t c = modgbs->modpolys[k]->cf_32[coef[k]][modgbs->nprimes - 1 ];

mpz_CRT_ui(dl->crt[k], dl->crt[k], mod_p,
c, newprime, prod_p, dl->tmp, 1);
}
mpz_set(mod_p, prod_p);
}
Expand All @@ -859,7 +864,7 @@ static inline void crt_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dlift,
modpolys_t *polys = modgbs->modpolys;

for(int32_t k = start; k < end; k++){
if(dlift->check1[k]){
if(dlift->check1[k] >= 1){
for(int32_t l = 0; l < polys[k]->len; l++){
for(uint32_t i = 0; i < modgbs->nprimes-1; i++){
modgbs->cf_64[i] = polys[k]->cf_32[l][i];
Expand Down Expand Up @@ -916,19 +921,19 @@ static inline int ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,

modpolys_t *polys = modgbs->modpolys;
set_recdata(dl, rd1, rd2, mod_p);

for(int32_t k = start; k < end; k++){
if(dl->check1[k]){
if(dl->check1[k] >= 1){
mpz_set_ui(dl->tmp, 1);
for(int32_t l = 0; l < polys[k]->len; l++){

if(ratreconwden(rnum, rden, polys[k]->cf_zz[l], mod_p, dl->den[k], rd1)){
mpz_set(polys[k]->cf_qq[2*l], rnum);
mpz_set(polys[k]->cf_qq[2*l + 1], rden);
mpz_lcm(dl->tmp, dl->tmp, rden);
}
else{
/* fprintf(stderr, "[%d/%d]", k, modgbs->ld - 1); */
mpz_set_ui(dl->gden, 1);
dl->check2[k] = 0;
mpz_clear(rnum);
mpz_clear(rden);
return k;
Expand All @@ -941,7 +946,7 @@ static inline int ratrecon_lift_modgbs(gb_modpoly_t modgbs, data_lift_t dl,
mpz_set_ui(polys[k]->cf_qq[2*l+1], 1);
}
mpz_mul(polys[k]->lm, polys[k]->lm, dl->tmp);
/* dl->S++; */
dl->check2[k] = 1;
}
else{
mpz_clear(rnum);
Expand All @@ -967,36 +972,72 @@ static inline int verif_coef(mpz_t num, mpz_t den, uint32_t prime, uint32_t coef
return (c==coef);
}

static inline int verif_lifted_basis(gb_modpoly_t modgbs, data_lift_t dl,
int thrds){
/* verification of the basis is performed at the very end of the computation */
if(dl->check1[dl->end-1] == 0){
return 1;
}
int b = 1;
mpz_t den;
mpz_init(den);
for(int32_t k = 0; k < modgbs->ld; k++){
if(dl->check1[k]>=1 && dl->check2[k] > 0 && dl->check2[k] < NBCHECK){
for(int i = 0; i < thrds; i++){
uint32_t prime = modgbs->primes[modgbs->nprimes - (thrds - i) ];
for(int32_t c = 0; c < modgbs->modpolys[k]->len; c++){
mpz_mul(den, modgbs->modpolys[k]->lm, modgbs->modpolys[k]->cf_qq[2*c+1]);
uint32_t coef = modgbs->modpolys[k]->cf_32[c][modgbs->nprimes - (thrds - i) ];
int b = verif_coef(modgbs->modpolys[k]->cf_qq[2*c], den, prime, coef);
if(!b){
dl->check2[k] = 0;
mpz_set_ui(dl->gden, 1);
b = 0;
mpz_clear(den);
return b;
}
}
dl->check2[k]++;
}
}
}
mpz_clear(den);
return b;
}

static inline int verif_lifted_rational_wcoef(gb_modpoly_t modgbs, data_lift_t dl,
int thrds){
for(int32_t k = dl->lstart; k < dl->lend; k++){
if(!dl->check1[k]){
if(dl->check1[k]==0){
/* too early to perform the verification */
return k;
}
for(int i = 0; i < thrds; i++){

uint32_t prime = modgbs->primes[modgbs->nprimes - (thrds - i) ];
uint32_t coef = modgbs->modpolys[k]->cf_32[dl->coef[k]][modgbs->nprimes - (thrds - i) ];
int b = verif_coef(dl->num[k], dl->den[k], prime, coef);
int boo = verif_coef(dl->num[k], dl->den[k], prime, coef);

if(!b){
dl->check1[k] = 0;
if(!boo){
for(int32_t kk = k; kk < dl->end; kk++){
dl->check1[kk] = 0;
}
dl->start = k;
return k;
}
dl->start = MIN(dl->start + 1, dl->lend);
dl->check1[k]++;
}

dl->start++;
}
return -1;
}

/* returns 0 iff rational number could not be lifted */
static inline int reconstructcoeff(data_lift_t dl, int32_t i, mpz_t mod_p,
rrec_data_t recdata1, rrec_data_t recdata2){

int b = ratreconwden(dl->num[i], dl->den[i],
dl->crt[i], mod_p, dl->gden, recdata1);

if(b){
mpz_mul(dl->den[i], dl->den[i], dl->gden);
mpz_gcd(dl->tmp, dl->den[i], dl->num[i]);
Expand All @@ -1022,9 +1063,19 @@ static void ratrecon_gb(gb_modpoly_t modgbs, data_lift_t dl,
rrec_data_t recdata1, rrec_data_t recdata2,
int thrds, double *st_crt, double *st_rrec){

double st = realtime();
verif_lifted_rational_wcoef(modgbs, dl, thrds);
int b = verif_lifted_basis(modgbs, dl, thrds);
if(b == 0){
for(int32_t i = 0; i < modgbs->ld; i++){
if(dl->check2[i] == 0){
dl->lstart = i;
break;
}
}
}
*st_rrec += realtime()-st;
if(dl->lstart != dl->start){
double st = realtime();
crt_lift_modgbs(modgbs, dl, dl->lstart, dl->start);
*st_crt += realtime() - st;
st = realtime();
Expand All @@ -1033,8 +1084,7 @@ static void ratrecon_gb(gb_modpoly_t modgbs, data_lift_t dl,
}
dl->lstart = dl->start;

double st = realtime();

st = realtime();
incremental_dlift_crt_full(modgbs, dl,
dl->coef, mod_p, prod_p,
thrds);
Expand All @@ -1056,10 +1106,11 @@ static void ratrecon_gb(gb_modpoly_t modgbs, data_lift_t dl,
int b = reconstructcoeff(dl, i, mod_p,
recdata1, recdata2);
if(!b){
dl->check1[i] = 0;
break;
}
else{
dl->check1[i] = 1;
dl->check1[i]++;
}
}
}
Expand Down Expand Up @@ -1197,6 +1248,8 @@ int msolve_gbtrace_qq(
prime = next_prime(rand() % (1303905301 - (1<<30) + 1) + (1<<30));
}

/* prime = 1246973177; */

primeinit = prime;
msd->lp->p[0] = primeinit;
if(gens->field_char){
Expand Down Expand Up @@ -1333,7 +1386,6 @@ int msolve_gbtrace_qq(
}
msd->lp->p[0] = prime;
}

int nthrds = 1; /* mono-threaded mult-mid comp */
for(len_t i = 1; i < nthrds/* st->nthrds */; i++){
prime = next_prime(prime);
Expand Down Expand Up @@ -1368,7 +1420,6 @@ int msolve_gbtrace_qq(

/* nprimes += st->nthrds; */
nprimes += 1; /* at the moment, multi-mod comp is mono-threaded */

if(nprimes == 1){
if(info_level>2){
fprintf(stderr, "------------------------------------------\n");
Expand Down Expand Up @@ -1408,12 +1459,13 @@ int msolve_gbtrace_qq(

int lstart = dlift->lstart;
double ost_rrec = st_rrec;
double ost_crt = st_crt;

if(!bad){
ratrecon_gb(modgbs, dlift, msd->mod_p, msd->prod_p, recdata1, recdata2,
nthrds/* st->nthrds */, &st_crt, &st_rrec);
}
if(/* (st_crt -ost_crt) + */ (st_rrec - ost_rrec) > dlift->rr * stf4){
if((st_crt -ost_crt) + (st_rrec - ost_rrec) > dlift->rr * stf4){
dlift->rr = 2*dlift->rr;
if(info_level){
fprintf(stderr, "(->%d)", dlift->rr);
Expand All @@ -1424,23 +1476,26 @@ int msolve_gbtrace_qq(
fprintf(stderr, "{%d}", nprimes);
}
}
if(dlift->lstart != lstart && dlift->lstart < modgbs->ld - 1){
if(info_level){
fprintf(stderr, "<%.2f%%>", 100* (float)(dlift->lstart + 1)/modgbs->ld);
apply = 0;
for(len_t i = 0; i < modgbs->ld; i++){
if(dlift->check2[i] < NBCHECK){
apply = 1;
break;
}
lstart = dlift->lstart;
}
if(dlift->lstart >= modgbs->ld){
if(dlift->lstart != lstart){
if(info_level){
fprintf(stderr, "<100%%>\n");
fprintf(stderr, "CRT time = %.2f, Rational reconstruction time = %.2f\n", st_crt, st_rrec);
fprintf(stderr, "<%.2f%%>", 100* (float)MIN((dlift->lstart + 1), modgbs->ld)/modgbs->ld);
}
apply = 0;
lstart = dlift->lstart;
}
/* this is where learn could be reset to 1 */
/* but then duplicated datas and others should be free-ed */
}
}
if(info_level){
fprintf(stderr, "\nCRT time = %.2f, Rational reconstruction time = %.2f\n", st_crt, st_rrec);
}
if(info_level){
long nbits = max_bit_size_gb(modgbs);
fprintf(stderr, "Maximum bit size of the coefficients: %ld\n", nbits);
Expand Down
4 changes: 2 additions & 2 deletions src/msolve/msolve.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1575,13 +1575,13 @@ static inline int new_rational_reconstruction(mpz_param_t mpz_param,
on verifie que mpz_pol / lc(mpz_pol) mod prime = nm_pol
renvoie 1 si il faut faire le modular check.
**/
static inline int check_unit_mpz_nmod_poly(const long len,
static inline int check_unit_mpz_nmod_poly(const deg_t len,
const mpz_upoly_t mpz_pol,
const nmod_poly_t nm_pol,
const int32_t prime){
uint32_t lc = mpz_fdiv_ui(mpz_pol->coeffs[len - 1], prime);
lc = mod_p_inverse_32(lc, prime);
for(long i = 0; i < len; i++){
for(deg_t i = 0; i < len; i++){
uint64_t c = mpz_fdiv_ui(mpz_pol->coeffs[i], prime);
c *= (uint64_t)lc;
c = c % prime;
Expand Down
Loading

0 comments on commit 0dd6ec8

Please sign in to comment.