zzxfactoring.cpp

数值算法库for Windows

void RecInitTab(TBL_T ***lookup_tab, long i, const vec_ulong& ratio, 
             long r, long k, unsigned long thresh1, long **shamt_tab,
             unsigned long sum, long card, long j)
{
   if (j >= i || card >= k-1) {
      if (card > 1) {
         long shamt = shamt_tab[i][card];
         unsigned long index1 = ((-sum) >> shamt);
         lookup_tab[i][card][index1 >> TBL_SHAMT] |= (1UL << (index1 & TBL_MSK));
         unsigned long index2 = ((-sum+thresh1) >> shamt);
         if (index1 != index2)
            lookup_tab[i][card][index2 >> TBL_SHAMT] |= (1UL << (index2 & TBL_MSK));

      }

      return;
   }


   RecInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab, sum, card, j+1);
   RecInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab, 
              sum+ratio[r-1-j], card+1, j+1);
}


static
void DoInitTab(TBL_T ***lookup_tab, long i, const vec_ulong& ratio, 
               long r, long k, unsigned long thresh1, long **shamt_tab)
{
   RecInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab, 0, 0, 0);
}

#else

// iterative version


static
void DoInitTab(TBL_T ***lookup_tab, long i, const vec_ulong& ratio,
               long r, long k, unsigned long thresh1, long **shamt_tab)
{
   vec_long sum_vec, card_vec, location_vec;
   sum_vec.SetLength(i+1);
   card_vec.SetLength(i+1);
   location_vec.SetLength(i+1);

   long j = 0;
   sum_vec[0] = 0;
   card_vec[0] = 0;

   unsigned long sum;
   long  card, location;

   location = 0;

   while (j >= 0) {
      sum = sum_vec[j];
      card = card_vec[j];

      switch (location) {

      case 0:

         if (j >= i || card >= k-1) {
            if (card > 1) {
               long shamt = shamt_tab[i][card];
               unsigned long index1 = ((-sum) >> shamt);
               lookup_tab[i][card][index1 >> TBL_SHAMT] |= (1UL << (index1 & TBL_MSK));
               unsigned long index2 = ((-sum+thresh1) >> shamt);
               if (index1 != index2)
                  lookup_tab[i][card][index2 >> TBL_SHAMT] |= (1UL << (index2 & TBL_MSK));
      
            }
      
            location = location_vec[j];
            j--;
            continue;
         }


         sum_vec[j+1] = sum;
         card_vec[j+1] = card;
         location_vec[j+1] = 1;
         j++;
         location = 0;
         continue;

      case 1:

         sum_vec[j+1] = sum+ratio[r-1-j];
         card_vec[j+1] = card+1;
         location_vec[j+1] = 2;
         j++;
         location = 0;
         continue;  

      case 2:

         location = location_vec[j];
         j--;
         continue;
      }
   }
}
         
#endif
   
   

static
void InitTab(TBL_T ***lookup_tab, const vec_ulong& ratio, long r, long k,
             unsigned long thresh1, long **shamt_tab, long pruning)
{
   long i, j, t;

   if (pruning) {
      for (i = 2; i <= pruning; i++) {
         long len = min(k-1, i);
         for (j = 2; j <= len; j++) {
            long ub = (((1L << (NTL_BITS_PER_LONG-shamt_tab[i][j])) 
                      + TBL_MSK) >> TBL_SHAMT); 
            for (t = 0; t < ub; t++)
               lookup_tab[i][j][t] = 0;
         }
   
         DoInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab);
      }
   }
}


