nonsysolve.c

IML package provides efficient routines to solve nonsingular systems of linear equations, certifie

	  }
      }

    /* increase k and reset mp_nb and mp_db */
    kincr = (long)(0.1*k) > 20 ? (long)(0.1*k) : 20;
    if (k+kincr >= maxk) 
      {
	kincr = maxk-k;
	k = maxk;
	mpz_set(mp_nb, mp_maxnb);
	mpz_set(mp_db, mp_maxdb);
	continue;
      }
    k += kincr;
    mpz_pow_ui(mp_m, mp_liftbasisprod, k);
    mpz_sub_ui(mp_m, mp_m, 1);
    mpz_divexact_ui(mp_m, mp_m, 2);
    mpz_sqrt(mp_nb, mp_m);
    mpz_set(mp_db, mp_nb);
  } while (1);

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
  printf("              total lifting steps: %ld\n", k);
  tt1 = tt1/CLOCKS_PER_SEC;
  tt2 = tt2/CLOCKS_PER_SEC;
  printf("              lifting time: %f\n", tt1);
  printf("              rational reconstruction time: %f\n", tt2);
  fflush(stdout);
#endif

  for (i = 0; i < k; i++) 
    {
      for (l = 0; l < liftbasislen; l++) { XFREE(C[i][l]); }
      XFREE(C[i]);
    }
  XFREE(C);
  for (i = 0; i < m*n; i++) { mpz_clear(mp_r[i]); } { XFREE(mp_r); }
  for (i = 0; i < 2; i++) { XFREE(extbasis[i]); } { XFREE(extbasis); }
  for (i = 0; i < liftbasislen; i++) { XFREE(AInv[i]); } { XFREE(AInv); }
  for (i = 0; i < extbasislen; i++) { XFREE(AExtRNS[i]); } { XFREE(AExtRNS); }
  { XFREE(liftbasis); XFREE(extbdcoeff); XFREE(cmliftbasis); }
  XFREE(liftbasisInv);
  { mpz_clear(mp_nb); mpz_clear(mp_db); mpz_clear(mp_maxnb); mpz_clear(mp_m); }
  { mpz_clear(mp_maxdb); mpz_clear(mp_beta); mpz_clear(mp_alpha); }
  { mpz_clear(mp_liftbasisprod); mpz_clear(mp_extbasisprod); }
  return;
}



/*
 * Calling Sequence:
 *   liftbasis <-- findLiftbasisSmall(n, mp_alpha, basislen)
 *
 * Summary:
 *   Compute the p-adic lifting basis
 *
 * Description:
 *   The function computes p-adic lifting basis liftbasis, such that after 
 *   every single lifting step, log[2](p) solution bits will be gained, where
 *   p is the product of the elements in the lifting basis. 
 *
 *   Let the lifting basis be 'liftbasis' and dimension be m, then 
 *   - elements in liftbasis are all primes
 *   - liftbasis[0] is the largest prime among all the primes less than 
 *     RNSbound(n)
 *   - liftbasis[i+1] is the next prime smaller than liftbasis[i] (i = 0..m-2)
 *   - product of first m-2 elements in liftbasis <= mp_maxInter
 *   - product of first m-1 elements in liftbasis > mp_maxInter
 *
 * Input:
 *          n: long, dimension of A
 *   mp_alpha: mpz_t, maximum magnitude of A
 *   basislen: pointer to a long integer, the dimension of lifting basis
 * 
 * Return:
 *   liftbasis: 1-dime FiniteField length *basislen, storing the lifting basis
 */

FiniteField *
findLiftbasisSmall (const long n, const mpz_t mp_alpha, long *basislen)
{
  long temp, len=0, count;
  FiniteField *liftbasis;
  mpz_t mp_temp, mp_bd, mp_prod;

  /* compute the upper bound of lifting basis */
  temp = RNSbound(n);
  mpz_init_set_ui(mp_temp, temp);

  /* compute n*mp_alpha */
  mpz_init_set_ui(mp_bd, n);
  mpz_mul(mp_bd, mp_bd, mp_alpha);
  mpz_init_set_ui(mp_prod, 1);
  liftbasis = NULL;
  while (mpz_cmp(mp_bd, mp_prod) > 0)
    {
      ++len;
      liftbasis = XREALLOC(FiniteField, liftbasis, len);
      while (mpz_probab_prime_p(mp_temp, 10) == 0) 
	mpz_sub_ui(mp_temp, mp_temp, 1); 
      liftbasis[len-1] = mpz_get_ui(mp_temp);
      mpz_sub_ui(mp_temp, mp_temp, 1);
      mpz_mul_ui(mp_prod, mp_prod, liftbasis[len-1]);
    } 

  /* increase length to 2*len+4 */  
  count = len+4;
  while (--count > 0)
    {
      ++len; 
      liftbasis = XREALLOC(FiniteField, liftbasis, len);
      while (mpz_probab_prime_p(mp_temp, 10) == 0) 
	mpz_sub_ui(mp_temp, mp_temp, 1); 
      liftbasis[len-1] = mpz_get_ui(mp_temp);
      mpz_sub_ui(mp_temp, mp_temp, 1);
    }
  *basislen = len;
  { mpz_clear(mp_temp); mpz_clear(mp_bd); mpz_clear(mp_prod); }
  return liftbasis;
}



