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

/*
 * Calling Sequence:
 *   sumliftCoeff(mp_basisprod, k, C, mp_sum)
 *
 * Summary:
 *   Recursively sum up the p-adic lifting coefficients
 * 
 * Description:
 *   Given p-adic lifing coefficients C after k lifting steps, the function
 *   recursively compute 
 *     A^(-1).B mod p^k = 
 *     C[0]+C[1]*mp_basisprod+ ... +C[k-1]*mp_basisprod^(k-1)
 *
 * Input: 
 *   mp_basisprod: mpz_t, product of lifting basis
 *              k: long, dimension of array C, number of lifting steps
 *              C: 1-dim mpz_t matrix length k, storing lifting coefficients
 *                 computed each lifting step
 *
 * Output:
 *   mp_sum: mpz_t, the sum of lifting coefficients
 *
 */

void
sumliftCoeff (const mpz_t mp_basisprod, const long k, mpz_t* C, mpz_t mp_sum)
{
  long i, t, splflag;
  long idx[1];
  mpz_t mp_right, *mp_pow, *mp_left;

  idx[0] = 0;

  /* precomputation mp_basisprod^log(k) */
  t = find2exp(k);
  mp_pow = XMALLOC(mpz_t, t+1);
  mpz_init_set(mp_pow[0], mp_basisprod);
  for (i = 1; i < t+1; i++) 
    { 
      mpz_init(mp_pow[i]);
      mpz_pow_ui(mp_pow[i], mp_pow[i-1], 2);
    }
  mpz_init(mp_right);
  if (t == 0) 
    { 
      mpz_set(mp_sum, C[0]); 
      for (i = 0; i < t+1; i++) { mpz_clear(mp_pow[i]); } { XFREE(mp_pow); }
      mpz_clear(mp_right);
      return; 
    }
  if ((1 << t) == k) { splflag = 1; }
  else { splflag = 0; }
  mp_left = XMALLOC(mpz_t, t);
  for (i = 0; i < t; i++) { mpz_init(mp_left[i]); }
  sumCoeff_rec(0, k, C, mp_pow, splflag, 0, idx, mp_left, mp_right);
  mpz_set(mp_sum, mp_left[0]);

  for (i = 0; i < t+1; i++) { mpz_clear(mp_pow[i]); } { XFREE(mp_pow); }
  for (i = 0; i < t; i++) { mpz_clear(mp_left[i]); } { XFREE(mp_left); }
  mpz_clear(mp_right);
  return;
}


/*
 * Calling Sequence:
 *   sumCoeff_rec(s, len, C, mp_pow, splflag, 0, idx, mp_left, mp_right)
 *
 * Summary:
 *   compute C[s]+C[s+1]*p+C[s+2]*p^2+...+C[s+len-1]*p^(len-1) recursively
 *
 * Input:
 *      start: long, recursion start index in array C, usually 0
 *        len: long, number of entries in C to use
 *          C: 1-dim mpz_t array, storing values for computation
 *     mp_pow: 1-dim mpz_t array length t+1, storing different powers of p,
 *             mp_pow[i] = p^(2^i), (i = 0..t), t = floor(log(len)/log(2))
 *    splflag: 1/0, input flag, 
 *           - if len=2^i and i >= 0, then splflag = 1
 *           - else splflag = 0
 *    savflag: 1/0, input flag, initially use 0
 *           - if savflag = 0, then current subproblem is in the left subtree 
 *             of recursion tree
 *           - else, current subproblem is in the right subtree of recursion 
 *             tree
 *        idx: pointer to a long integer, the index of available entries in
 *             array mp_left, initially *idx = 0
 *    mp_left: 1-dim mpz_t array length t, storing the intermidiate result of 
 *             left subtree
 *   mp_right: mpz_t, storing the result of right subtree
 *
 * Output: 
 *   mp_left[0] will store the computation result
 *
 * Idea: the result of subproblem on the right subtree of recursion tree is 
 *       overwritten into mp_right
 *       the result of subproblem on the left subtree of recursion tree is 
 *       stored into mp_left[*idx], and make the *idx always be availabe 
 *       index in mp_left, e.g., the space of mp_left used in right subtree can
 *       can be reused by its siblings and parent
 *
 */

void 
sumCoeff_rec (long start, long len, mpz_t *C, mpz_t *mp_pow, long splflag,\
	      long savflag, long *idx, mpz_t *mp_left, mpz_t mp_right)
{
  long lpow, llen, tidx;

  /* basecase */
  if (len == 1) 
    {
      if (savflag == 0) 
	{
	  mpz_set(mp_left[*idx], C[start]); 
	  ++*idx;
	}
      else { mpz_set(mp_right, C[start]); }
      return;
    }
  /* if len is power of 2, then directly divide len by 2 to
     obtain two subproblems*/
  if (splflag == 1)
    {
      lpow = find2exp(len)-1;
      llen = 1 << lpow;
      sumCoeff_rec(start, llen, C, mp_pow, 1, 0, idx, mp_left, mp_right);
      tidx = *idx-1;
      sumCoeff_rec(start+llen, len-llen, C, mp_pow, 1, 1, idx, mp_left, \
		   mp_right);
    }
  /* otherwise, find maximal lpow, 2^lpow < len, to make two subproblems */
  else
    {
      lpow = find2exp(len);
      llen = 1 << lpow;
      sumCoeff_rec(start, llen, C, mp_pow, 1, 0, idx, mp_left, mp_right);
      tidx = *idx-1;

      /* if len-2^lpow is power of 2, then not compute the second subproblem */
      if (llen == len)
	{
	  mpz_set(mp_right, mp_left[tidx]);
	  return;
	}
      sumCoeff_rec(start+llen, len-llen, C, mp_pow, 0, 1, idx, mp_left, \
		   mp_right);
    }

  /* combine subproblems by position of two subproblems in the recursion tree*/
  /* if subproblem lies in the left subtree, store result in a new mpz_t int */
  if (savflag == 0)
    mpz_addmul(mp_left[tidx], mp_right, mp_pow[lpow]);

  /* if in the right subtree, store result in the pre-allocated mp_right */
  else
    {
      mpz_mul(mp_right, mp_right, mp_pow[lpow]);
      mpz_add(mp_right, mp_left[tidx], mp_right);
    }
  *idx = tidx+1;
  return;
}


/*
 * Calling Sequence:
 *   i <-- find2exp(len)
 * 
 * Summary:
 *   Compute floor(log[2](len))
 *
 * Description:
 *   Compute i >= 0 such that 2^i <= len and 2^(i+1) > len
 *
 */

long
find2exp (const long len)  
{
  long i=0, l=1;
  while ((l <<= 1) <= len) { ++i; }
  return i;
}


/*
 * Calling Sequence:
 *   1/0 <-- iratrecon(mp_u, mp_m, mp_nb, mp_db, mp_N, mp_D)
 *
 * Summary:
 *   Perform rational number reconstruction
 *
 * Description:
 *   Given image mp_u, modulus mp_m, the function reconstruct numerator mp_N 
 *   and denominator mp_D, such that 
 *   mp_N/mp_D = mp_u (mod mp_m), and
 *   abs(mp_N) <= mp_nb, 0 < mp_D <= mp_db
 *
 * Input: