rnsop.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 "RNSop.h"

/*
 *
 * Calling Sequence:
 *   basisExt(len, n, p, basis, cmbasis, cumprod, bdcoeff, R, RE)
 *
 * Summary:
 *   Given a representation of a matrix/vector in some RNS, extend to compute 
 *   the representation in another positive integer
 *
 * Description:
 *   Let R be the representation of a matrix/vector M in a residue basis
 *   'basis', i.e., R[i] = mod(M, basis[i]) (i = 0..len-1). The function 
 *   computes the representation of M in another positive integer p, 
 *   RE = mod(M, p) using Garner's algorithm. 'mod' represents positive modular
 *   operation.
 *
 *   Let q be product of all elements in the RNS basis. Every entry m in M 
 *   satisfies -(q-1)/2 <= m <= (q-1)/2. That is, M has both positive entries
 *   and negative entries.
 *
 *   To avoid repeat computations, the function takes some precomputed 
 *   informations as input, which are listed below.
 *
 * Input: 
 *       len: long, dimension of RNS basis
 *         n: long, length of array RE
 *         p: FiniteField, modulus
 *     basis: 1-dim FiniteField array length len, RNS basis
 *   cmbasis: 1-dim FiniteField array length len, computed by function 
 *            combBasis, inverses of special combination of RNS basis
 *   cumprod: double, computed by function cumProd, 
 *            (-basis[0]*basis[1]*...*basis[len-1]) mod p 
 *   bdcoeff: 1-dim FiniteField array length len, computed by function repBound
 *         R: 1-dim Double array length n*len, representation of a len x n 
 *            matrix, R[i]=mod(M, basis[i]) (i=0..len-1)
 * 
 * Output:
 *   RE: 1-dim Double array length n, RE = mod(M, p), the space of RE 
 *       should be allocated before calling the function.
 *
 * Precondition:
 *   t <= 2^53-1, where t is the maximal intermidiate value arised in this 
 *   function,
 *   t = max(2*(basis[len-1]-1)^2, (p-1)^2+basis[len-1]-1)
 *           
 */

void
basisExt (const long len, const long n, const FiniteField p, \
	  const FiniteField *basis, const FiniteField *cmbasis, \
	  const double cumprod, const FiniteField *bdcoeff, Double **R, \
	  Double *RE)
{
  long i, j;
  double temp;
  Double **U;
  const FiniteField *q, *qinv;

  q = basis;
  qinv = cmbasis;

  /* if p = q[i] then just copy the corresponding column to RE */
  for (i = 0; i < len; i++) 
    { if (p == q[i]) { cblas_dcopy(n, R[i], 1, RE, 1); return; } }
  U = XMALLOC(Double *, len);
  for (i = 0; i < len; i++) { U[i] = XMALLOC(Double, n); }

  /* compute the coefficients of mix radix in positive representation by 
     inplacing modular matrix U */
  cblas_dcopy(n, R[0], 1, U[0], 1);
  for (i = 1; i < len; i++)
    {
      cblas_dcopy(n, U[i-1], 1, U[i], 1);
      for (j = i-2; j >= 0; j--)
	{
	  catlas_daxpby(n, 1.0, U[j], 1, (double)(q[j] % q[i]), U[i], 1);
	  Dmod((double)q[i], U[i], 1, n, n);
	}
      temp = (double)qinv[i]*(double)(q[i]-1);
      temp = fmod(temp, (double)q[i]);
      catlas_daxpby(n, (double)qinv[i], R[i], 1, temp, U[i], 1);
      Dmod((double)q[i], U[i], 1, n, n);
    }

  /* compute mod(r, p) in positive representation and store into RE */
  cblas_dcopy(n, U[len-1], 1, RE, 1);
  Dmod((double)p, RE, 1, n, n);
  for (i = len-2; i >= 0; i--)
    {
      catlas_daxpby(n, 1.0, U[i], 1, (double)(q[i] % p), RE, 1);
      Dmod((double)p, RE, 1, n, n);
    }

  /* convert to symmetric representation */
  for (i = 0; i < n; i++)
    for (j = len-1; j >= 0; j--)
      {
	if (U[j][i] > bdcoeff[j])
	  {
	    RE[i] = fmod(RE[i]+cumprod, (double)p);
	    break;
	  }
	else if (U[j][i] < bdcoeff[j]) { break; }
      }
  for (i = 0; i < len; i++) { XFREE(U[i]); } { XFREE(U); }

  return;
}


