nonsysolve.c

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

/* ---------------------------------------------------------------------
 *
 * -- Integer Matrix Library (IML)
 *    (C) Copyright 2004 All Rights Reserved
 *
 * -- IML routines -- Version 0.1.0 -- September, 2004
 *
 * Author         : Zhuliang Chen
 * Contributor(s) : Arne Storjohann
 * University of Waterloo -- School of Computer Science
 * Waterloo, Ontario, N2L3G1 Canada
 *
 * ---------------------------------------------------------------------
 *
 * -- Copyright notice and Licensing terms:
 *
 *  Redistribution  and  use in  source and binary forms, with or without
 *  modification, are  permitted provided  that the following  conditions
 *  are met:
 *
 * 1. Redistributions  of  source  code  must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce  the above copyright
 *    notice,  this list of conditions, and the  following disclaimer in
 *    the documentation and/or other materials provided with the distri-
 *    bution.
 * 3. The name of the University,  the IML group,  or the names of its
 *    contributors  may not be used to endorse or promote products deri-
 *    ved from this software without specific written permission.
 *
 * -- Disclaimer:
 *
 * THIS  SOFTWARE  IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,  INCLUDING,  BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,  INDIRECT, INCIDENTAL, SPE-
 * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
 * TO,  PROCUREMENT  OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
 * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  (IN-
 * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */


#include "nonsysolve.h"

/*
 * Calling Sequence:
 *   nonsingSolvMM(solupos, n, m, A, mp_B, mp_N, mp_D)
 *
 * Summary:
 *   Solve nonsingular system of linear equations, where the left hand side 
 *   input matrix is a signed long matrix.
 *
 * Description:
 *   Given a n x n nonsingular signed long matrix A, 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
 *         n: long, dimension of A
 *         m: long, column or row dimension of mp_B depending on solupos
 *         A: 1-dim signed long array length n*n, representing the n x n
 *            left hand side input matrix
 *      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.
 *
 * 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
nonsingSolvMM (const enum SOLU_POS solupos, const long n, const long m, \
	       const long *A, mpz_t *mp_B, mpz_t *mp_N, mpz_t mp_D)
{
  long alpha, 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; 
  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_temp); }
  { mpz_init(mp_basisprod); mpz_init(mp_extbasisprod); }

  /* find lifting basis such that the product of the basis elements is 
     approximately n*alpha */
  alpha = maxMagnLong(A, n, n, n);
  mpz_set_ui(mp_alpha, alpha);
  liftbasis = findLiftbasisSmall(n, mp_alpha, &basislen);

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

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

  minv = liftInit(basislen, liftbasis, n, 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 = liftInit(basislen, liftbasis, n, A, mp_basisprod, \
		      mp_extbasisprod, &extbasislen, &cmbasis, &extbdcoeff, \
		      &liftbasisInv, AInv, &extbasis, &ARNS);
    }
  printf("liftinginit2\n");
#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", basislen);
  for (i = 0; i < basislen; 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_basisprod, 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, basislen, extbasislen, mp_basisprod, \
	      mp_extbasisprod, liftbasis, cmbasis, extbdcoeff, liftbasisInv, \
	      mp_r, extbasis, AInv, ARNS);

#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, basislen, n, m, mp_basisprod, liftbasis, \
		   cmbasis, 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 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_si(mp_nb, mp_nb, 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);

#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 < 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:
 *   nonsingSolvLlhsMM(solupos, n, m, mp_A, mp_B, mp_N, mp_D)
 *
 * Summary:
 *   Solve nonsingular system of linear equations, where the left hand side 
 *   input matrix is a mpz_t matrix.
 *
 * Description:
 *   Given a n x n nonsingular mpz_t matrix A, a n x m or m x n mpz_t 
 *   matrix mp_B, this function will compute the solution of the system
 *   1. (mp_A)X = mp_B 
 *   2. X(mp_A) = 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 mp_Amp_N = mp_D*mp_B or mp_Nmp_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
 *         n: long, dimension of A
 *         m: long, column or row dimension of mp_B depending on solupos
 *      mp_A: 1-dim mpz_t array length n*n, representing the n x n left hand
 *            side input matrix
 *      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: 
 *   mp_A must be a nonsingular matrix.
 *
 * Note:
 *   - It is necessary to make sure the input parameters are correct,
 *     expecially the dimension, since there is no parameter checks in the