static
void RatioInit1(vec_ulong& ratio, const vec_ZZ_pX& W, const ZZ_p& lc,
                long pruning, TBL_T ***lookup_tab, 
                vec_vec_ulong& pair_ratio, long k, unsigned long thresh1, 
                long **shamt_tab)
{
   long r = W.length();
   long i, j;

   ZZ_p a;

   ZZ p;
   p = ZZ_p::modulus();

   ZZ aa;

   for (i = 0; i < r; i++) {
      long m = deg(W[i]);
      mul(a, W[i].rep[m-1], lc);
      LeftShift(aa, rep(a), NTL_BITS_PER_LONG);
      div(aa, aa, p);
      ratio[i] = to_ulong(aa);
   }

   InitTab(lookup_tab, ratio, r, k, thresh1, shamt_tab, pruning);

   for (i = 0; i < r; i++)
      for (j = 0; j < i; j++) {
         mul(a, W[i].rep[deg(W[i])-1], W[j].rep[deg(W[j])-1]);
         mul(a, a, lc);
         LeftShift(aa, rep(a), NTL_BITS_PER_LONG);
         div(aa, aa, p);
         pair_ratio[i][j] = to_ulong(aa);
      }

   for (i = 0; i < r; i++) {
      long m = deg(W[i]);
      if (m >= 2) {
         mul(a, W[i].rep[m-2], lc);
         LeftShift(aa, rep(a), NTL_BITS_PER_LONG);
         div(aa, aa, p);
         pair_ratio[i][i] = to_ulong(aa);
      }
      else
         pair_ratio[i][i] = 0;
   }
}

static 
long SecondOrderTest(const vec_long& I_vec, const vec_vec_ulong& pair_ratio_vec,
                     vec_ulong& sum_stack_vec, long& SumLen)
{
   long k = I_vec.length();
   const long *I = I_vec.elts();
   unsigned long *sum_stack = sum_stack_vec.elts();

   unsigned long sum, thresh1;

   if (SumLen == 0) {
      unsigned long epsilon = (1UL << (NTL_BITS_PER_LONG-ZZX_OVERLIFT));
      unsigned long delta = (unsigned long) ((k*(k+1)) >> 1);
      unsigned long thresh = epsilon + delta;
      thresh1 = (epsilon << 1) + delta;

      sum = thresh;
      sum_stack[k] = thresh1;
   }
   else {
      sum = sum_stack[SumLen-1];
      thresh1 = sum_stack[k];
   }

   long i, j;

   for (i = SumLen; i < k; i++) {
      const unsigned long *p = pair_ratio_vec[I[i]].elts();
      for (j = 0; j <= i; j++) {
         sum += p[I[j]];
      }

      sum_stack[i] = sum;
   }

   SumLen = k-1;

   return (sum <= thresh1);
}


static
ZZ choose_fn(long r, long k)
{
   ZZ a, b;

   a = 1; 
   b = 1;

   long i;
   for (i = 0; i < k; i++) {
      a *= r-i;
      b *= k-i;
   }

   return a/b;
}

static
void PrintInfo(const char *s, const ZZ& a, const ZZ& b)
{
   cerr << s << a << " / " << b << " = ";
   
   double x = to_double(a)/to_double(b);

   if (x == 0) 
      cerr << "0"; 
   else {
      int n;
      double f;

      f = frexp(x, &n);
      cerr << f << "*2^" << n;
   }

   cerr << "\n";
}

static
void RemoveFactors1(vec_long& W, const vec_long& I, long r)
{
   long k = I.length();
   long i, j;

   i = 0;
   for (j = 0; j < r; j++) {
      if (i < k && j == I[i])
         i++;
      else
         swap(W[j-i], W[j]); 
   }
}

static
void RemoveFactors1(vec_vec_long& W, const vec_long& I, long r)
{
   long k = I.length();
   long i, j;

   i = 0;
   for (j = 0; j < r; j++) {
      if (i < k && j == I[i])
         i++;
      else
         swap(W[j-i], W[j]); 
   }

   for (i = 0; i < r-k; i++)
      RemoveFactors1(W[i], I, r);
}


// should this swap go in tools.h?
// Maybe not...I don't want to pollute the interface too much more.

inline void swap(unsigned long& a, unsigned long& b)  
   { unsigned long t;  t = a; a = b; b = t; }

static
void RemoveFactors1(vec_ulong& W, const vec_long& I, long r)
{
   long k = I.length();
   long i, j;

   i = 0;
   for (j = 0; j < r; j++) {
      if (i < k && j == I[i])
         i++;
      else
         swap(W[j-i], W[j]); 
   }
}

