certsolve.c

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

		 compute the bound after compression in specialminSolveRNS */
	      mpz_init(mp_bd);
	      compressBoundMP(r, m, Pt, mp_A, mp_bd);
	      ver = specialminSolveRNS(certflag, 1, nullcol, r, m,\
				       basislen, mp_bd, basis, CRNS, mp_b1, \
				       mp_N, mp_D, mp_NZ, mp_DZ);
	      mpz_clear(mp_bd);
	      if (ver == 0) 
		{
		  for (i = 0; i < n; i++) { mpz_clear(mp_b1[i]); } 
		  for (i = 0; i < basislen; i++) { XFREE(CRNS[i]); } 
		  { XFREE(CRNS); XFREE(mp_b1); }
		  { XFREE(P); XFREE(rp); XFREE(Pt); XFREE(rpt); }

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
		  printf("certificate checking fails, repeat!\n");
#endif

		  continue; 
		}
	    }
	  for (i = 0; i < basislen; i++) { XFREE(CRNS[i]); } { XFREE(CRNS); }

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
	  tt = (double)(clock()-ti)/CLOCKS_PER_SEC;
	  printf("Minimial Denominator Solution Solving Time: %f\n", tt);
	  ti = clock();
#endif

	  /* not in full row rank, need to check whether <A_31|A_32>x = b_3 */
	  if (r < n)
	    {
	      C1RNS = XMALLOC(Double *, basislen);
	      for (i = 0; i < basislen; i++)
		{
		  C1RNS[i] = XMALLOC(Double, (n-r)*m);
		  for (j = 0; j < n-r; j++)
		    for (l = 0; l < m; l++)
		      C1RNS[i][j*m+l] = \
			(Double)mpz_fdiv_ui(mp_A[Pt[r+j]*m+rpt[l]], basis[i]);
		}

	      cmp = LVecSMatMulCmp(RightMul, basislen, n-r, m, basis, C1RNS, \
				   mp_D, mp_N, mp_b1+r);
	      for (i = 0; i < basislen; i++) { XFREE(C1RNS[i]); } 
	      XFREE(C1RNS);

	      /* if compare fails, restart */
	      if (cmp == 0) 
		{ 
		  for (i = 0; i < n; i++) { mpz_clear(mp_b1[i]); } 
		  { XFREE(P); XFREE(rp); XFREE(Pt); XFREE(rpt); XFREE(mp_b1); }
#if HAVE_VERBOSE_MODE && HAVE_TIME_H
		  printf("checking lower n-r rows fails, repeat!\n");
#endif

		  continue; 
		}
	    }
	  { XFREE(Pt); XFREE(rpt); }
	  for (i = 0; i < n; i++) { mpz_clear(mp_b1[i]); } { XFREE(mp_b1); }
	  { XFREE(basiscmb[0]); XFREE(basiscmb[1]); XFREE(basiscmb); }

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
	  tt = (double)(clock()-ti)/CLOCKS_PER_SEC;
	  printf("Solution Checking Time: %f\n", tt);
