rnsop.c

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

  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]);
    }
  XFREE(U);

  return;
}



/*
 *
 * Calling Sequence:
 *   cmbasis <-- combBasis(basislen, basis)
 * 
 * Summary:
 *   Compute the special combination of a RNS basis
 * 
 * Description:
 *   Let 'basis' be RNS basis. The function computes an array cmbasis 
 *   satisfying 
 *   cmbasis[0] = 0, cmbasis[i] = mod(1/(basis[0]*...*basis[i-1]), basis[i])
 *                   (i = 1..basislen-1)
 * 
 * Input:
 *   basislen: long, dimension of RNS basis
 *      basis: 1-dim FiniteField array length basislen, RNS basis
 * 
 * Return:
 *   cmbasis: 1-dim FiniteField array length basislen, shown as above
 *
 */

FiniteField *
combBasis (const long basislen, const FiniteField *basis)
{
  long i, j;
  double dtemp;
  mpz_t mp_prod, mp_q;
  FiniteField *cmbasis;

  cmbasis = XMALLOC(FiniteField, basislen);
  cmbasis[0] = 0;
  mpz_init(mp_prod);
  mpz_init(mp_q);
  for (i = 1; i < basislen; i++)
    {
      dtemp = fmod((double)basis[0], (double)basis[i]);
      for (j = 1; j <= i-1; j++)
	dtemp = fmod(dtemp*(double)basis[j], (double)basis[i]);
      mpz_set_ui(mp_q, basis[i]);
      mpz_set_d(mp_prod, dtemp);
      mpz_invert(mp_prod, mp_prod, mp_q);
      cmbasis[i] = mpz_get_ui(mp_prod);
    }
  mpz_clear(mp_prod);
  mpz_clear(mp_q);

  return cmbasis;
}



/*
 *
 * Calling Sequence:
 *   cumprod <-- cumProd(basislen, basis, extbasislen, extbasis)
 *
 * Summary:
 *   Compute the representation of the combination of elements of one RNS basis
 *   in another RNS basis
 * 
 * Description:
 *   Let 'basis' be one RNS basis with dimension basislen, and 'extbasis' be 
 *   another RNS basis with dimension extbasislen. The function computes an 
 *   array cumprod length extbasislen satisfying
 *   cumprod[i] = modp(-basis[0]*...*basis[basislen-1], extbasis[i]),
 *   i = 0..extbasislen-1
 *
 * Input: 
 *      basislen: long, dimension of RNS basis 'basis'
 *         basis: 1-dim FiniteField array length basislen, one RNS basis
 *   extbasislen: long, dimension of RNS basis 'extbasis'
 *      extbasis: 1-dim FiniteField array length basislen, another RNS basis
 * 
 * Return: 
 *   cumprod: 1-dim double array length extbasislen, shown above
 *
 */

double *
cumProd (const long basislen, const FiniteField *basis, \
	 const long extbasislen, const FiniteField *extbasis)
{
  long i, j;
  double dtemp, dextbasis;
  double *cumprod;

  cumprod = XMALLOC(double, extbasislen);
  for (i = 0; i < extbasislen; i++)
    {
      dextbasis = (double)extbasis[i];
      cumprod[i] = fmod((double)basis[0], dextbasis);
      for (j = 1; j < basislen; j++)
	{
	  dtemp = fmod((double)basis[j], dextbasis);
	  cumprod[i] = fmod(cumprod[i]*dtemp, dextbasis);
	}
      cumprod[i] = dextbasis-cumprod[i];
    }

  return cumprod;
}


/*
 *
 * Calling Sequence:
 *   basiscmb <-- findRNS(RNS_bound, mp_maxInter, len)
 *
 * Summary:
 *   Find a RNS basis and its special combination
 *   
 * Description:
 *   Given RNS_bound, the upper bound of the RNS basis, and mp_maxInter, the
 *   function finds a best RNS basis and a combination of that basis.
 *   
 *   The RNS basis 'basis' has the property:
 *   - its elements are all primes
 *   - basis[0] is the largest prime among all the primes at most RNS_bound
 *   - basis[i+1] is the next prime smaller than basis[i] (i = 0..len-2)
 *   - basis[0]*basis[1]*...*basis[len-1] >= mp_maxInter
 *
 *   After finding 'basis', the functions also computes the combination of
 *   'basis' as the operations in function combBasis.
 *
 * Input:  
 *     RNS_bound: FiniteField, the upper bound of the RNS basis
 *   mp_maxInter: mpz_t, the lower bound for the product of elements of basis
 *
 * Return:
 *   basiscmb: 2-dim FiniteField array, dimension 2 x len, where
 *           - basiscmb[0] represents the RNS basis
 *           - basiscmb[1] represents the special combination of basis 
 *
 * Output:
 *   len: pointer to a long int, storing the dimension of the computed
 *        RNS basis
 *
 */