/*
 *
 * Calling Sequence:
 *   basisExtPos(len, n, p, basis, cmbasis, R, RE)
 *
 * Summary:
 *   Given a representation of a non-negative matrix/vector in some RNS, 
 *   extend to compute the representation in another positive integer
 *
 * Description:
 *   Let R be the representation of a matrix/vector M in a residue basis
 *   'basis', i.e., R[i] = mod(M, basis[i]) (i = 0..len-1). The function 
 *   computes the representation of M in another positive integer p, 
 *   RE = mod(M, p) using Garner's algorithm. 'mod' represents positive modular
 *   operation.
 *
 *   Let q be product of all elements in the RNS basis. Every entry m in M 
 *   satisfies 0 <= m <= q-1. That is, M only contains non-negative entries.
 *
 *   To avoid repeat computations, the function takes some precomputed 
 *   informations as input, which are listed below.
 *
 * Input:
 *       len: long, dimension of RNS basis
 *         n: long, length of array RE
 *         p: FiniteField, modulus
 *     basis: 1-dim FiniteField array length len, RNS basis
 *   cmbasis: 1-dim FiniteField array length len, computed by function 
 *            combBasis, inverses of special combination of RNS basis
 *         R: 1-dim Double array length n*len, representation of a len x n 
 *            matrix, R[i]=mod(M, basis[i]) (i=0..len-1)
 * 
 * Output:
 *   RE: 1-dim Double array length n, RE = mod(M, p), the space of RE 
 *       should be allocated before calling the function.
 *
 * Precondition: 
 *   t <= 2^53-1, where t is the maximal intermidiate value arised in this 
 *   function, t = max(2*(basis[len-1]-1)^2, (p-1)^2+basis[len-1]-1)
 *           
 */

void
basisExtPos (const long len, const long n, const FiniteField p, \
	     const FiniteField *basis, const FiniteField *cmbasis, \
	     Double **R, Double *RE)
{
  long i, j;
  double temp;
  Double **U;
  const FiniteField *q, *qinv;

  q = basis;
  qinv = cmbasis;

  /* if p = q[i] then just copy the corresponding column to RE */
  for (i = 0; i < len; i++) 
    { if (p == q[i]) { cblas_dcopy(n, R[i], 1, RE, 1); return; } }
  U = XMALLOC(Double *, len);
  for (i = 0; i < len; i++) { U[i] = XMALLOC(Double, n); }

  /* compute the coefficients of mix radix in positive representation by 
     inplacing modular matrix U */
  cblas_dcopy(n, R[0], 1, U[0], 1);
  for (i = 1; i < len; i++)
    {
      cblas_dcopy(n, U[i-1], 1, U[i], 1);
      for (j = i-2; j >= 0; j--)
	{
	  catlas_daxpby(n, 1.0, U[j], 1, (double)(q[j] % q[i]), U[i], 1);
	  Dmod((double)q[i], U[i], 1, n, n);
	}
      temp = (double)qinv[i]*(double)(q[i]-1);
      temp = fmod(temp, (double)q[i]);
      catlas_daxpby(n, (double)qinv[i], R[i], 1, temp, U[i], 1);
      Dmod((double)q[i], U[i], 1, n, n);
    }

  /* compute mod(r, p) in positive representation and store into RE */
  cblas_dcopy(n, U[len-1], 1, RE, 1);
  Dmod((double)p, RE, 1, n, n);
  for (i = len-2; i >= 0; i--)
    {
      catlas_daxpby(n, 1.0, U[i], 1, (double)(q[i] % p), RE, 1);
      Dmod((double)p, RE, 1, n, n);
    }
  for (i = 0; i < len; i++) { XFREE(U[i]); } { XFREE(U); }

  return;
}



/*
 *
 * Calling Sequence:
 *   basisProd(len, basis, mp_prod)
 *
 * Summary:
 *   Compute the product of elements of a RNS basis
 *
 * Description:
 *   Let a RNS basis be 'basis'. The function computes the product of its 
 *   elements basis[0]*basis[1]*...*basis[len-1].
 *
 * Input:
 *     len: long, dimension of RNS basis
 *   basis: 1-dim FiniteField array length len, RNS basis
 * 
 * Output:
 *   mp_prod: mpz_t, product of elements in 'basis'
 *
 */

void
basisProd (const long len, const FiniteField *basis, mpz_t mp_prod)
{
  long i;

  mpz_set_ui(mp_prod, basis[0]);
  for (i = 1; i < len; i++) { mpz_mul_ui(mp_prod, mp_prod, basis[i]); }
}