#endif

	  /* recover the solution vector and the certificate vector */
	  if (r < m) 
	    {
	      /* system has more than one solution, the first case */
	      /* compute Qy */
	      for (i = r; i > 0; i--)
		if (rp[i] != i) { mpz_swap(mp_N[i-1], mp_N[rp[i]-1]); }
	      
	      /* set last n-r entries in z be 0 and compute zP */
	      if (certflag == 1)
		{
		  for (i = r; i < n; i++) { mpz_set_ui(mp_NZ[i], 0); }
		  for (i = r; i > 0; i--)
		    if (P[i] != i) { mpz_swap(mp_NZ[i-1], mp_NZ[P[i]-1]); }
		}
	      { XFREE(P); XFREE(rp); }
	      return 1;
	    }
	  else
	    { 
	      /* system has a unique solution, the second case */
	      { XFREE(P); XFREE(rp); }
	      return 2; 
	    }
	}
      else 
	{
	  if ((r == m) && (certflag == 0))
	    {
	      { XFREE(basiscmb[0]); XFREE(basiscmb[1]); XFREE(basiscmb); }
	      { XFREE(P); XFREE(rp); XFREE(Pt); XFREE(rpt); }
	      return 3;
	    }
	  P[r+1] = idx+1;      /* update P for permutation of row r+1 */

	  /* compute u = A_21.A_11^(-1) */
	  mpz_init(mp_Du);
	  mp_Nu = XMALLOC(mpz_t, r);
	  for (i = 0; i < r; i++) { mpz_init(mp_Nu[i]); }
	  A11RNS = XMALLOC(Double *, basislen);
	  for (i = 0; i < basislen; i++)
	    {
	      A11RNS[i] = XMALLOC(Double, r*r);
	      for (j = 0; j < r; j++)
		for (l = 0; l < r; l++)
		  A11RNS[i][j*r+l] = \
		    (Double)mpz_fdiv_ui(mp_A[Pt[j]*m+rpt[l]], basis[i]);
	    }
	  mp_A21 = XMALLOC(mpz_t, r);
	  for (i = 0; i < r; i++)
	    mpz_init_set(mp_A21[i], mp_A[Pt[idx]*m+rpt[i]]);
	  nonsingSolvRNSMM(LeftSolu, r, 1, basislen, basis, A11RNS, mp_A21, \
			   mp_Nu, mp_Du);	  
	  for (i = 0; i < basislen; i++)  { XFREE(A11RNS[i]); } 
	  for (i = 0; i < r; i++) { mpz_clear(mp_A21[i]); }
	  { XFREE(mp_A21); XFREE(A11RNS); }
	  if (r < m)
	    {
	      /* compare u.A_12 with A_22 */
	      A12RNS = XMALLOC(Double *, basislen);
	      for (i = 0; i < basislen; i++)
		{
		  A12RNS[i] = XMALLOC(Double, r*(m-r));
		  for (j = 0; j < r; j++)
		    for (l = 0; l < m-r; l++)
		      A12RNS[i][j*(m-r)+l] = \
			(Double)mpz_fdiv_ui(mp_A[Pt[j]*m+rpt[r+l]], basis[i]);
		}
	      mp_A22 = XMALLOC(mpz_t, m-r);
	      for (i = 0; i < m-r; i++)
		mpz_init_set(mp_A22[i], mp_A[Pt[idx]*m+rpt[r+i]]);
	      cmp = LVecSMatMulCmp(LeftMul, basislen, r, m-r, basis, A12RNS, \
				   mp_Du, mp_Nu, mp_A22);
	      for (i = 0; i < basislen; i++) { XFREE(A12RNS[i]); } 
	      for (i = 0; i < m-r; i++) { mpz_clear(mp_A22[i]); }
	      { XFREE(A12RNS); XFREE(mp_A22); }

	      /* comparison fails, restart */
	      if (cmp == 0) 
		{ 
		  mpz_clear(mp_Du);
		  for (i = 0; i < r; i++) { mpz_clear(mp_Nu[i]); }
		  { XFREE(P); XFREE(rp); XFREE(Pt); XFREE(rpt); XFREE(mp_Nu); }

#if HAVE_VERBOSE_MODE && HAVE_TIME_H
		  printf("checking no solution case fails, repeat!\n");
#endif

		  continue;
		}
	    }
	  if (certflag == 1)
	    {
	      /* set q = (u, -1, 0, ..., 0), which is store in mp_NZ */
	      mpz_set(mp_DZ, mp_Du); 
	      for (i = 0; i < r; i++) { mpz_set(mp_NZ[i], mp_Nu[i]); }
	      mpz_set(mp_NZ[r], mp_Du);
	      mpz_neg(mp_NZ[r], mp_NZ[r]);
	      for (i = r+1; i < n; i++) { mpz_set_ui(mp_NZ[i], 0); }

	      /* compute qP */
	      for (i = r+1; i > 0; i--)
		if (P[i] != i) { mpz_swap(mp_NZ[i-1], mp_NZ[P[i]-1]); }
	    }

	  mpz_clear(mp_Du);
	  { XFREE(basiscmb[0]); XFREE(basiscmb[1]); XFREE(basiscmb); }
	  for (i = 0; i < r; i++) { mpz_clear(mp_Nu[i]); } { XFREE(mp_Nu); }
	  { XFREE(P); XFREE(rp); XFREE(Pt); XFREE(rpt); }

	  /* system has no solution, the third case */
	  return 3;
	}
    }
}