static
void RemoveFactors1(vec_vec_ulong& W, const vec_long& I, long r)
{
   long k = I.length();
   long i, j;

   i = 0;
   for (j = 0; j < r; j++) {
      if (i < k && j == I[i])
         i++;
      else
         swap(W[j-i], W[j]); 
   }

   for (i = 0; i < r-k; i++)
      RemoveFactors1(W[i], I, r);
}


static
void RemoveFactors1(vec_ZZ_p& W, const vec_long& I, long r)
{
   long k = I.length();
   long i, j;

   i = 0;
   for (j = 0; j < r; j++) {
      if (i < k && j == I[i])
         i++;
      else
         swap(W[j-i], W[j]);
   }
}

static
void SumCoeffs(ZZ& sum, const ZZX& a)
{
   ZZ res;
   res = 0;
   long i;
   long n = a.rep.length();
   for (i = 0; i < n; i++)
      res += a.rep[i];

   sum = res;
}

static
void SumCoeffs(ZZ_p& sum, const ZZ_pX& a)
{
   ZZ_p res;
   res = 0;
   long i;
   long n = a.rep.length();
   for (i = 0; i < n; i++)
      res += a.rep[i];

   sum = res;
}


static
long ConstTermTest(const vec_ZZ_p& W, 
                  const vec_long& I,
                  const ZZ& ct,
                  const ZZ_p& lc,
                  vec_ZZ_p& prod,
                  long& ProdLen) 
{
   long k = I.length();
   ZZ_p t;
   ZZ t1, t2;
   long i;

   if (ProdLen == 0) {
      mul(prod[0], lc, W[I[0]]);
      ProdLen++;
   }

   for (i = ProdLen; i < k; i++)
      mul(prod[i], prod[i-1], W[I[i]]);

   ProdLen = k-1;

   // should make this a routine in ZZ_p
   t1 = rep(prod[k-1]);
   RightShift(t2, ZZ_p::modulus(), 1);
   if (t1 > t2)
      sub(t1, t1, ZZ_p::modulus());

   return divide(ct, t1);
}


long ZZXFac_MaxPrune = 10;



static
long pruning_bnd(long r, long k)
{
   double x = 0; 

   long i;
   for (i = 0; i < k; i++) {
      x += log(double(r-i)/double(k-i));
   }

   return long((x/log(2.0)) * 0.75);
}

static
long shamt_tab_init(long pos, long card, long pruning, long thresh1_len)
{
   double x = 1;
   long i;

   for (i = 0; i < card; i++) {
      x *= double(pos-i)/double(card-i);
   }

   x *= pruning;  // this can be adjusted to control the density
   if (pos <= 6) x *= 2;  // a little boost that costs very little
      

   long t = long(ceil(log(x)/log(2.0)));

   t = max(t, TBL_SHAMT); 

   t = min(t, NTL_BITS_PER_LONG-thresh1_len);


   return NTL_BITS_PER_LONG-t;
}

// The following routine should only be called for k > 1,
// and is only worth calling for k > 2.