/*
 * Calling Sequence:
 *   liftbasis <-- findLiftbasisLarge(n, mp_alpha, basislen)
 *
 * Summary:
 *   Compute the p-adic lifting basis
 *
 * Description:
 *   The function computes p-adic lifting basis liftbasis, such that after 
 *   every single lifting step, log[2](p) solution bits will be gained, where
 *   p is the product of the elements in the lifting basis. 
 *
 *   Let the lifting basis be 'liftbasis' and dimension be m, then 
 *   - elements in liftbasis are all primes
 *   - liftbasis[0] is the largest prime among all the primes less than 
 *     RNSbound(n)
 *   - liftbasis[i+1] is the next prime smaller than liftbasis[i] (i = 0..m-2)
 *   - product of first m-2 elements in liftbasis <= mp_maxInter
 *   - product of first m-1 elements in liftbasis > mp_maxInter
 *
 * Input:
 *          n: long, dimension of A
 *   mp_alpha: mpz_t, maximum magnitude of A
 *   basislen: pointer to a long integer, the dimension of lifting basis
 * 
 * Return:
 *   liftbasis: 1-dime FiniteField length *basislen, storing the lifting basis
 */

FiniteField *
findLiftbasisLarge (const long n, const mpz_t mp_alpha, long *basislen)
{
  long temp, len=0, count;
  FiniteField *liftbasis;
  mpz_t mp_temp, mp_bd, mp_prod;

  /* compute the upper bound of lifting basis */
  temp = RNSbound(n);
  mpz_init_set_ui(mp_temp, temp);

  /* compute n*mp_alpha */
  mpz_init_set_ui(mp_bd, n);
  mpz_mul(mp_bd, mp_bd, mp_alpha);
  mpz_init_set_ui(mp_prod, 1);
  liftbasis = NULL;
  while (mpz_cmp(mp_bd, mp_prod) > 0)
    {
      ++len;
      liftbasis = XREALLOC(FiniteField, liftbasis, len);
      while (mpz_probab_prime_p(mp_temp, 10) == 0) 
	mpz_sub_ui(mp_temp, mp_temp, 1); 
      liftbasis[len-1] = mpz_get_ui(mp_temp);
      mpz_sub_ui(mp_temp, mp_temp, 1);
      mpz_mul_ui(mp_prod, mp_prod, liftbasis[len-1]);
    } 

  count = 3;
  while (--count > 0)
    {
      ++len; 
      liftbasis = XREALLOC(FiniteField, liftbasis, len);
      while (mpz_probab_prime_p(mp_temp, 10) == 0) 
	mpz_sub_ui(mp_temp, mp_temp, 1); 
      liftbasis[len-1] = mpz_get_ui(mp_temp);
      mpz_sub_ui(mp_temp, mp_temp, 1);
      } 


  *basislen = len;
  { mpz_clear(mp_temp); mpz_clear(mp_bd); mpz_clear(mp_prod); }
  return liftbasis;
}



/*
 * Calling Sequence:
 *   adBasis(idx, basislen, liftbasis)
 *
 * Summary:
 *   Adjust the lifting basis if some element in the lifting basis is bad.
 *
 * Description:
 *   If A^(-1) mod liftbasis[idx] does not exist, then use this function to
 *   delete liftbasis[idx] and add the previous prime of liftbasis[basislen-1] 
 *   into the lifting basis. Meanwhile, move the elements in lifting basis to
 *   maintain the decreasing order.
 * 
 * Input: 
 *         idx: long, index in liftbasis where A^(-1) mod liftbasis[idx] 
 *              does not exist
 *    basislen: long, dimension of basis
 *   liftbasis: 1-dim FiniteField array length basislen, lifting basis
 *
 * Note: 
 *   liftbasis is update inplace.
 *
 */