/*
 *
 * Calling Sequence:
 *   At <-- revseq(r, m, A)
 * 
 * Summary:
 *   Perform operations on the permutation vector
 *
 * Description:
 *   Let A be a vector length m+1 to record the permutations over a matrix M. 
 *   A[i] represents the switch of row/column i-1 of M with row/column A[i]-1
 *   of M, i = 1..r. The permutation order is A[r].A[r-1]. ... .A[1].M. 
 *
 *   The function at first generates a vector At length m, At[i] = i, 
 *   i = 0..m-1. Then apply A[r]. ... .A[1] in order on At. The outputed At 
 *   means that row/column At[i] of matrix M will be changed to row/column i 
 *   of M after applying all the permutations to M in order.
 *
 * Input:
 *   r: long, permutation happens in the first r row/column of M
 *   m: long, length of A
 *   A: 1-dim long array length m+1, record of permutations on matrix M
 *
 * Return value:
 *   At: 1-dim long array length m, explained above
 *
 */

long *
revseq(const long r, const long m, const long *A)
{
  long i, t, *At;

  At = XMALLOC(long, m);
  for (i = 0; i < m; i++) { At[i] = i; }
  for (i = 1; i < r+1; i++)
    if (A[i] != i)
      { t = At[i-1];  At[i-1] = At[A[i]-1];  At[A[i]-1] = t; }

  return At;
}




/*
 * Calling Sequence:
 *   compressBoundLong(n, m, Pt, A, mp_bd)
 * 
 * Summary:
 *   Compute the maximum magnitude of a compressed signed long submatrix
 * 
 * Description:
 *   Let A be a n1 x m matrix, B be a n x m submatrix of A. Pt[i] represents 
 *   that row i of B comes from row Pt[i] of A. The function computes the 
 *   maximum magnitude of entries in the compressed matrix B.P, where P is an
 *   m x n+k 0-1 matrix.
 *
 * Input: 
 *      n: long, row dimension of B
 *      m: long, column dimension of B
 *     Pt: 1-dim long array length n1+1, storing row relations between B and A
 *   mp_A: 1-dim long array length n1*m, n1 x m matrix A
 *
 * Output:
 *   mp_bd: mpz_t, maximum magnitude of entries in B.P
 *
 */

void 
compressBoundLong (const long n, const long m, const long *Pt, const long *A,\
		   mpz_t mp_bd)
{
  long i, j, temp;
  mpz_t mp_temp;

  mpz_init(mp_temp);
  mpz_set_ui(mp_bd, 0);
  for (i = 0; i < n; i++) 
    {
      mpz_set_ui(mp_temp, 0);
      for (j = 0; j < m; j++)
	{
	  temp = labs(A[Pt[i]*m+j]);
	  mpz_add_ui(mp_temp, mp_temp, temp);
	}
      if (mpz_cmp(mp_bd, mp_temp) < 0) { mpz_set(mp_bd, mp_temp); }
    }

  mpz_clear(mp_temp);
}



/*
 * Calling Sequence:
 *   compressBoundMP(n, m, Pt, mp_A, mp_bd)
 *
 * Summary: 
 *   Compute the maximum magnitude of a compressed mpz_t submatrix
 * 
 * Description:
 *   Let A be a n1 x m matrix, B be a n x m submatrix of A. Pt[i] represents 
 *   that row i of B comes from row Pt[i] of A. The function computes the 
 *   maximum magnitude of entries in the compressed matrix B.P, where P is an
 *   m x n+k 0-1 matrix.
 *
 * Input: 
 *      n: long, row dimension of B
 *      m: long, column dimension of B
 *     Pt: 1-dim long array length n1+1, storing row relations between B and A
 *   mp_A: 1-dim mpz_t array length n1*m, n1 x m matrix A
 *
 * Output:
 *   mp_bd: mpz_t, maximum magnitude of entries in B.P
 *
 */

