nonsysolve.c

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

 *     function.
 *   - Input and output matrices are row majored and represented by
 *     one-dimension array.
 *   - It is needed to preallocate the memory space of mp_N and mp_D. 
 *
 */

void
nonsingSolvLlhsMM (const enum SOLU_POS solupos, const long n, \
		   const long m, mpz_t *mp_A, mpz_t *mp_B, mpz_t *mp_N, \
		   mpz_t mp_D)
{
  long basislen, extbasislen, k=0, maxk=0, i, rt, kincr, ks=0, j, l, minv;
  mpz_t mp_m, mp_alpha, mp_temp, mp_basisprod, mp_extbasisprod, mp_beta, \
    mp_nb, mp_db, mp_maxnb, mp_maxdb;
  mpz_t *mp_r;
  FiniteField *liftbasis, *extbdcoeff, *cmbasis, **extbasis;
  Double *liftbasisInv, **ARNS, **AInv, ***C, ***C1; 

  { mpz_init(mp_nb); mpz_init(mp_db); mpz_init(mp_maxnb); mpz_init(mp_maxdb); }
  { mpz_init(mp_m); mpz_init(mp_beta); mpz_init(mp_alpha); mpz_init(mp_temp); }
  { mpz_init(mp_basisprod); mpz_init(mp_extbasisprod); }

  /* find lifting basis such that the product of the basis elements is 
     approximately n*mp_alpha */
  maxMagnMP(mp_A, n, n, n, mp_alpha);
  liftbasis = findLiftbasisLarge(n, mp_alpha, &basislen);

  /* initialize for lifting */
  AInv = XMALLOC(Double *, basislen);
  for (i = 0; i < basislen; i++) 
    AInv[i] = XMALLOC(Double, n*n);
  minv = liftInitLlhs(basislen, liftbasis, n, mp_A, mp_basisprod, \
		      mp_extbasisprod, &extbasislen, &cmbasis,  &extbdcoeff, \
		      &liftbasisInv, AInv, &extbasis, &ARNS);

  /* if A^(-1) mod liftbasis[minv] does not exist, adjust liftbasis */
  while (minv != -1)
    {
      adBasis(minv, basislen, liftbasis);
      minv = liftInitLlhs(basislen, liftbasis, n, mp_A, mp_basisprod, \
			  mp_extbasisprod,  &extbasislen, &cmbasis, \
			  &extbdcoeff, &liftbasisInv, AInv, &extbasis, &ARNS);
    }
  mp_r = XMALLOC(mpz_t, m*n);
  for (i = 0; i < m*n; i++) { mpz_init(mp_r[i]); }
  mpz_set_ui(mp_D, 1);
  if (solupos == RightSolu) 
    {
      maxMagnMP(mp_B, n, m, m, mp_beta); 
      for (i = 0; i < m*n; i++) { mpz_set(mp_r[i], mp_B[i]); }
    }
  else if (solupos == LeftSolu) 
    {
      maxMagnMP(mp_B, m, n, n, mp_beta); 
      for (i = 0; i < m; i++)
	for (j = 0; j < n; j++) 
	  mpz_set(mp_r[j*m+i], mp_B[i*n+j]);
    }
  /* set up k and bound of numerator and denominator */
  liftbd(mp_basisprod, n, mp_alpha, mp_beta, &maxk, mp_maxnb, mp_maxdb, \
	 &k, mp_nb, mp_db);
  kincr = k;
  ks = 0;
  C = NULL;
  do {
    /* lifting kincr more steps */
    C1 = lift(solupos, kincr, n, m, basislen, extbasislen, mp_basisprod, \
	      mp_extbasisprod, liftbasis, cmbasis, extbdcoeff, liftbasisInv, \
	      mp_r, extbasis, AInv, ARNS);

    /* update the lifting coefficients */
    C = XREALLOC(Double **, C, k);
    for (i = 0; i < kincr; i++) { C[i+ks] = C1[i]; }
    XFREE(C1);
    ks = k;

    /* rational reconstruction */
    rt = soluRecon(solupos, k, basislen, n, m, mp_basisprod, liftbasis, \
		   cmbasis, C, mp_nb, mp_db, mp_N, mp_D);

    /* break the loop  when maximum step reached */
    if (k == maxk) { break; }

    /* rational reconstruction succeeds, check the result by magnitude bound */
    if (rt == 1)
      {
	if (solupos == RightSolu) { maxMagnMP(mp_N, n, m, m, mp_nb); }
	else if (solupos == LeftSolu) { maxMagnMP(mp_N, m, n, n, mp_nb); }
	mpz_mul(mp_nb, mp_nb, mp_alpha);
	mpz_mul_ui(mp_nb, mp_nb, n);
	mpz_pow_ui(mp_m, mp_basisprod, k);
	mpz_sub_ui(mp_m, mp_m, 1);
	mpz_divexact_ui(mp_m, mp_m, 2);
	if (mpz_cmp(mp_nb, mp_m) < 0)
	  {
	    mpz_mul(mp_db, mp_D, mp_beta);
	    if (mpz_cmp(mp_db, mp_m) < 0) 
	      break;
	  }
      }
    /* 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_basisprod, 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);

  for (i = 0; i < k; i++) 
    {
      for (l = 0; l < basislen; 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 < basislen; i++) { XFREE(AInv[i]); } { XFREE(AInv); }
  for (i = 0; i < extbasislen; i++) { XFREE(ARNS[i]); } { XFREE(ARNS); }
  { XFREE(liftbasis); XFREE(extbdcoeff); XFREE(cmbasis); 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_basisprod); mpz_clear(mp_extbasisprod); mpz_clear(mp_temp); }

  return;
}



/*
 * Calling Sequence:
 *   nonsingSolvRNSMM(solupos, basislen, n, m, basis, ARNS, mp_B, mp_N, mp_D)
 *
 * Summary:
 *   Solve nonsingular system of linear equations, where the left hand side 
 *   input matrix is represented in a RNS.
 *
 * Description:
 *   Given a n x n nonsingular matrix A represented in a RNS, a n x m or m x n
 *   mpz_t matrix mp_B, this function will compute the solution of the system
 *   1. AX = mp_B 
 *   2. XA = mp_B. 
 *   The parameter solupos controls whether the system is in the type of 1
 *   or 2.
 *   
 *   Since the unique solution X is a rational matrix, the function will
 *   output the numerator matrix mp_N and denominator mp_D respectively, 
 *   such that A(mp_N) = mp_D*mp_B or (mp_N)A = mp_D*mp_B. 
 *
 * Input: 
 *    solupos: enumerate, flag to indicate the system to be solved
 *           - solupos = LeftSolu: solve XA = mp_B
 *           - solupos = RightSolu: solve AX = mp_B
 *   basislen: long, dimension of RNS basis
 *          n: long, dimension of A
 *          m: long, column or row dimension of mp_B depending on solupos
 *      basis: 1-dim FiniteField array length basislen, RNS basis
 *       ARNS: 2-dim Double array, dimension basislen x n*n, representation of
 *             n x n input matrix A in RNS, where ARNS[i] = A mod basis[i]
 *       mp_B: 1-dim mpz_t array length n*m, representing the right hand side
 *             matrix of the system
 *           - solupos = LeftSolu: mp_B a m x n matrix
 *           - solupos = RightSolu: mp_B a n x m matrix
 * 
 * Output:
 *   mp_N: 1-dim mpz_t array length n*m, representing the numerator matrix 
 *         of the solution
 *       - solupos = LeftSolu: mp_N a m x n matrix
 *       - solupos = RightSolu: mp_N a n x m matrix
 *   mp_D: mpz_t, denominator of the solution
 *
 * Precondition: 
 *   - A must be a nonsingular matrix.
 *   - Any element p in RNS basis must satisfy 2*(p-1)^2 <= 2^53-1.
 *
 * Note:
 *   - It is necessary to make sure the input parameters are correct,
 *     expecially the dimension, since there is no parameter checks in the 
 *     function.
 *   - Input and output matrices are row majored and represented by
 *     one-dimension array.
 *   - It is needed to preallocate the memory space of mp_N and mp_D. 
 *
 */

void
nonsingSolvRNSMM (const enum SOLU_POS solupos, const long n, const long m,\
		  const long basislen, const FiniteField *basis, \
		  Double **ARNS, mpz_t *mp_B, mpz_t *mp_N, mpz_t mp_D)
{
  long liftbasislen, extbasislen, k=0, maxk=0, i, rt, kincr, ks=0, j, l, minv;
  mpz_t mp_m, mp_alpha, mp_liftbasisprod, mp_extbasisprod, mp_beta, \
    mp_nb, mp_db, mp_maxnb, mp_maxdb;
  mpz_t *mp_r;
  FiniteField *liftbasis, *extbdcoeff, *cmliftbasis, **extbasis;
  Double *liftbasisInv, **AExtRNS, **AInv, ***C, ***C1; 
  double tt, tt1=0, tt2=0;

#if HAVE_TIME_H
  clock_t ti;
#endif


  { mpz_init(mp_nb); mpz_init(mp_db); mpz_init(mp_maxnb); mpz_init(mp_maxdb); }
  { mpz_init(mp_m); mpz_init(mp_beta); mpz_init(mp_alpha); }
  { mpz_init(mp_liftbasisprod); mpz_init(mp_extbasisprod); }

  /* find lifting basis such that the product of the basis elements is 
     approximately n*mp_alpha */
  basisProd(basislen, basis, mp_alpha);
  liftbasis = findLiftbasisLarge(n, mp_alpha, &liftbasislen);

  /* initialization step before lifting */
  AInv = XMALLOC(Double *, liftbasislen);
  for (i = 0; i < liftbasislen; i++) 
    AInv[i] = XMALLOC(Double, n*n);

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
  ti = clock();
#endif

  minv = liftInitRNS(liftbasislen, liftbasis, basislen, basis, n, ARNS, \
		     mp_liftbasisprod, mp_extbasisprod, &extbasislen, \
		     &cmliftbasis, &extbdcoeff, &liftbasisInv, AInv, \
		     &extbasis, &AExtRNS);

  /* if A^(-1) mod liftbasis[minv] does not exist, adjust liftbasis */
  while (minv != -1)
    {
      adBasis(minv, liftbasislen, liftbasis);
      minv = liftInitRNS(liftbasislen, liftbasis, basislen, basis, n, ARNS, \
			 mp_liftbasisprod, mp_extbasisprod, &extbasislen, \
			 &cmliftbasis, &extbdcoeff, &liftbasisInv, AInv, \
			 &extbasis, &AExtRNS);
    }

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
  tt = (double)(clock()-ti)/CLOCKS_PER_SEC;
  printf("              lifting initialization time: %f\n", tt);
  printf("              lifting basis length: %ld\n", liftbasislen);
  for (i = 0; i < liftbasislen; i++) 
    printf("          %ld  ", liftbasis[i]);
  printf("\n");
  printf("              extended lifting basis length: %ld\n", extbasislen);
  for (i = 0; i < extbasislen; i++)
    printf("           %ld  ", extbasis[0][i]);
  printf("\n");
  fflush(stdout);
#endif

  mp_r = XMALLOC(mpz_t, m*n);
  for (i = 0; i < m*n; i++) { mpz_init(mp_r[i]); }
  mpz_set_ui(mp_D, 1);
  if (solupos == RightSolu) 
    {
      maxMagnMP(mp_B, n, m, m, mp_beta); 
      for (i = 0; i < m*n; i++) { mpz_set(mp_r[i], mp_B[i]); }
    }
  else if (solupos == LeftSolu) 
    {
      maxMagnMP(mp_B, m, n, n, mp_beta); 
      for (i = 0; i < m; i++)
	for (j = 0; j < n; j++) 
	  mpz_set(mp_r[j*m+i], mp_B[i*n+j]);
    }
  /* set up k and bound of numerator and denominator */
  liftbd(mp_liftbasisprod, n, mp_alpha, mp_beta, &maxk, mp_maxnb, mp_maxdb, \
	 &k, mp_nb, mp_db);
  kincr = k;
  ks = 0;
  C = NULL;
  do {

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
    ti = clock();
#endif

    /* lifting kincr more steps */
    C1 = lift(solupos, kincr, n, m, liftbasislen, extbasislen, \
	      mp_liftbasisprod,  mp_extbasisprod, liftbasis, cmliftbasis, \
	      extbdcoeff, liftbasisInv, mp_r, extbasis, AInv, AExtRNS);

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
    tt1 += (double)(clock()-ti);
#endif 

    /* update the lifting coefficients */
    C = XREALLOC(Double **, C, k);
    for (i = 0; i < kincr; i++) { C[i+ks] = C1[i]; }
    XFREE(C1);
    ks = k;

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
    ti = clock();
#endif

    /* rational reconstruction */
    rt = soluRecon(solupos, k, liftbasislen, n, m, mp_liftbasisprod, \
		   liftbasis, cmliftbasis, C, mp_nb, mp_db, mp_N, mp_D);

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
    tt2 += (double)(clock()-ti);
#endif

    /* break the loop  when maximum lifting step reached */
    if (k == maxk) { break; }

    /* rational reconstruction succeed, check the result by magnitude bound */
    if (rt == 1)
      {
	if (solupos == RightSolu) { maxMagnMP(mp_N, n, m, m, mp_nb); }
	else if (solupos == LeftSolu) { maxMagnMP(mp_N, m, n, n, mp_nb); }
	mpz_mul(mp_nb, mp_nb, mp_alpha);
	mpz_mul_ui(mp_nb, mp_nb, n);
	mpz_pow_ui(mp_m, mp_liftbasisprod, k);
	mpz_sub_ui(mp_m, mp_m, 1);
	mpz_divexact_ui(mp_m, mp_m, 2);
	if (mpz_cmp(mp_nb, mp_m) < 0)
	  {
	    mpz_mul(mp_db, mp_D, mp_beta);
	    if (mpz_cmp(mp_db, mp_m) < 0) 
	      break;