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

/*
 * Calling Sequence:
 *   Dmod(p, A, n, m, lda)
 * 
 * Summary:
 *   Perform mod p operation inplace over a double matrix/submatrix
 *
 * Description:
 *   Given a n x m matrix or submatrix A and an integer p, this function uses
 *   fmod function to do operation A mod p inplace. Although the type of input
 *   modulus p in this function is Double, it is a cast of an integer.
 *
 * Input:
 *     p: Double, module
 *     A: 1-dim Double array length n*m, representation of a n x m double 
 *        matrix/submatrix
 *     n: long, row dimension of A
 *     m: long, column dimension of A
 *   lda: long, number of entries between two continuous rows of A, help track
 *        the pointer if A is a submatrix
 *
 */

void
Dmod (const Double p, Double *A, const long n, const long m, const long lda)
{
  long i, j;
  Double temp;

  for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
      *(A+i*lda+j) = (temp = fmod(*(A+i*lda+j), p)) >= 0 ? temp : p+temp;
}

  
/*
 * Calling Sequence:
 *   DCopy(n, m, A, lda, B, ldb)
 * 
 * Summary:
 *   Copy Double matrix/submatrix A to another matrix/submatrix B
 *
 * Input:
 *     n: long, row dimension of matrix/submatrix A and B
 *     m: long, column dimension of matrix/submatrix A and B
 *     A: 1-dim Double array length n*m, representation of a n x m Double 
 *        matrix/submatrix
 *   lda: long, stride of two continuous rows of A
 *     B: 1-dim Double array length n*m, representation of a n x m Double 
 *        matrix/submatrix
 *   ldb: long, stride of two continuous rows of B
 *
 */

void 
DCopy (const long n, const long m, const Double* A, const long lda, \
       Double* B, const long ldb)
{
  long i, j;

  for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
      *(B+i*ldb+j) = *(A+i*lda+j);
}



/*
 * Calling Sequence:
 *   randomDb(n, m, bd)
 * 
 * Summary:
 *   Generate a random Double dense matrix
 * 
 * Description:
 *   Given the bound bd, this function generates a random Double dense matrix
 *   such that any entry e satisfies 0 <= e <= 2^bd-1
 *
 * Input:
 *    n: long, row dimension of the generated matrix
 *    m: long, column dimension of the generated matrix
 *   bd: long, magnitude bound of entries
 *
 * Return:
 *   M: 1-dim Double array length n*m, representation of a n x m random matrix
 *
 */

Double *
randomDb (const long n, const long m, const long bd)
{
  long i, j;
  mpz_t mp_rand;
  gmp_randstate_t state;
  unsigned long seed;
  FILE* devrandom;
  Double* M;

#if HAVE_TIME_H
  time_t t1;
#endif

  M = XMALLOC(Double, n*m);
  mpz_init(mp_rand);
  gmp_randinit_default(state);
  mpz_set_ui(mp_rand, 5);

  /* generate random seed using /dev/random */
  if ((devrandom = fopen("/dev/urandom", "r")) != NULL)
    {
      for (i = 0; i < n*m; i++)
	{
	  fread(&seed, sizeof(seed), 1, devrandom);
	  gmp_randseed_ui(state, seed);
	  mpz_urandomb(mp_rand, state, bd);
	  M[i] = mpz_get_d(mp_rand); 
	}
      fclose(devrandom);
    }
  else
    {
      for (i = 0; i < n*m; i++)
	{
#if HAVE_TIME_H
	  time(&t1);
	  seed = (unsigned long)t1;
#endif
	  gmp_randseed_ui(state, seed);
	  mpz_urandomb(mp_rand, state, bd);
	  M[i] = mpz_get_d(mp_rand); 
	}
    }
  mpz_clear(mp_rand);
  gmp_randclear(state);
  return M;
}



/*
 * Calling Sequence:
 *   RandomPrime(lb, hb)
 *
 * Summary:
 *   Generate a random prime p satisfying 2^lb <= p <= 2^hb-1
 *
 * Input:
 *   lb: FiniteField, lower bound of the prime
 *   hb: FiniteField, upper bound of the prime
 *
 * Return:
 *   p: FiniteField, randomly generated prime
 *
 * Note: 
 *   lb and hb can not exceed the bit-length of unsigned long.
 *
 */