static
void CardinalitySearch1(vec_ZZX& factors, ZZX& f, 
                       vec_ZZ_pX& W, 
                       LocalInfoT& LocalInfo, 
                       long k,
                       long bnd,
                       long verbose)
{
   double start_time, end_time;

   if (verbose) {
      start_time = GetTime();
      cerr << "\n************ ";
      cerr << "start cardinality " << k << "\n";
   }

   if (k <= 1) Error("internal error: call CardinalitySearch");

   // This test is needed to ensure correcntes of "n-2" test
   if (NumBits(k) > NTL_BITS_PER_LONG/2-2)
      Error("Cardinality Search: k too large...");

   vec_ZZ pdeg;
   CalcPossibleDegrees(pdeg, W, k);
   ZZ pd;

   bit_and(pd, pdeg[0], LocalInfo.PossibleDegrees);
   if (pd == 0) {
      if (verbose) cerr << "skipping\n";
      return;
   }

   vec_long I, D;
   I.SetLength(k);
   D.SetLength(k);

   long r = W.length();

   long initial_r = r;

   vec_ulong ratio, ratio_sum;
   ratio.SetLength(r);
   ratio_sum.SetLength(k);

   unsigned long epsilon = (1UL << (NTL_BITS_PER_LONG-ZZX_OVERLIFT));
   unsigned long delta = (unsigned long) k;
   unsigned long thresh = epsilon + delta;
   unsigned long thresh1 = (epsilon << 1) + delta;

   long thresh1_len = NumBits(long(thresh1)); 

   long pruning;

   pruning = min(r/2, ZZXFac_MaxPrune);
   pruning = min(pruning, pruning_bnd(r, k));
   pruning = min(pruning, NTL_BITS_PER_LONG-EXTRA_BITS-thresh1_len);

   if (pruning <= 4) pruning = 0;

   long init_pruning = pruning;

   TBL_T ***lookup_tab = 0;

   long **shamt_tab = 0;

   if (pruning) {
      typedef long *long_p;

      long i, j;

      shamt_tab = NTL_NEW_OP long_p[pruning+1];
      if (!shamt_tab) Error("out of mem");
      shamt_tab[0] = shamt_tab[1] = 0;

      for (i = 2; i <= pruning; i++) {
         long len = min(k-1, i);
         shamt_tab[i] = NTL_NEW_OP long[len+1];
         if (!shamt_tab[i]) Error("out of mem");
         shamt_tab[i][0] = shamt_tab[i][1] = 0;

         for (j = 2; j <= len; j++)
            shamt_tab[i][j] = shamt_tab_init(i, j, pruning, thresh1_len);
      }

      typedef  TBL_T *TBL_T_p;
      typedef  TBL_T **TBL_T_pp;

      lookup_tab = NTL_NEW_OP TBL_T_pp[pruning+1];
      if (!lookup_tab) Error("out of mem");

      lookup_tab[0] = lookup_tab[1] = 0;

      for (i = 2; i <= pruning; i++) {
         long len = min(k-1, i);
         lookup_tab[i] = NTL_NEW_OP TBL_T_p[len+1];
         if (!lookup_tab[i]) Error("out of mem");

         lookup_tab[i][0] = lookup_tab[i][1] = 0;

         for (j = 2; j <= len; j++) {
            lookup_tab[i][j] = NTL_NEW_OP TBL_T[((1L << (NTL_BITS_PER_LONG-shamt_tab[i][j]))+TBL_MSK) >> TBL_SHAMT];
            if (!lookup_tab[i][j]) Error("out of mem");
         }
      }
   }

   if (verbose) {
      cerr << "pruning = " << pruning << "\n";
   }

   vec_ZZ_p prod;
   prod.SetLength(k);
   long ProdLen;

   vec_ZZ_p prod1;
   prod1.SetLength(k);
   long ProdLen1;

   vec_ulong sum_stack;
   sum_stack.SetLength(k+1);
   long SumLen;

   vec_long upd;

   long i, state;

   long cnt = 0;

   ZZ ct;
   mul(ct, ConstTerm(f), LeadCoeff(f));

   ZZ_p lc;
   conv(lc, LeadCoeff(f));

   vec_vec_ulong pair_ratio;
   pair_ratio.SetLength(r);
   for (i = 0; i < r; i++)
      pair_ratio[i].SetLength(r);

   RatioInit1(ratio, W, lc, pruning, lookup_tab, pair_ratio, k, thresh1, shamt_tab);

   ZZ c1;
   SumCoeffs(c1, f);
   mul(c1, c1, LeadCoeff(f));

   vec_ZZ_p sum_coeffs;
   sum_coeffs.SetLength(r);
   for (i = 0; i < r; i++)
      SumCoeffs(sum_coeffs[i], W[i]);

   vec_long degv;
   degv.SetLength(r);

   for (i = 0; i < r; i++)
      degv[i] = deg(W[i]);

   ZZ_pX gg;