/*
 *
 * Calling Sequence:
 *   ChineseRemainder(len, mp_prod, basis, cmbasis, bdcoeff, Ac, mp_Ac)
 *
 * Summary:
 *   Given a representation of an integer in some RNS, use Chinese Remainder 
 *   Algorithm to reconstruct the integer
 * 
 * Description:
 *   Let A be an integer, and Ac contains the representation of A in a RNS
 *   basis 'basis', i.e. Ac[i] = mod(A, basis[i]), (i = 0..len). Here 'mod' 
 *   is in positive representation. The function reconstructs the integer A 
 *   given the RNS basis 'basis' and Ac. 
 *
 *   To avoid repeat computations, the function takes some precomputed 
 *   informations as input, which are listed below.
 *
 * Input: 
 *       len: long, dimension of RNS basis
 *   mp_prod: mpz_t, computed by function basisProd, product of RNS basis
 *     basis: 1-dim FiniteField array length len, RNS basis
 *   cmbasis: 1-dim FiniteField array length len, computed by function 
 *            combBasis, inverses of special combination of RNS basis
 *   bdcoeff: 1-dim FiniteField array length len, computed by function repBound
 *        Ac: 1-dim Double array length n, representation of A in RNS 
 *
 * Output:
 *   mp_Ac: mpz_t, reconstructed integer A
 *
 * Precondition:
 *   Let q be product of all elements in the RNS basis. Then A must satisfy 
 *   -(q-1)/2 <= A <= (q-1)/2. 
 *
 */

void
ChineseRemainder (const long len, const mpz_t mp_prod, \
		  const FiniteField *basis, const FiniteField *cmbasis, \
		  const FiniteField *bdcoeff, Double *Ac, mpz_t mp_Ac) 
{
  long i, j;
  double temp, tempq, tempqinv;
  Double *U;

  U = XMALLOC(Double, len);

  /* compute the coefficients of mix radix in positive representation by 
     inplacing modular matrix U */
  U[0] = Ac[0];
  for (i = 1; i < len; i++)
    {
      U[i] = U[i-1];
      tempq = (double)basis[i];
      tempqinv = (double)cmbasis[i];
      for (j = i-2; j >= 0; j--)
	{
	  U[i] = U[j] + U[i]*fmod((double)basis[j], tempq);
	  U[i] = fmod(U[i], tempq);
	}
      temp = fmod(tempqinv*(double)(basis[i]-1), tempq);
      U[i] = fmod(tempqinv*Ac[i]+temp*U[i], tempq);
    }
  /* compute Ac in positive representation */
  mpz_set_d(mp_Ac, U[len-1]); 
  for (j = len-2; j >= 0; j--)
    {
      mpz_mul_ui(mp_Ac, mp_Ac, basis[j]);
      mpz_add_ui(mp_Ac, mp_Ac, (FiniteField)U[j]);
    }
  /* transfer from positive representation to symmetric representation */
  for (j = len-1; j >= 0; j--)
    {
      if (U[j] > bdcoeff[j])
	{
	  mpz_sub(mp_Ac, mp_Ac, mp_prod);
	  break;
	}
      else if (U[j] < bdcoeff[j]) { break; }
    }
  XFREE(U);

  return;
}



/*
 *
 * Calling Sequence:
 *   ChineseRemainderPos(len, basis, cmbasis, Ac, mp_Ac)
 *
 * Summary:
 *   Given a representation of a non-negative integer in some RNS, use Chinese 
 *   Remainder Algorithm to reconstruct the integer
 *
 * Description:
 *   Let A be a non-negative integer, and Ac contains the representation of A
 *   in a RNS basis 'basis', i.e. Ac[i] = mod(A, basis[i]), (i = 0..len). 
 *   Here 'mod' is in positive representation. The function reconstructs the 
 *   integer A given the RNS basis 'basis' and Ac. 
 *
 *   To avoid repeat computations, the function takes some precomputed 
 *   informations as input, which are listed below.
 *
 * Input: 
 *       len: long, dimension of RNS basis
 *     basis: 1-dim FiniteField array length len, RNS basis
 *   cmbasis: 1-dim FiniteField array length len, computed by function 
 *            combBasis, inverses of special combination of RNS basis
 *        Ac: 1-dim Double array length n, representation of A in RNS 
 *
 * Output:
 *   mp_Ac: mpz_t, reconstructed integer A
 *
 * Precondition:
 *   Let q be product of all elements in the RNS basis. Then A must satisfy 
 *   0 <= A <= q-1. 
 *
 */

void
ChineseRemainderPos (const long len, const FiniteField *basis, \
		     const FiniteField *cmbasis, Double *Ac, mpz_t mp_Ac) 
{
  long i, j;
  double temp, tempq, tempqinv;
  Double *U;

  U = XMALLOC(Double, len);

  /* compute the coefficients of mix radix in positive representation by 
     inplacing modular matrix U */
  U[0] = Ac[0];