FiniteField **
findRNS (const FiniteField RNS_bound, const mpz_t mp_maxInter, long *length)
{
  long i, j, len=0;
  double prod;
  mpz_t mp_l, mp_prod, mp_q;
  FiniteField **qqinv;

  mpz_init_set_ui(mp_prod, 1);
  mpz_init_set_ui(mp_l, RNS_bound);
  qqinv = XMALLOC(FiniteField *, 2);
  qqinv[0] = NULL;
  while (mpz_cmp(mp_maxInter, mp_prod) > 0)
    {
      ++len;
      qqinv[0] = XREALLOC(FiniteField, qqinv[0], len);
      while (mpz_probab_prime_p(mp_l, 10) == 0) { mpz_sub_ui(mp_l, mp_l, 1); }
      qqinv[0][len-1] = mpz_get_ui(mp_l);
      mpz_sub_ui(mp_l, mp_l, 1);
      mpz_mul_ui(mp_prod, mp_prod, qqinv[0][len-1]);
    }
  mpz_clear(mp_prod);
  mpz_clear(mp_l);
  qqinv[1] = XMALLOC(FiniteField, len);
  qqinv[1][0] = 0;
  mpz_init(mp_prod);
  mpz_init(mp_q);
  for (i = 1; i < len; i++)
    {  
      prod = (double)(qqinv[0][0] % qqinv[0][i]);
      for (j = 1; j <= i-1; j++) 
	prod = fmod(prod*(double)qqinv[0][j], (double)qqinv[0][i]);
      mpz_set_ui(mp_q, qqinv[0][i]);
      mpz_set_d(mp_prod, prod);
      mpz_invert(mp_prod, mp_prod, mp_q);
      qqinv[1][i] = mpz_get_ui(mp_prod);
    }
  mpz_clear(mp_prod);
  mpz_clear(mp_q);  
  *length = len;

  return qqinv;
}



/*
 *
 * Calling Sequence:
 *   maxInter(mp_prod, mp_alpha, n, mp_b)
 *
 * Summary:
 *   Compute the maximum interval of positive and negative results of a
 *   matrix-matrix or matrix-vector product
 * 
 * Description:
 *   Let mp_alpha be the maximum magnitude of a m x n matrix A, mp_prod-1 be
 *   the maximum magnitude of a n x k matrix C. The function computes the 
 *   maximum interval of positive and negative entries of A.C. That is, the 
 *   function computes mp_b satisfying
 *   (mp_b-1)/2 = n*mp_alpha*(mp_prod-1)
 * 
 * Input:
 *    mp_prod: mpz_t, mp_prod-1 be the maximum magnitude of matrix C
 *   mp_alpha: mpz_t, maximum magnitude of matrix A
 *          n: long, column dimension of A
 *
 * Output:
 *   mp_b: mpz_t, shown above
 *
 */

void
maxInter (const mpz_t mp_prod, const mpz_t mp_alpha, const long n, mpz_t mp_b)
{
  mpz_t mp_temp;

  mpz_init(mp_temp);
  mpz_sub_ui(mp_temp, mp_prod, 1);
  mpz_set(mp_b, mp_alpha);
  mpz_mul_ui(mp_b, mp_b, n);
  mpz_mul(mp_b, mp_b, mp_temp);
  mpz_mul_ui(mp_b, mp_b, 2);
  mpz_add_ui(mp_b, mp_b, 1);
  mpz_clear(mp_temp);
}



/*
 *
 * Calling Sequence:
 *   maxExtInter(mp_alpha, n, mp_b)
 *
 * Summary:
 *   Compute the maximum interval of positive and negative results for 
 *   lifting
 * 
 * Description:
 *   Let mp_alpha be the maximum magnitude of a m x n matrix A, 
 *   The function computes the mp_b satisfying
 *   (mp_b-1)/2 = n*mp_alpha+1
 * 
 * Input:
 *   mp_alpha: mpz_t, maximum magnitude of matrix A
 *          n: long, column dimension of A
 *
 * Output:
 *   mp_b: mpz_t, shown above
 *
 */

