padiclift.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 "padiclift.h"

/*
 * Calling Sequence:
 *   -1/i <-- liftInit(liftbasislen, liftbasis, n, A, mp_basisprod, 
 *            mp_extbasisprod, extbasislen, cmbasis, extbdcoeff, 
 *            liftbasisInv, AInv, extbasis, ARNS)
 *
 * Summary:
 *   Perform initialization operations before lifting, where the left hand 
 *   side input matrix is a signed long matrix
 *
 * Description:
 *   Initialization step before lifting computes necessory data and stores
 *   them to avoid recomputations if continuous lifting is needed. Seperating
 *   initialization step and lifting step makes it possible to lift adaptively.
 *   Pointers are passed into the calling sequence to store the outputs of 
 *   initialization.
 *
 * Input:
 *   liftbasislen: long, dimension of lifting basis
 *      liftbasis: 1-dim FiniteField array length liftbasislen, lifting basis
 *              n: long, dimension of input matrix A
 *              A: 1-dim long array length n*n, representation of n x n input
 *                 matrix
 *
 * Return:
 *   The first index i such that A^(-1) mod liftbasis[i] doesn't exist, where 
 *   i starts from 0. Otherwise, return -1 if A^(-1) mod liftbasis[i] exists 
 *   for any i
 *
 * Output: 
 *     mp_basisiprod: mpz_t, product of lifting basis
 *   mp_extbasisprod: mpz_t, product of extended RNS basis 
 *       extbasislen: pointer to long int, storing the dimension of extended 
 *                    RNS basis. Extended basis is used to compute the 
 *                    matrix-vector product AC_i during lifting
 *           cmbasis: pointer to a 1-dim FiniteField array, storing the 
 *                    special combination of lifting basis computed by 
 *                    function combBasis
 *        extbdcoeff: pointer to a 1-dim FiniteField array, storing the mix 
 *                    radix coefficients of a special integer, computed by 
 *                    function repBound, in extended RNS basis 
 *      liftbasisInv: pointer to a 1-dim Double array, storing 
 *                    (1/mp_basisprod) mod extbasis[i]
 *              AInv: pointer to a 1-dim Double array, storing the modp matrix
 *                    A^(-1) mod liftbasis[i] 
 *          extbasis: pointer to a 2-dim FiniteField array, where
 *                  - (*extbasis)[0] = extended RNS basis
 *                  - (*extbasis)[1] = the special combination of extended RNS
 *                    basis computed by function combBasis
 *              ARNS: pointer to a 2-dim Double array, where the last 
 *                    dimension (*ARNS)[i] stores A mod ith element in 
 *                    extended RNS basis
 */

long
liftInit (const long liftbasislen, const FiniteField *liftbasis, \
	  const long n, const long *A, mpz_t mp_basisprod, \
	  mpz_t mp_extbasisprod, long *extbasislen, FiniteField **cmbasis, \
	  FiniteField **extbdcoeff, Double **liftbasisInv, Double **AInv, \
	  FiniteField ***extbasis, Double ***ARNS)
{
  long i, j, alpha, minv, p, temp, len=0;
  mpz_t mp_maxInter, mp_alpha, mp_temp;
  FiniteField *q, *qinv;

  for (i = 0; i < liftbasislen; i++)
    {
      p = (long)liftbasis[i]; 
      for (j = 0; j < n*n; j++) 
	AInv[i][j] = (double)((temp = (A[j] % p)) >= 0 \
			      ? temp : p+temp);
      minv = mInverse(liftbasis[i], AInv[i], n);

      /* if fail to find inverse of A mod liftbasis[i] */
      if (minv == 0) { return i; }
    }
  *cmbasis = combBasis(liftbasislen, liftbasis);
  basisProd(liftbasislen, liftbasis, mp_basisprod);

  /* compute maximum intermediate result mp_maxInter */
  alpha = maxMagnLong(A, n, n, n);
  mpz_init_set_ui(mp_alpha, alpha);
  mpz_init(mp_maxInter);
  maxExtInter(mp_alpha, n, mp_maxInter);
  mpz_clear(mp_alpha);

  *extbasis = findRNS(liftbasis[liftbasislen-1]-1, mp_maxInter, &len);
  mpz_clear(mp_maxInter);
  *extbasislen = len;
  q = *(*extbasis);
  qinv = *((*extbasis)+1);
  *liftbasisInv = invBasis(len, q, mp_basisprod);
  basisProd(len, q, mp_extbasisprod);
  *extbdcoeff = repBound(len, q, qinv);
  *ARNS = XMALLOC(Double *, len);
  for (i = 0; i < len; i++)
    {
      p = (long)q[i];
      (*ARNS)[i] = XMALLOC(Double, n*n);
      for (j = 0; j < n*n; j++) 
	(*ARNS)[i][j] = (double)((temp = (A[j] % p)) >= 0 ? \
				 temp : p+temp);
    }
  return -1;
}