void
adBasis (const long idx, const long basislen, FiniteField *liftbasis)
{
  long i;
  mpz_t mp_temp;
  mpz_init(mp_temp);

  /* move the elements forward */
  for (i = idx+1; i < basislen; i++) { liftbasis[i-1] = liftbasis[i]; }
  mpz_set_ui(mp_temp, liftbasis[basislen-1]);
  mpz_sub_ui(mp_temp, mp_temp, 1);
  while (mpz_probab_prime_p(mp_temp, 10) == 0) 
    mpz_sub_ui(mp_temp, mp_temp, 1); 
  liftbasis[basislen-1] = mpz_get_ui(mp_temp);

  mpz_clear(mp_temp);
  return;
}


/*
 * Calling Sequence:
 *   liftbd(mp_basisprod, n, mp_alpha, mp_beta, maxk, mp_maxnb, mp_maxdb, k,
 *          mp_nb, mp_db)
 *
 * Summary:
 *   Compute the initial number of lifting steps, initial solution bounds, 
 *   maximum number of lifting steps and maximum solution bounds.
 *
 * Description:
 *   Let the system of linear equations be Ax = b. The function computes
 *   the worst bound(Hadamard's bound) of denominator and numerators and 
 *   corresponding number of lifting steps maxk, which satisfies
 *   mp_basisprod^maxk >= 2*N*D and mp_basisprod^(maxk-1) < 2*N*D, where 
 *   D = n^ceil(n/2)*mp_alpha^n, worst case bound for denominator, 
 *   N = n^ceil(n/2)*mp_alpha^(n-1)*mp_beta, worst case bound for denominators.
 *
 *   The function also computes the initial number of lifting step k = 20 and
 *   the corresponding bound of denominator and numerators, such that
 *   mp_basisprod^k >= 2*N1*D1, mp_basisprod^(k-1) < 2*N1*D1, and N1 = D1, 
 *   where N1 is the estimate bound for numerators and D1 is the estimate
 *   bound for denominator.
 *
 * Input:
 *   mp_basisprod: mpz_t, product of elements of lifting basis
 *              n: long, dimension of matrix A
 *       mp_alpha: mpz_t, maximum magnitude of A
 *        mp_beta: mpz_t, maximum magnitude of b
 *
 * Output:
 *       maxk: pointer to long int, maximum number of lifting steps
 *   mp_maxnb: mpz_t, the worst numerator bound N
 *   mp_maxdb: mpz_t, the worst denominator bound D
 *          k: pointer to long int, estimated number of lifting steps
 *      mp_nb: mpz_t, estimated numerator bound N1
 *      mp_db: mpz_t, estimated denominator bound D1
 *
 */

void
liftbd (const mpz_t mp_basisprod, const long n, const mpz_t mp_alpha, \
	const mpz_t mp_beta, long *maxk, mpz_t mp_maxnb, mpz_t mp_maxdb, \
	long *k, mpz_t mp_nb, mpz_t mp_db)
{
  long n1;
  mpz_t mp_t1, mp_t2, mp_prod;

  if ((n % 2) == 0) { n1 = n/2; }
  else { n1 = n/2+1; }

  /* compute the worst case bound */
  mpz_init(mp_t1);
  mpz_init(mp_t2);
  mpz_ui_pow_ui(mp_t1, n, n1);
  mpz_pow_ui(mp_db, mp_alpha, n);
  mpz_mul(mp_maxdb, mp_db, mp_t1);
  mpz_pow_ui(mp_nb, mp_alpha, n-1);
  mpz_mul(mp_nb, mp_nb, mp_beta);
  mpz_mul(mp_maxnb, mp_nb, mp_t1);
  mpz_init_set(mp_prod, mp_maxdb);
  mpz_mul(mp_prod, mp_prod, mp_maxnb);
  mpz_mul_ui(mp_prod, mp_prod, 2);
  mpz_add_ui(mp_prod, mp_prod, 1);

  /* compute maxk */
  *maxk = 1;
  mpz_set(mp_t1, mp_basisprod);
  while (mpz_cmp(mp_t1, mp_prod) < 0) 
    {
      mpz_mul(mp_t1, mp_t1, mp_basisprod); 
      ++(*maxk);
    }

  /* compute k and estimate bound */
  *k = 20;
  mpz_pow_ui(mp_prod, mp_basisprod, *k);
  mpz_sub_ui(mp_prod, mp_prod, 1);
  mpz_divexact_ui(mp_prod, mp_prod, 2);
  mpz_sqrt(mp_nb, mp_prod);
  mpz_set(mp_db, mp_nb);
  if (*k >= *maxk) { mpz_set(mp_nb, mp_maxnb); mpz_set(mp_db, mp_maxdb); }
  mpz_clear(mp_t1);
  mpz_clear(mp_t2);
  mpz_clear(mp_prod);
  return;
}