FiniteField
RandPrime (const FiniteField lb, const FiniteField hb)
{
  mpz_t mp_rand, mp_temp, mp_lb, mp_hb;
  gmp_randstate_t state;
  FiniteField p;
  unsigned long seed;
  FILE* devrandom;

#if HAVE_TIME_H
  time_t t1;
#endif

  { mpz_init(mp_rand); mpz_init(mp_temp); } 
  { mpz_init(mp_lb);  mpz_init(mp_hb); }
  mpz_ui_pow_ui(mp_lb, 2, lb);
  mpz_ui_pow_ui(mp_hb, 2, hb);
  mpz_sub(mp_temp, mp_hb, mp_lb);
  gmp_randinit_default(state);

  /* generate random seed using /dev/urandom */
  if ((devrandom = fopen("/dev/urandom", "r")) != NULL)
    {
      fread(&seed, sizeof(seed), 1, devrandom);
      fclose(devrandom);
    }
  else
    {
#if HAVE_TIME_H
      time(&t1);
      seed = (unsigned long)t1;
#endif
    }
  gmp_randseed_ui(state, seed);
  mpz_urandomm(mp_rand, state, mp_temp); /* 0 <= mp_rand <= 2^hb-2^lb-1 */
  mpz_add(mp_rand, mp_rand, mp_lb);      /* 2^lb <= mp_rand <= 2^hb-1 */
  while (mpz_probab_prime_p(mp_rand, 10) == 0) 
    mpz_sub_ui(mp_rand, mp_rand, 1); 
  p = mpz_get_ui(mp_rand);

  { mpz_clear(mp_rand); mpz_clear(mp_temp); }
  { mpz_clear(mp_lb); mpz_clear(mp_hb); }
  gmp_randclear(state);
  return p;
}



/*
 * Calling Sequence:
 *   max <-- maxMagnLong(A, n, m, lda)
 *
 * Summary:
 *   Compute maximum magnitude of a n x m signed long matrix/submatrix A
 *
 * Input:
 *     A: 1-dim long array length n*m, representation of a n x m 
 *        matrix/submatrix
 *     n: long, row dimension of A
 *     m: long, column dimension of A
 *   lda: long, number of entries of two continuous rows of A (normally m)
 *
 * Return: 
 *   max: long, maximum magnitude of A
 *
 */

long 
maxMagnLong (const long *A, const long n, const long m, const long lda)
{
  long i, j, temp, max=0;

  for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
      if ((temp = labs(A[i*lda+j])) > max) { max = temp; }
  return max;
}



/*
 *
 * Calling Sequence:
 *   maxMagnMP(mp_A, n, m, lda, mp_max)
 *
 * Summary:
 *   Compute maximum magnitude of a n x m mpz_t matrix/submatrix mp_A
 *
 * Input:
 *   mp_A: 1-dim mpz_t array length n*m, representation of a n x m 
 *         matrix/submatrix
 *      n: long, row dimension of mp_A
 *      m: long, column dimension of mp_A
 *    lda: long, number of entries of two continuous rows of mp_A (normally m)
 *
 * Output: 
 *   mp_max: mpz_t, maximum magnitude of mp_A
 *
 */

void
maxMagnMP (mpz_t *mp_A, const long n, const long m, const long lda, \
	   mpz_t mp_max)
{
  long i, j;

  mpz_set_ui(mp_max, 0);
  for (i = 0; i < n; i++) 
    for (j = 0; j < m; j++)
      if (mpz_cmpabs(mp_A[i*lda+j], mp_max) > 0) 
	mpz_abs(mp_max, mp_A[i*lda+j]);
  return;
}


/*
 * Calling Sequence:
 *   scalCpMP(n, m, lda, ldm, mp_a, mp_A, mp_M)
 *
 * Summary:
 *   Copy a submatrix with scale from one matrix to the other one
 *
 * Description:
 *   Copy with scalar n x m submatrix from A to M. The submatrix starts
 *   from position pointed by mp_A with stride lda. Then the submatrix is 
 *   copied to position pointed by mp_M with stride ldm, meanwhile, the 
 *   entries of the submatrix are scaled by factor mp_a.
 * 
 * Input: 
 *      n: long, row dimension of the submatrix
 *      m: long, column dimension of the submatrix
 *    lda: long, stride of two continuous rows of the submatrix in A
 *    ldm: long, stride of two continuous rows of the submatrix in M
 *   mp_a: mpz_t, scalar factor
 *   mp_A: mpz_t pointer, start point of the submatrix in A
 *   mp_M: mpz_t pointer, start point of the submatrix in M
 *
 */

void 
scalCpMP (const long n, const long m, const long lda, const long ldm, \
	  const mpz_t mp_a, mpz_t *mp_A, mpz_t *mp_M)
{
  long i, j;

  if (mpz_cmp_ui(mp_a, 1) != 0)
    for (i = 0; i < n; i++)
      for (j = 0; j < m; j++)
	{
	  mpz_set(mp_M[i*ldm+j], mp_A[i*lda+j]);
	  mpz_mul(mp_M[i*ldm+j], mp_M[i*ldm+j], mp_a);
	}
  else
    for (i = 0; i < n; i++)
      for (j = 0; j < m; j++)
	mpz_set(mp_M[i*ldm+j], mp_A[i*lda+j]);
}