/*
 * Calling Sequence:
 *   -1/i <-- liftInitLlhs(liftbasislen, liftbasis, n, mp_A, mp_basisprod, 
 *                         mp_extbasisprod, extbasislen, cmbasis, extbdcoeff, 
 *                         liftbasisInv, AInv, extbasis, ARNS)
 *
 * Summary:
 *   Perform initialization operations before lifting, where the left hand 
 *   side input matrix is a mpz_t matrix
 *
 * Description:
 *   Initialization step before lifting computes necessory data and stores
 *   them to avoid recomputations if continuous lifting is needed. Seperating
 *   initialization step and lifting step makes it possible to lift adaptively.
 *   Pointers are passed into the calling sequence to store the outputs of 
 *   initialization.
 *
 * Input:
 *   liftbasislen: long, dimension of lifting basis
 *      liftbasis: 1-dim FiniteField array length liftbasislen, lifting basis
 *              n: long, dimension of input matrix A
 *           mp_A: 1-dim mpz_t array length n*n, representation of n x n input
 *                 matrix
 *
 * Return:
 *   The first index i such that A^(-1) mod liftbasis[i] doesn't exist, where 
 *   i starts from 0. Otherwise, return -1 if A^(-1) mod liftbasis[i] exists 
 *   for any i
 *
 * Output: 
 *     mp_basisiprod: mpz_t, product of lifting basis
 *   mp_extbasisprod: mpz_t, product of extended RNS basis 
 *       extbasislen: pointer to long int, storing the dimension of extended 
 *                    RNS basis. Extended basis is used to compute the 
 *                    matrix-vector product AC_i during lifting
 *           cmbasis: pointer to a 1-dim FiniteField array, storing the 
 *                    special combination of lifting basis computed by 
 *                    function combBasis
 *        extbdcoeff: pointer to a 1-dim FiniteField array, storing the mix 
 *                    radix coefficients of a special integer, computed by 
 *                    function repBound, in extended RNS basis 
 *      liftbasisInv: pointer to a 1-dim Double array, storing 
 *                    (1/mp_basisprod) mod extbasis[i]
 *              AInv: pointer to a 1-dim Double array, storing the modp matrix
 *                    A^(-1) mod liftbasis[i] 
 *          extbasis: pointer to a 2-dim FiniteField array, where
 *                  - (*extbasis)[0] = extended RNS basis
 *                  - (*extbasis)[1] = the special combination of extended RNS
 *                    basis computed by function combBasis
 *              ARNS: pointer to a 2-dim Double array, where the last 
 *                    dimension (*ARNS)[i] stores A mod ith element in 
 *                    extended RNS basis
 */

long
liftInitLlhs (const long liftbasislen, const FiniteField *liftbasis, \
	      const long n, mpz_t *mp_A, mpz_t mp_basisprod, \
	      mpz_t mp_extbasisprod, long *extbasislen, \
	      FiniteField **cmbasis, FiniteField **extbdcoeff, \
	      Double **liftbasisInv, Double **AInv, FiniteField ***extbasis, \
	      Double ***ARNS)
{
  long i, j, minv, len=0;
  double dtemp;
  mpz_t mp_maxInter, mp_alpha;
  FiniteField *q, *qinv;

  for (i = 0; i < liftbasislen; i++)
    {
      for (j = 0; j < n*n; j++) 
	AInv[i][j] = (Double)mpz_fdiv_ui(mp_A[j], liftbasis[i]);
      minv = mInverse(liftbasis[i], AInv[i], n);

      /* if fail to find inverse of A mod liftbasis[i] */
      if (minv == 0) { return i; }
    }
  *cmbasis = combBasis(liftbasislen, liftbasis);
  basisProd(liftbasislen, liftbasis, mp_basisprod);

  /* compute maximum intermediate result mp_maxInter */
  mpz_init(mp_alpha);
  maxMagnMP(mp_A, n, n, n, mp_alpha);
  mpz_init(mp_maxInter);
  maxExtInter(mp_alpha, n, mp_maxInter);
  mpz_clear(mp_alpha);

  *extbasis = findRNS(liftbasis[liftbasislen-1]-1, mp_maxInter, &len);
  mpz_clear(mp_maxInter);
  *extbasislen = len;
  q = *(*extbasis);
  qinv = *((*extbasis)+1);
  *liftbasisInv = invBasis(len, q, mp_basisprod);
  basisProd(len, q, mp_extbasisprod);
  *extbdcoeff = repBound(len, q, qinv);
  *ARNS = XMALLOC(Double *, len);
  for (i = 0; i < len; i++)
    {
      (*ARNS)[i] = XMALLOC(Double, n*n);
      for (j = 0; j < n*n; j++) 
	(*ARNS)[i][j] = (Double)mpz_fdiv_ui(mp_A[j], q[i]);
    }

  return -1;
}


/*
 * Calling Sequence:
 * -1/i <-- liftInitRNS(liftbasislen, liftbasis, basislen, basis, n, ARNS, 
 *          mp_liftbasisprod, mp_extbasisprod, extbasislen, cmliftbasis, 
 *	    liftbasisInv, extbdcoeff, AInv, extbasis, AExtRNS)
 *
 * Summary:
 *   Perform initialization operations before lifting where the left hand 
 *   side input matrix is represented in a RNS
 *   
 * Description:
 *   Initialization step before lifting computes necessory data and stores
 *   them to avoid recomputations if continuous lifting is needed. Seperating
 *   initialization step and lifting step makes it possible to lift adaptively.
 *   Pointers are passed into the calling sequence to store the outputs of 
 *   initialization.
 *
 * Input:
 *   liftbasislen: long, dimension of lifting basis
 *      liftbasis: 1-dim FiniteField array length liftbasislen, lifting basis
 *       basislen: long, dimension of RNS basis used to represent A
 *          basis: 1-dim FiniteField array length basislen, RNS basis used to
 *                 represent A
 *              n: long, dimension of A
 *           ARNS: 2-dim Double array, dimension basislen x n^2,
 *                 representation of A in RNS, ARNS[i] = A mod basis[i]
 *
 * Return:
 *   The first index i such that A^(-1) mod liftbasis[i] doesn't exist, where