void
maxExtInter (const mpz_t mp_alpha, const long n, mpz_t mp_b)
{
  mpz_set_ui(mp_b, 1);
  mpz_addmul_ui(mp_b, mp_alpha, n);
  mpz_mul_ui(mp_b, mp_b, 2);
  mpz_add_ui(mp_b, mp_b, 1);
}



/*
 *
 * Calling Sequence:
 *   bdcoeff <-- repBound(len, basis, cmbasis)
 *
 * Summary:
 *   Compute the mix radix coefficients of a special integer in a RNS basis
 *
 * Description:
 *   Given a RNS basis, suppose the product of elements in the basis be q, 
 *   then this RNS basis is able to represent integers lying in 
 *   [-(q-1)/2, (q-1)/2] and [0, q-1] respectively with symmetric 
 *   representation and positive representation. To transfer the result from
 *   positive representation to symmetric representation, the function 
 *   computes the mix radix coefficients of the boundary value (q-1)/2 in the
 *   positive representation.
 *
 *   Let RNS basis be P. The function computes coefficient array U, such that
 * (q-1)/2 = U[0] + U[1]*P[0] + U[2]*P[0]*P[1] +...+ U[len-1]*P[0]*...*P[len-2]
 *
 * 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
 *
 * Output:
 *   bdcoeff: 1-dim FiniteField array length len, the coefficient array U above
 *
 */

FiniteField *
repBound (const long len, const FiniteField *basis, const FiniteField *cmbasis)
{
  long i, j;
  double dtemp;
  mpz_t mp_bd, mp_prod;
  FiniteField *bdcoeff;
  const FiniteField *q, *qinv;

  q = basis;
  qinv = cmbasis;

  /* set the bound of transformation from positive to negative */
  mpz_init(mp_prod);
  basisProd(len, q, mp_prod);
  mpz_init(mp_bd);
  mpz_sub_ui(mp_bd, mp_prod, 1);
  mpz_divexact_ui(mp_bd, mp_bd, 2);

  /* compute the coeffcients of bound of mix radix and store in bdcoeff */
  bdcoeff = XMALLOC(FiniteField, len);
  bdcoeff[0] = mpz_fdiv_ui(mp_bd, q[0]);
  for (i = 1; i < len; i++)
    {
      dtemp = (double)bdcoeff[i-1];
      for (j = i-2; j >= 0; j--)
	{ 
	  dtemp = fmod(dtemp*q[j], (double)q[i]);
	  dtemp = fmod(dtemp+bdcoeff[j], (double)q[i]);
	}
      bdcoeff[i] = mpz_fdiv_ui(mp_bd, q[i]);
      dtemp = fmod((double)bdcoeff[i]-dtemp, (double)q[i]);
      if (dtemp < 0) { dtemp = q[i]+dtemp; }
      bdcoeff[i] = (FiniteField)fmod(dtemp*qinv[i], (double)q[i]);
    }
  mpz_clear(mp_bd);
  mpz_clear(mp_prod);

  return bdcoeff;
}



/*
 * Calling Sequence:
 *   bd <-- RNSbound(n)
 *
 * Summary:
 *   Compute the upper bound of a RNS basis
 *
 * Description:
 *   Given a m x n mod p matrix A, and a n x k mod p matrix B, the maximum 
 *   magnitude of A.B is n*(p-1)^2. To use BLAS, it is needed that 
 *   n*(p-1)^2 <= 2^53-1 to make the result of product correct.
 *
 *   The function computes an integer bd, such that 
 *      n*(bd-1)^2 <= 2^53-1 and n*((bd+1)-1)^2 > 2^53-1
 * 
 * Input:
 *   n: long, column dimension of matrix A
 *
 * Output:
 *   bd: FiniteField, shown above
 *
 */

FiniteField
RNSbound (const long n)
{
  FiniteField bd;
  mpz_t mp_n, mp_d, mp_q;

  mpz_init(mp_n);
  mpz_init_set_ui(mp_d, n);
  mpz_init(mp_q);
  mpz_ui_pow_ui(mp_n, 2, 53);
  mpz_sub_ui(mp_n, mp_n, 1);
  mpz_fdiv_q(mp_q, mp_n, mp_d);
  mpz_sqrt(mp_q, mp_q);
  bd = mpz_get_ui(mp_q)+1;
  mpz_clear(mp_n);
  mpz_clear(mp_d);
  mpz_clear(mp_q);

  return bd;
}