kChinese2.cc

Go to the documentation of this file.

#include "misc/auxiliary.h"
 
#include "misc/intvec.h"
#include "polys/monomials/p_polys.h"
#include "polys/matpol.h"
#include "polys/simpleideals.h"
 
#include <gmp.h>
#include <stdlib.h>
#include <stdio.h>
 
// Garner's Algorithm.
// See Algorithm 14.71, Handbook of Cryptography.
 
//    x - result    v residuals    m - primes   t-size of vectors
//    prod: product of m_i
// prepare:
static void prepareCRT(mpz_ptr C[], mpz_ptr *m, mpz_t prod, int t)
{
  mpz_t u;
  int i, j;
 
  mpz_init(u);
  mpz_init_set(prod,m[0]);
  for (i=1; i<t; i++)
  {
    mpz_mul(prod,prod,m[i]);
    mpz_init(C[i]);
    mpz_set_ui(C[i], 1);
    for (j=0; j<i; j++)
    {
      mpz_invert(u, m[j], m[i]);
      mpz_mul(C[i], C[i], u);
      mpz_mod(C[i], C[i], m[i]);
    }
  }
  mpz_clear(u);
}
 
// the proper CRT
//    x - result    v residuals    m - primes   t-size of vectors
static void CRT(mpz_t x, mpz_ptr *v, mpz_ptr* C, mpz_ptr *m, int t)
{
  mpz_t u;
  int i, j;
 
  mpz_init(u);
  mpz_set(u, v[0]);
  mpz_set(x, u);
  for (i=1; i<t; i++)
  {
    mpz_sub(u, v[i], x);
    mpz_mul(u, u, C[i]);
    mpz_mod(u, u, m[i]);
    for (j=0; j<i; j++)
    {
      mpz_mul(u, u, m[j]);
    }
    mpz_add(x, x, u);
  }
  mpz_clear(u);
}
 
/*2
* xx,q: arrays of length 0..rl-1
* xx[i]: SB mod q[i]
* assume: char=0
* assume: q[i]!=0
* destroys xx
*/
poly p_ChineseRemainder(poly *xx, mpz_ptr *x,mpz_ptr *q, int rl, mpz_ptr *C, mpz_t prod2, mpz_t prod,const ring R)
{
  poly r,h,hh;
  int j;
  poly  res_p=NULL;
  loop
  {
    /* search the lead term */
    r=NULL;
    for(j=rl-1;j>=0;j--)
    {
      h=xx[j];
      if ((h!=NULL)
      &&((r==NULL)||(p_LmCmp(r,h,R)==-1)))
        r=h;
    }
    /* nothing found -> return */
    if (r==NULL) break;
    /* create the monomial in h */
    h=p_Head(r,R);
    /* collect the coeffs in x[..]*/
    for(j=rl-1;j>=0;j--)
    {
      hh=xx[j];
      if ((hh!=NULL) && (p_LmCmp(h,hh,R)==0))
      {
        mpz_set_si(x[j],n_Int(pGetCoeff(hh),R->cf));
        hh=p_LmFreeAndNext(hh,R);
        xx[j]=hh;
      }
      else
        mpz_set_si(x[j],0);
    }
    mpz_t nn; mpz_init(nn);
    CRT(nn,x,C,q,rl);
    if (mpz_cmp(nn,prod2)>0) mpz_sub(nn,nn,prod);
    number n=n_InitMPZ(nn,R->cf);
    if (n_IsZero(n,R->cf)) p_Delete(&h,R);
    else
    {
      //Print("new mon:");pWrite(h);
      p_SetCoeff(h,n,R);
      pNext(h)=res_p;
      res_p=h; // building res_p in reverse order!
    }
  }
  res_p=pReverse(res_p);
  p_Test(res_p, R);
  return res_p;
}
 
ideal id_ChineseRemainder_0(ideal *xx, number *q, int rl, const ring r)
{
  int cnt=0;int rw=0; int cl=0;
  int i,j;
  // find max. size of xx[.]:
  for(j=rl-1;j>=0;j--)
  {
    i=IDELEMS(xx[j])*xx[j]->nrows;
    if (i>cnt) cnt=i;
    if (xx[j]->nrows >rw) rw=xx[j]->nrows; // for lifting matrices
    if (xx[j]->ncols >cl) cl=xx[j]->ncols; // for lifting matrices
  }
  if (rw*cl !=cnt)
  {
    WerrorS("format mismatch in CRT");
    return NULL;
  }
  ideal result=idInit(cnt,xx[0]->rank);
  result->nrows=rw; // for lifting matrices
  result->ncols=cl; // for lifting matrices
  mpz_ptr x[rl]; // work space
  mpz_ptr qq[rl];// moduli
  mpz_ptr C[rl];// Cache
  for(j=rl-1;j>=0;j--)
  {
    x[j]=(mpz_ptr)omAlloc(sizeof(mpz_t));
    mpz_init(x[j]);
    qq[j]=(mpz_ptr)omAlloc(sizeof(mpz_t));
    n_MPZ(qq[j],q[j],r->cf);
    C[j]=(mpz_ptr)omAlloc(sizeof(mpz_t));
  }
  poly *p=(poly *)omAlloc(rl*sizeof(poly));
  //    x - result    v residuals    m - primes   t-size of vectors
  // prepare:
  mpz_t prod,prod2;
  prepareCRT(C, qq, prod, rl);
  mpz_init_set(prod2,prod);
  mpz_tdiv_ui(prod2,2);
  for(i=cnt-1;i>=0;i--)
  {
    for(j=rl-1;j>=0;j--)
    {
      if(i>=IDELEMS(xx[j])*xx[j]->nrows) // out of range of this ideal
        p[j]=NULL;
      else
        p[j]=xx[j]->m[i];
    }
    result->m[i]=p_ChineseRemainder(p,x,qq,rl,C,prod2,prod,r);
    for(j=rl-1;j>=0;j--)
    {
      if(i<IDELEMS(xx[j])*xx[j]->nrows) xx[j]->m[i]=p[j];
    }
  }
  mpz_clear(prod2);
  mpz_clear(prod);
  omFreeSize(p,rl*sizeof(poly));
  for(j=rl-1;j>=0;j--)
  {
    mpz_clear(qq[j]);omFree(qq[j]);
    mpz_clear(x[j]);omFree(x[j]);
  }
  for(j=rl-1;j>0;j--)
  {
    mpz_clear(C[j]); omFree(C[j]);
  }
  for(i=rl-1;i>=0;i--) id_Delete(&(xx[i]),r);
  omFreeSize(xx,rl*sizeof(ideal));
  return result;
}