void 
compressBoundMP (const long n, const long m, const long *Pt, mpz_t *mp_A, \
		 mpz_t mp_bd)
{
  long i, j;
  mpz_t mp_temp, mp_temp1;

  { mpz_init(mp_temp); mpz_init(mp_temp1); }
  mpz_set_ui(mp_bd, 0);
  for (i = 0; i < n; i++) 
    {
      mpz_set_ui(mp_temp, 0);
      for (j = 0; j < m; j++)
	{
	  mpz_abs(mp_temp1, mp_A[Pt[i]*m+j]);
	  mpz_add(mp_temp, mp_temp, mp_temp1);
	}
      if (mpz_cmp(mp_bd, mp_temp) < 0) { mpz_set(mp_bd, mp_temp); }
    }

  { mpz_clear(mp_temp); mpz_clear(mp_temp1); }
}






/*
 * Calling Sequence:
 *   1/0 <-- LVecSMatMulCmp(mulpos, basislen, n, m, basis, ARNS, mp_s, 
 *                          mp_V, mp_b)
 *
 * Summary:
 *   Verify the equality of a matrix-vector product and a scalar-vector 
 *   product
 *
 * Description:
 *   The function verifies whether a matrix-vector product A.V or V.A is equal
 *   to a scalar-vector product s.b, where A is a n x m matrix, b and V is a 
 *   vector, and s is a scalar. The flag mulpos is used to specify whether the
 *   matrix-vector product is A.V or V.A. The comparison result is output 
 *   sensitive, i.e., the comparsion will stop as long as one failure occurs.
 *   The input matrix A is represented in a RNS.
 *
 * Input:
 *     mulpos: LeftMul/RightMul, flag to indicate whether A.V or V.A is 
 *             performed. 
 *             If mulpos = LeftMul ==> V.A. If mulpos = RightMul ==> A.V
 *   basislen: long, dimension of the RNS basis
 *          n: long, row dimension of A
 *          m: long, column dimension of A
 *      basis: 1-dim FiniteField array length basislen, RNS basis
 *       ARNS: 2-dim Double array, dimension basislen x n*m, representation of
 *             A in the RNS. ARNS[i] = A mod basis[i], i = 0..basislen-1
 *       mp_s: mpz_t, scalar s
 *       mp_V: 1-dim mpz_t array length n or m depending on mulpos, vector V
 *       mp_b: 1-dim mpz_t array length n or m depending on mulpos, vector b
 *
 * Return:
 *   1: comparison succeeds 
 *   0: comparison fails
 *     
 */

long 
LVecSMatMulCmp (const enum MULT_POS mulpos, const long basislen, \
		const long n, const long m, const FiniteField *basis, \
		Double **ARNS, mpz_t mp_s, mpz_t *mp_V, mpz_t *mp_b)
{
  long i, j, l, sharednum, unsharednum;
  FiniteField *bdcoeff, *cmbasis;
  double temp;
  mpz_t mp_AV, mp_temp, mp_AV1;
  Double **U;

  if (mulpos == LeftMul) { sharednum = n;  unsharednum = m; }
  else if (mulpos == RightMul) { sharednum = m;  unsharednum = n; }

  /* allocate matrix U to store mix radix coefficients of one row/column 
     of matrix A */
  U = XMALLOC(Double *, basislen);
  for (i = 0; i < basislen; i++) 
    U[i] = XMALLOC(Double, sharednum);
  cmbasis = combBasis(basislen, basis);
  bdcoeff = repBound(basislen, basis, cmbasis);
  mpz_init(mp_AV);
  mpz_init(mp_AV1);
  mpz_init(mp_temp);
  for (i = 0; i < unsharednum; i++)
    {
      /* apply Garne