source: git/dyn_modules/syzextra/syzextra.cc @ bfeed6

spielwiese
Last change on this file since bfeed6 was bfeed6, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Additional initializations for buckets (for unusual entry points) TODO: remove these 'if'-s
  • Property mode set to 100644
File size: 47.6 KB
Line 
1// -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2/*****************************************************************************\
3 * Computer Algebra System SINGULAR   
4\*****************************************************************************/
5/** @file syzextra.cc
6 *
7 * Here we implement the Computation of Syzygies
8 *
9 * ABSTRACT: Computation of Syzygies due to Schreyer
10 *
11 * @author Oleksandr Motsak
12 *
13 **/
14/*****************************************************************************/
15
16// include header file
17#include <kernel/mod2.h>
18
19#include "syzextra.h"
20
21#include "DebugPrint.h"
22
23#include <omalloc/omalloc.h>
24
25#include <misc/intvec.h>
26#include <misc/options.h>
27
28#include <coeffs/coeffs.h>
29
30#include <polys/monomials/p_polys.h>
31#include <polys/monomials/ring.h>
32#include <polys/simpleideals.h>
33
34#include <polys/kbuckets.h> // for kBucket*
35#include <polys/sbuckets.h> // for sBucket*
36//#include <polys/nc/summator.h> // for CPolynomialSummator
37#include <polys/operations/p_Mult_q.h> // for MIN_LENGTH_BUCKET
38
39#include <kernel/GBEngine/kstd1.h>
40#include <kernel/polys.h>
41#include <kernel/GBEngine/syz.h>
42#include <kernel/ideals.h>
43
44#include <kernel/oswrapper/timer.h>
45
46
47#include <Singular/tok.h>
48#include <Singular/ipid.h>
49#include <Singular/lists.h>
50#include <Singular/attrib.h>
51
52#include <Singular/ipid.h> 
53#include <Singular/ipshell.h> // For iiAddCproc
54
55#include <stdio.h>
56#include <stdlib.h>
57#include <string.h>
58
59// USING_NAMESPACE_SINGULARXX;
60USING_NAMESPACE( SINGULARXXNAME :: DEBUG )
61
62
63BEGIN_NAMESPACE_SINGULARXX     BEGIN_NAMESPACE(SYZEXTRA)
64
65
66BEGIN_NAMESPACE(SORT_c_ds)
67
68
69#ifdef _GNU_SOURCE
70static int cmp_c_ds(const void *p1, const void *p2, void *R)
71{
72#else
73static int cmp_c_ds(const void *p1, const void *p2)
74{
75  void *R = currRing; 
76#endif
77
78  const int YES = 1;
79  const int NO = -1;
80
81  const ring r =  (const ring) R; // TODO/NOTE: the structure is known: C, lp!!!
82
83  assume( r == currRing );
84
85  const poly a = *(const poly*)p1;
86  const poly b = *(const poly*)p2;
87
88  assume( a != NULL );
89  assume( b != NULL );
90
91  assume( p_LmTest(a, r) );
92  assume( p_LmTest(b, r) );
93
94
95  const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r);
96
97  // TODO: test this!!!!!!!!!!!!!!!!
98
99  //return -( compare (c, qsorts) )
100
101#ifndef NDEBUG
102  const int __DEBUG__ = 0;
103  if( __DEBUG__ )
104  {
105    PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2);
106    PrintS("b: "); dPrint(b, r, r, 2);
107    PrintLn();
108  }
109#endif
110
111
112  if( iCompDiff > 0 )
113    return YES;
114
115  if( iCompDiff < 0 )
116    return  NO;
117
118  assume( iCompDiff == 0 );
119
120  const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r);
121
122  if( iDegDiff > 0 )
123    return YES;
124
125  if( iDegDiff < 0 )
126    return  NO;
127
128  assume( iDegDiff == 0 );
129
130#ifndef NDEBUG
131  if( __DEBUG__ )
132  {
133    PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn();
134  }
135#endif
136
137  for (int v = rVar(r); v > 0; v--)
138  {
139    assume( v > 0 );
140    assume( v <= rVar(r) );
141
142    const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r);
143
144    if( d > 0 )
145      return YES;
146
147    if( d < 0 )
148      return NO;
149
150    assume( d == 0 );
151  }
152
153  return 0; 
154}
155
156END_NAMESPACE
157/* namespace SORT_c_ds */
158
159/// return a new term: leading coeff * leading monomial of p
160/// with 0 leading component!
161poly leadmonom(const poly p, const ring r, const bool bSetZeroComp)
162{
163  poly m = NULL;
164
165  if( p != NULL )
166  {
167    assume( p != NULL );
168    assume( p_LmTest(p, r) );
169
170    m = p_LmInit(p, r);
171    p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);
172
173    if( bSetZeroComp )
174      p_SetComp(m, 0, r);
175    p_Setm(m, r);
176
177   
178    assume( m != NULL );
179    assume( pNext(m) == NULL );
180    assume( p_LmTest(m, r) );
181   
182    if( bSetZeroComp )
183      assume( p_GetComp(m, r) == 0 );
184  }
185
186  return m;
187}
188
189
190   
191poly p_Tail(const poly p, const ring r)
192{
193  if( p == NULL)
194    return NULL;
195  else
196    return p_Copy( pNext(p), r );
197}
198
199
200ideal id_Tail(const ideal id, const ring r)
201{
202  if( id == NULL)
203    return NULL;
204
205  const ideal newid = idInit(IDELEMS(id),id->rank);
206 
207  for (int i=IDELEMS(id) - 1; i >= 0; i--)
208    newid->m[i] = p_Tail( id->m[i], r );
209
210  newid->rank = id_RankFreeModule(newid, currRing);
211
212  return newid; 
213}
214
215
216
217void Sort_c_ds(const ideal id, const ring r)
218{
219  const int sizeNew = IDELEMS(id);
220
221#ifdef _GNU_SOURCE
222#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r)
223#else
224#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp)
225#endif
226 
227  if( sizeNew >= 2 )
228    qsort_my(id->m, sizeNew, sizeof(poly), r, FROM_NAMESPACE(SORT_c_ds, cmp_c_ds));
229
230#undef qsort_my
231
232  id->rank = id_RankFreeModule(id, r);
233}
234
235/// Clean up all the accumulated data
236void SchreyerSyzygyComputation::CleanUp()
237{
238  extern void id_Delete (ideal*, const ring);
239   
240  id_Delete(const_cast<ideal*>(&m_idTails), m_rBaseRing); // TODO!!!
241
242  if( m_sum_bucket != NULL )
243  {
244    sBucketDestroy(&m_sum_bucket);
245    m_sum_bucket = NULL;
246  }
247
248  if( m_spoly_bucket != NULL )
249  {
250    kBucketDestroy(&m_spoly_bucket);
251    m_spoly_bucket = NULL;
252  } 
253}
254  /*
255  for( TTailTerms::const_iterator it = m_idTailTerms.begin(); it != m_idTailTerms.end(); it++ )
256  {
257    const TTail& v = *it;
258    for(TTail::const_iterator vit = v.begin(); vit != v.end(); vit++ )
259      delete const_cast<CTailTerm*>(*vit);
260  }
261  */
262
263
264
265int CReducerFinder::PreProcessTerm(const poly t, CReducerFinder& syzChecker) const
266{
267  assume( t != NULL );
268
269  if( __DEBUG__ & __TAILREDSYZ__ )
270    assume( !IsDivisible(t) ); // each input term should NOT be in <L>
271
272  const ring r = m_rBaseRing;
273
274
275  if( __TAILREDSYZ__ )
276    if( p_LmIsConstant(t, r) ) // most basic case of baing coprime with L, whatever that is...
277      return 1; // TODO: prove this...?
278
279  //   return false; // appears to be fine
280
281  const long comp = p_GetComp(t, r);
282
283  CReducersHash::const_iterator itr = m_hash.find(comp);
284
285  if ( itr == m_hash.end() ) 
286    return 2; // no such leading component!!!
287
288  const bool bIdealCase = (comp == 0);   
289  const bool bSyzCheck = syzChecker.IsNonempty(); // need to check even in ideal case????? proof?  "&& !bIdealCase"
290
291  if( __TAILREDSYZ__ && (bIdealCase || bSyzCheck) )
292  {
293    const TReducers& v = itr->second;
294    const int N = rVar(r);
295    // TODO: extract exps of t beforehand?!
296    bool coprime = true;
297    for(TReducers::const_iterator vit = v.begin(); (vit != v.end()) && coprime; ++vit )
298    {
299      assume( m_L->m[(*vit)->m_label] == (*vit)->m_lt );
300
301      const poly p = (*vit)->m_lt;
302
303      assume( p_GetComp(p, r) == comp ); 
304
305      // TODO: check if coprime with Leads... if __TAILREDSYZ__ !
306      for( int var = N; var > 0; --var )
307        if( (p_GetExp(p, var, r) != 0) && (p_GetExp(t, var, r) != 0) )
308        {       
309#ifndef NDEBUG
310          if( __DEBUG__ | 0)
311          {         
312            PrintS("CReducerFinder::PreProcessTerm, 't' is NOT co-prime with the following leading term: \n");
313            dPrint(p, r, r, 1);
314          }
315#endif     
316          coprime = false; // t not coprime with p!
317          break;
318        }
319
320      if( bSyzCheck && coprime )
321      {
322        poly ss = p_LmInit(t, r);
323        p_SetCoeff0(ss, n_Init(1, r), r); // for delete & printout only!...
324        p_SetComp(ss, (*vit)->m_label + 1, r); // coeff?
325        p_Setm(ss, r);
326
327        coprime = ( syzChecker.IsDivisible(ss) );
328
329#ifndef NDEBUG
330        if( __DEBUG__ && !coprime)
331        {         
332          PrintS("CReducerFinder::PreProcessTerm, 't' is co-prime with p but may lead to NOT divisible syz.term: \n");
333          dPrint(ss, r, r, 1);
334        }
335#endif
336         
337        p_LmDelete(&ss, r); // deletes coeff as well???
338      }
339
340    }
341
342#ifndef NDEBUG
343    if( __DEBUG__ && coprime )
344      PrintS("CReducerFinder::PreProcessTerm, the following 't' is 'co-prime' with all of leading terms! \n");
345#endif
346     
347    return coprime? 3: 0; // t was coprime with all of leading terms!!!
348
349  }
350  //   return true; // delete the term
351
352  return 0;
353}
354 
355   
356void SchreyerSyzygyComputation::SetUpTailTerms()
357{
358  const ideal idTails = m_idTails; 
359  assume( idTails != NULL );
360  assume( idTails->m != NULL );
361  const ring r = m_rBaseRing;
362
363#ifndef NDEBUG
364  if( __DEBUG__ | 0)
365  {
366    PrintS("SchreyerSyzygyComputation::SetUpTailTerms(): Tails: \n");
367    dPrint(idTails, r, r, 0);
368  }
369
370  unsigned long pp[4] = {0,0,0,0}; // count preprocessed terms...
371#endif
372
373  for( int p = IDELEMS(idTails) - 1; p >= 0; --p )
374    for( poly* tt = &(idTails->m[p]); (*tt) != NULL;  )
375    {   
376      const poly t = *tt;
377      const int k = m_div.PreProcessTerm(t, m_checker); // 0..3
378      assume( 0 <= k && k <= 3 );
379
380#ifndef NDEBUG
381      pp[k]++;
382#endif
383       
384      if( k )
385      {
386#ifndef NDEBUG
387        if( __DEBUG__)
388        {         
389          Print("SchreyerSyzygyComputation::SetUpTailTerms(): PP (%d) the following TT: \n", k);
390          dPrint(t, r, r, 1);
391        }
392#endif
393
394        (*tt) = p_LmDeleteAndNext(t, r); // delete the lead and next...
395      }   
396      else 
397        tt = &pNext(t); // go next?
398
399    }
400
401#ifndef NDEBUG
402  if( TEST_OPT_PROT | 1)
403    Print("      **!!**      SchreyerSyzygyComputation::SetUpTailTerms()::PreProcessing(): X: {c: %lu, C: %lu, P: %lu} + %lu\n", pp[1], pp[2], pp[3], pp[0]);
404#endif
405   
406#ifndef NDEBUG
407  if( __DEBUG__ | 0)
408  {
409    PrintS("SchreyerSyzygyComputation::SetUpTailTerms(): Preprocessed Tails: \n");
410    dPrint(idTails, r, r, 0);
411  }
412#endif
413   
414}
415/* 
416  m_idTailTerms.resize( IDELEMS(idTails) );
417 
418  for( unsigned int p = IDELEMS(idTails) - 1; p >= 0; p -- )
419  {
420    TTail& v = m_idTailTerms[p];
421    poly t = idTails->m[p];
422    v.resize( pLength(t) );
423
424    unsigned int pp = 0;
425   
426    while( t != NULL )
427    {
428      assume( t != NULL );
429      // TODO: compute L:t!
430//      ideal reducers;     
431//      CReducerFinder m_reducers
432     
433      CTailTerm* d = v[pp] = new CTailTerm();
434     
435      ++pp; pIter(t);
436    }
437  }
438*/
439
440
441
442ideal SchreyerSyzygyComputation::Compute1LeadingSyzygyTerms()
443{
444  const ideal& id = m_idLeads;
445  const ring& r = m_rBaseRing;
446//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
447
448//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
449//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
450//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
451//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
452//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
453
454  assume(!__LEAD2SYZ__);
455
456  // 1. set of components S?
457  // 2. for each component c from S: set of indices of leading terms
458  // with this component?
459  // 3. short exp. vectors for each leading term?
460
461  const int size = IDELEMS(id);
462
463  if( size < 2 )
464  {
465    const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal...       
466    return newid;
467  }
468
469  // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!??
470
471  // components should come in groups: count elements in each group
472  // && estimate the real size!!!
473
474
475  // use just a vector instead???
476  const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
477
478  int k = 0;
479
480  for (int j = 0; j < size; j++)
481  {
482    const poly p = id->m[j];
483    assume( p != NULL );
484    const int  c = p_GetComp(p, r);
485
486    for (int i = j - 1; i >= 0; i--)
487    {
488      const poly pp = id->m[i];
489      assume( pp != NULL );
490      const int  cc = p_GetComp(pp, r);
491
492      if( c != cc )
493        continue;
494
495      const poly m = p_Init(r); // p_New???
496
497      // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)!
498      for (int v = rVar(r); v > 0; v--)
499      {
500        assume( v > 0 );
501        assume( v <= rVar(r) );
502
503        const short e1 = p_GetExp(p , v, r);
504        const short e2 = p_GetExp(pp, v, r);
505
506        if( e1 >= e2 )
507          p_SetExp(m, v, 0, r);
508        else
509          p_SetExp(m, v, e2 - e1, r);
510
511      }
512
513      assume( (j > i) && (i >= 0) );
514
515      p_SetComp(m, j + 1, r);
516      pNext(m) = NULL;
517      p_SetCoeff0(m, n_Init(1, r->cf), r); // for later...
518
519      p_Setm(m, r); // should not do anything!!!
520
521      newid->m[k++] = m;
522    }
523  }
524
525//   if( __DEBUG__ && FALSE )
526//   {
527//     PrintS("ComputeLeadingSyzygyTerms::Temp0: \n");
528//     dPrint(newid, r, r, 1);
529//   }
530
531  // the rest of newid is assumed to be zeroes...
532
533  // simplify(newid, 2 + 32)??
534  // sort(newid, "C,ds")[1]???
535  id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
536
537//   if( __DEBUG__ && FALSE )
538//   {
539//     PrintS("ComputeLeadingSyzygyTerms::Temp1: \n");
540//     dPrint(newid, r, r, 1);
541//   }
542
543  idSkipZeroes(newid); // #define SIMPL_NULL 2
544
545//   if( __DEBUG__ )
546//   {
547//     PrintS("ComputeLeadingSyzygyTerms::Output: \n");
548//     dPrint(newid, r, r, 1);
549//   }
550
551  Sort_c_ds(newid, r);
552
553  return newid;
554}
555
556ideal SchreyerSyzygyComputation::Compute2LeadingSyzygyTerms()
557{
558  const ideal& id = m_idLeads;
559  const ring& r = m_rBaseRing;
560//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
561 
562//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
563//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
564//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
565//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
566//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
567 
568
569  // 1. set of components S?
570  // 2. for each component c from S: set of indices of leading terms
571  // with this component?
572  // 3. short exp. vectors for each leading term?
573
574  const int size = IDELEMS(id);
575
576  if( size < 2 )
577  {
578    const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...   
579    return newid;
580  }
581
582
583  // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
584 
585  // components should come in groups: count elements in each group
586  // && estimate the real size!!!
587
588
589  // use just a vector instead???
590  ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
591
592  int k = 0;
593
594  for (int j = 0; j < size; j++)
595  {
596    const poly p = id->m[j];
597    assume( p != NULL );
598    const int  c = p_GetComp(p, r);
599
600    for (int i = j - 1; i >= 0; i--)
601    {
602      const poly pp = id->m[i];
603      assume( pp != NULL );
604      const int  cc = p_GetComp(pp, r);
605
606      if( c != cc )
607        continue;
608
609        // allocate memory & zero it out!
610      const poly m = p_Init(r); const poly mm = p_Init(r);
611
612
613        // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
614        // TODO: optimize: knowing the ring structure: (C/lp)!
615
616      for (int v = rVar(r); v > 0; v--)
617      {
618        assume( v > 0 );
619        assume( v <= rVar(r) );
620
621        const short e1 = p_GetExp(p , v, r);
622        const short e2 = p_GetExp(pp, v, r);
623
624        if( e1 >= e2 )
625          p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
626        else
627          p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
628
629      }
630
631      assume( (j > i) && (i >= 0) );
632
633      p_SetComp(m, j + 1, r);
634      p_SetComp(mm, i + 1, r);
635
636      const number& lc1 = p_GetCoeff(p , r);
637      const number& lc2 = p_GetCoeff(pp, r);
638
639      number g = n_Lcm( lc1, lc2, r );
640
641      p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
642      p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);
643
644      n_Delete(&g, r);
645
646      p_Setm(m, r); // should not do anything!!!
647      p_Setm(mm, r); // should not do anything!!!
648
649      pNext(m) = mm; //        pNext(mm) = NULL;
650
651      newid->m[k++] = m;
652    }
653  }
654
655//   if( __DEBUG__ && FALSE )
656//   {
657//     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
658//     dPrint(newid, r, r, 1);
659//   }
660
661  if( !__TAILREDSYZ__ )
662  {
663      // simplify(newid, 2 + 32)??
664      // sort(newid, "C,ds")[1]???
665    id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
666
667//     if( __DEBUG__ && FALSE )
668//     {
669//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
670//       dPrint(newid, r, r, 1);
671//     }
672  }
673  else
674  {
675      //      option(redSB); option(redTail);
676      //      TEST_OPT_REDSB
677      //      TEST_OPT_REDTAIL
678    assume( r == currRing );
679 
680    BITSET _save_test; SI_SAVE_OPT1(_save_test);
681    SI_RESTORE_OPT1(Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB) | _save_test);
682
683    intvec* w=new intvec(IDELEMS(newid));
684    ideal tmp = kStd(newid, currQuotient, isHomog, &w);
685    delete w;
686
687    SI_RESTORE_OPT1(_save_test)
688
689    id_Delete(&newid, r);
690    newid = tmp;
691
692//     if( __DEBUG__ && FALSE )
693//     {
694//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
695//       dPrint(newid, r, r, 1);
696//     }
697
698  }
699
700  idSkipZeroes(newid);
701
702  Sort_c_ds(newid, r);
703 
704  return newid;
705}
706
707poly SchreyerSyzygyComputation::TraverseNF(const poly a, const poly a2) const
708{
709  const ideal& L = m_idLeads;
710  const ideal& T = m_idTails;
711
712  const ring& R = m_rBaseRing;
713
714  const int r = p_GetComp(a, R) - 1; 
715
716  assume( r >= 0 && r < IDELEMS(T) );
717  assume( r >= 0 && r < IDELEMS(L) );
718
719  poly aa = leadmonom(a, R); assume( aa != NULL); // :(
720
721 
722  poly t = TraverseTail(aa, r);
723
724  if( a2 != NULL )
725  {
726    assume( __LEAD2SYZ__ );
727
728    const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
729
730    assume( r2 >= 0 && r2 < IDELEMS(T) );
731
732    t = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, r2), R), R);
733
734    p_Delete(&aa2, R);
735  } else
736    t = p_Add_q(t, ReduceTerm(aa, L->m[r], a), R);
737
738  p_Delete(&aa, R);
739
740  return t;
741}
742
743
744void SchreyerSyzygyComputation::ComputeSyzygy()
745{
746  assume( m_idLeads != NULL );
747  assume( m_idTails != NULL );
748
749  const ideal& L = m_idLeads;
750  const ideal& T = m_idTails;
751
752  ideal& TT = m_syzTails;
753  const ring& R = m_rBaseRing;
754
755  if( m_sum_bucket == NULL )
756    m_sum_bucket = sBucketCreate(R);
757
758  if( m_spoly_bucket == NULL )
759    m_spoly_bucket = kBucketCreate(R);
760
761
762  assume( IDELEMS(L) == IDELEMS(T) );
763#ifndef NDEBUG
764  int t, r; 
765#endif
766   
767  if( m_syzLeads == NULL )
768  {   
769#ifndef NDEBUG
770    if( TEST_OPT_PROT | 1)
771    {
772      t = getTimer(); r = getRTimer();
773      Print("%5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::ComputeLeadingSyzygyTerms: t: %d, r: %d\n", r, t, r);
774    }
775#endif     
776    ComputeLeadingSyzygyTerms( __LEAD2SYZ__ && !__IGNORETAILS__ ); // 2 terms OR 1 term!
777#ifndef NDEBUG
778    if( TEST_OPT_PROT | 1)
779    {
780      t = getTimer() - t; r = getRTimer() - r;
781      Print("%5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::ComputeLeadingSyzygyTerms: dt: %d, dr: %d\n", getRTimer(), t, r);
782    }
783#endif
784     
785  }
786
787  assume( m_syzLeads != NULL );
788  ideal& LL = m_syzLeads;
789  const int size = IDELEMS(LL);
790
791  TT = idInit(size, 0);
792
793  if( size == 1 && LL->m[0] == NULL )
794    return;
795
796  // use hybrid method?
797  const bool method = (__HYBRIDNF__ == 1) || (__HYBRIDNF__ == 2 && __SYZNUMBER__ < 3);
798
799  if(  !__IGNORETAILS__)
800  {
801    if( T != NULL )
802    {
803#ifndef NDEBUG
804      if( TEST_OPT_PROT | 1 )
805      {
806        t = getTimer(); r = getRTimer();
807        Print("%5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SetUpTailTerms(): t: %d, r: %d\n", r, t, r);
808      }
809#endif
810       
811      SetUpTailTerms();
812#ifndef NDEBUG
813      if( TEST_OPT_PROT | 1)
814      {
815        t = getTimer() - t; r = getRTimer() - r;
816        Print("%5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SetUpTailTerms(): dt: %d, dr: %d\n", getRTimer(), t, r);
817      }
818#endif
819    }     
820  }
821
822#ifndef NDEBUG
823  if( TEST_OPT_PROT | 1)
824  {
825    t = getTimer(); r = getRTimer();
826    Print("%5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SyzygyLift: t: %d, r: %d\n", r, t, r);
827  }
828#endif
829   
830  for( int k = size - 1; k >= 0; k-- )
831  {
832    const poly a = LL->m[k]; assume( a != NULL );
833
834    poly a2 = pNext(a);   
835
836    // Splitting 2-terms Leading syzygy module
837    if( a2 != NULL )
838      pNext(a) = NULL;
839
840    if( __IGNORETAILS__ )
841    {
842      TT->m[k] = NULL;
843
844      assume( a2 != NULL );
845
846      if( a2 != NULL )
847        p_Delete(&a2, R);
848
849      continue;
850    }
851
852    //    TT->m[k] = a2;
853
854    if( method )
855      TT->m[k] = SchreyerSyzygyNF(a, a2);
856    else
857      TT->m[k] = TraverseNF(a, a2);
858  }
859
860#ifndef NDEBUG
861  if( TEST_OPT_PROT | 1)
862  {
863    t = getTimer() - t; r = getRTimer() - r;
864    Print("%5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SyzygyLift: dt: %d, dr: %d\n", getRTimer(), t, r);
865  }
866#endif   
867
868  TT->rank = id_RankFreeModule(TT, R); 
869}
870
871void SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms(bool bComputeSecondTerms)
872{
873//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
874
875//  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
876//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
877
878  assume( m_syzLeads == NULL );
879
880  if( bComputeSecondTerms )
881  {
882    assume( __LEAD2SYZ__ );
883//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _Compute2LeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
884    m_syzLeads = Compute2LeadingSyzygyTerms();
885  }
886  else
887  {
888    assume( !__LEAD2SYZ__ );
889   
890    m_syzLeads = Compute1LeadingSyzygyTerms();
891  }
892//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _ComputeLeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
893 
894  // NOTE: set m_LS if tails are to be reduced!
895  assume( m_syzLeads!= NULL );
896
897  if (__TAILREDSYZ__ && !__IGNORETAILS__ && (IDELEMS(m_syzLeads) > 0) && !((IDELEMS(m_syzLeads) == 1) && (m_syzLeads->m[0] == NULL)))
898  {
899    m_LS = m_syzLeads;
900    m_checker.Initialize(m_syzLeads);
901#ifndef NDEBUG   
902    if( __DEBUG__ )
903    {
904      const ring& r = m_rBaseRing;
905      PrintS("SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms: \n");
906      PrintS("m_syzLeads: \n");
907      dPrint(m_syzLeads, r, r, 1);
908      PrintS("m_checker.Initialize(m_syzLeads) => \n");
909      m_checker.DebugPrint();
910    }
911#endif
912    assume( m_checker.IsNonempty() ); // TODO: this always fails... BUG????
913  }
914}
915
916#define NOPRODUCT 1
917
918poly SchreyerSyzygyComputation::SchreyerSyzygyNF(const poly syz_lead, poly syz_2) const
919{
920 
921  assume( !__IGNORETAILS__ );
922
923  const ideal& L = m_idLeads;
924  const ideal& T = m_idTails;
925  const ring& r = m_rBaseRing;
926
927  assume( syz_lead != NULL );
928
929  if( syz_2 == NULL )
930  {
931    const int rr = p_GetComp(syz_lead, r) - 1; 
932
933    assume( rr >= 0 && rr < IDELEMS(T) );
934    assume( rr >= 0 && rr < IDELEMS(L) );
935
936
937#if NOPRODUCT
938    syz_2 = m_div.FindReducer(syz_lead, L->m[rr], syz_lead, m_checker);
939#else   
940    poly aa = leadmonom(syz_lead, r); assume( aa != NULL); // :(
941    aa = p_Mult_mm(aa, L->m[rr], r);
942
943    syz_2 = m_div.FindReducer(aa, syz_lead, m_checker);
944
945    p_Delete(&aa, r);
946#endif
947   
948    assume( syz_2 != NULL ); // by construction of S-Polynomial   
949  }
950
951
952 
953  assume( syz_2 != NULL );
954
955  assume( L != NULL );
956  assume( T != NULL );
957
958  assume( IDELEMS(L) == IDELEMS(T) );
959
960  int  c = p_GetComp(syz_lead, r) - 1;
961
962  assume( c >= 0 && c < IDELEMS(T) );
963
964  if( m_sum_bucket == NULL )
965    m_sum_bucket = sBucketCreate(r);
966
967  if( m_spoly_bucket == NULL )
968    m_spoly_bucket = kBucketCreate(r);
969
970  sBucket_pt& tail   = m_sum_bucket; assume( tail != NULL );
971  kBucket_pt& bucket = m_spoly_bucket; assume( bucket != NULL ); kbTest(bucket);
972
973
974//  kBucketInit(bucket, NULL, 0); // not needed!?
975     
976  poly p = leadmonom(syz_lead, r); // :( 
977//  poly spoly = pp_Mult_qq(p, T->m[c], r); 
978  kBucket_Plus_mm_Mult_pp(bucket, p, T->m[c], 0); // TODO: store pLength(T->m[c]) separately!?
979  p_Delete(&p, r);
980 
981  kbTest(bucket); 
982
983  c = p_GetComp(syz_2, r) - 1;
984  assume( c >= 0 && c < IDELEMS(T) );
985
986  p = leadmonom(syz_2, r); // :(
987//  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
988  kBucket_Plus_mm_Mult_pp(bucket, p, T->m[c], 0); // pLength(T->m[c])?!
989  kbTest(bucket);
990  p_Delete(&p, r);
991
992  // TODO: use bucket!?
993//  const bool bUsePolynomial = TEST_OPT_NOT_BUCKETS; //  || (pLength(spoly) < MIN_LENGTH_BUCKET);
994//  CPolynomialSummator tail(r, bUsePolynomial);                       
995  sBucket_Add_p(tail, syz_2, 1); // tail.AddAndDelete(syz_2, 1); 
996 
997  kbTest(bucket);
998  for( poly spoly = kBucketExtractLm(bucket); spoly != NULL; p_LmDelete(&spoly, r), spoly = kBucketExtractLm(bucket))
999  {   
1000    kbTest(bucket);
1001    poly t = m_div.FindReducer(spoly, NULL, m_checker);   
1002
1003    if( t != NULL )
1004    {
1005      p = leadmonom(t, r); // :(
1006      c = p_GetComp(t, r) - 1;
1007
1008      assume( c >= 0 && c < IDELEMS(T) );
1009
1010      kBucket_Plus_mm_Mult_pp(bucket, p, T->m[c], 0); // pLength(T->m[c])?
1011//      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
1012
1013      p_Delete(&p, r);
1014
1015      sBucket_Add_p(tail, t, 1); // tail.AddAndDelete(t, 1);
1016    } // otherwise discard that leading term altogether!
1017    kbTest(bucket);
1018  }
1019
1020  // now bucket must be empty!
1021  kbTest(bucket);
1022  assume( kBucketClear(bucket) == NULL );
1023
1024  poly result; int len;
1025  sBucketClearAdd(tail, &result, &len); // TODO: use Merge with sBucket???
1026  assume( pLength(result) == len );
1027     
1028  return result;
1029}
1030
1031poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, const int tail) const
1032{
1033  // TODO: store (multiplier, tail) -.-^-.-^-.--> !
1034  assume(m_idTails !=  NULL && m_idTails->m != NULL);
1035  assume( tail >= 0 && tail < IDELEMS(m_idTails) );
1036
1037  const poly t = m_idTails->m[tail]; // !!!
1038
1039  if(t != NULL)
1040    return TraverseTail(multiplier, t);
1041
1042  return NULL;
1043}
1044
1045
1046poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
1047{
1048  assume( !__IGNORETAILS__ );
1049
1050  const ideal& L = m_idLeads;
1051  const ideal& T = m_idTails;
1052  const ring& r = m_rBaseRing;
1053
1054  assume( multiplier != NULL );
1055
1056  assume( L != NULL );
1057  assume( T != NULL );
1058
1059
1060  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
1061  {
1062//    const bool bUsePolynomial = TEST_OPT_NOT_BUCKETS; //  || (pLength(tail) < MIN_LENGTH_BUCKET);
1063    if( m_sum_bucket == NULL )
1064      m_sum_bucket = sBucketCreate(r);
1065
1066    sBucket_pt& sum   = m_sum_bucket; assume( sum != NULL );
1067    poly s; int len;
1068   
1069//    CPolynomialSummator sum(r, bUsePolynomial);
1070//    poly s = NULL;
1071    for(poly p = tail; p != NULL; p = pNext(p))   // iterate over the tail
1072    {
1073      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1074      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1075      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1076      const poly rt = ReduceTerm(multiplier, p, NULL); // TODO: also return/store length?
1077      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1078      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1079      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1080     
1081      const int lp = pLength(rt);
1082      if( rt != NULL && lp != 0 )
1083        sBucket_Add_p(sum, rt, lp); 
1084    }
1085
1086    sBucketClearAdd(sum, &s, &len);
1087    assume( pLength(s) == len );
1088    return s;
1089  }
1090
1091  return NULL;
1092
1093}
1094
1095
1096
1097
1098poly SchreyerSyzygyComputation::ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck) const
1099{
1100  assume( !__IGNORETAILS__ );
1101
1102  const ideal& L = m_idLeads;
1103  const ideal& T = m_idTails;
1104  const ring& r = m_rBaseRing;
1105
1106  assume( multiplier != NULL );
1107  assume( term4reduction != NULL );
1108
1109
1110  assume( L != NULL );
1111  assume( T != NULL );
1112
1113  // simple implementation with FindReducer:
1114  poly s = NULL;
1115
1116  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
1117  {
1118#if NOPRODUCT
1119    s = m_div.FindReducer(multiplier, term4reduction, syztermCheck, m_checker);
1120#else   
1121    // NOTE: only LT(term4reduction) should be used in the following:
1122    poly product = pp_Mult_mm(multiplier, term4reduction, r);
1123    s = m_div.FindReducer(product, syztermCheck, m_checker);
1124    p_Delete(&product, r);
1125#endif
1126  }
1127
1128  if( s == NULL ) // No Reducer?
1129    return s;
1130
1131  poly b = leadmonom(s, r);
1132
1133  const int c = p_GetComp(s, r) - 1;
1134  assume( c >= 0 && c < IDELEMS(T) );
1135
1136  const poly t = TraverseTail(b, c); // T->m[c];
1137
1138  if( t != NULL )
1139    s = p_Add_q(s, t, r); 
1140
1141  return s;
1142}
1143
1144
1145
1146
1147
1148BEGIN_NAMESPACE_NONAME
1149   
1150static inline int atGetInt(idhdl rootRingHdl, const char* attribute, long def)
1151{
1152  return ((int)(long)(atGet(rootRingHdl, attribute, INT_CMD, (void*)def)));
1153}
1154
1155END_NAMESPACE   
1156
1157SchreyerSyzygyComputationFlags::SchreyerSyzygyComputationFlags(idhdl rootRingHdl):
1158#ifndef NDEBUG
1159     __DEBUG__( atGetInt(rootRingHdl,"DEBUG", 0) ),
1160#else
1161    __DEBUG__( atGetInt(rootRingHdl,"DEBUG", 0) ),
1162#endif
1163//    __SYZCHECK__( (BOOLEAN)atGetInt(rootRingHdl, "SYZCHECK", __DEBUG__) ),
1164    __LEAD2SYZ__( atGetInt(rootRingHdl, "LEAD2SYZ", 1) ), 
1165    __TAILREDSYZ__( atGetInt(rootRingHdl, "TAILREDSYZ", 1) ), 
1166    __HYBRIDNF__( atGetInt(rootRingHdl, "HYBRIDNF", 2) ),
1167    __IGNORETAILS__( atGetInt(rootRingHdl, "IGNORETAILS", 0) ),
1168    __SYZNUMBER__( atGetInt(rootRingHdl, "SYZNUMBER", 0) ),
1169    m_rBaseRing( rootRingHdl->data.uring )
1170{   
1171#ifndef NDEBUG
1172  if( __DEBUG__ )
1173  {
1174    PrintS("SchreyerSyzygyComputationFlags: \n");
1175    Print("        DEBUG: \t%d\n", __DEBUG__);
1176//    Print("   SYZCHECK  : \t%d\n", __SYZCHECK__);
1177    Print("     LEAD2SYZ: \t%d\n", __LEAD2SYZ__);
1178    Print("   TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
1179    Print("  IGNORETAILS: \t%d\n", __IGNORETAILS__);
1180  }
1181#endif
1182   
1183  // TODO: just current setting!
1184  assume( rootRingHdl == currRingHdl );
1185  assume( rootRingHdl->typ == RING_CMD );
1186  assume( m_rBaseRing == currRing );
1187  // move the global ring here inside???
1188}
1189
1190   
1191
1192CLeadingTerm::CLeadingTerm(unsigned int _label,  const poly _lt, const ring R):
1193    m_sev( p_GetShortExpVector(_lt, R) ),  m_label( _label ),  m_lt( _lt )
1194{ }
1195
1196
1197CReducerFinder::~CReducerFinder()
1198{
1199  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++ )
1200  {
1201    const TReducers& v = it->second;
1202    for(TReducers::const_iterator vit = v.begin(); vit != v.end(); vit++ )
1203      delete const_cast<CLeadingTerm*>(*vit);
1204  }
1205}
1206
1207
1208void CReducerFinder::Initialize(const ideal L)
1209{
1210  assume( m_L == NULL || m_L == L );
1211  if( m_L == NULL )
1212    m_L = L;
1213
1214  assume( m_L == L );
1215 
1216  if( L != NULL )
1217  {
1218    const ring& R = m_rBaseRing;
1219    assume( R != NULL );
1220   
1221    for( int k = IDELEMS(L) - 1; k >= 0; k-- )
1222    {
1223      const poly a = L->m[k]; // assume( a != NULL );
1224
1225      // NOTE: label is k \in 0 ... |L|-1!!!
1226      if( a != NULL )
1227        m_hash[p_GetComp(a, R)].push_back( new CLeadingTerm(k, a, R) );
1228    }
1229  }
1230}
1231
1232CReducerFinder::CReducerFinder(const ideal L, const SchreyerSyzygyComputationFlags& flags):
1233    SchreyerSyzygyComputationFlags(flags),
1234    m_L(const_cast<ideal>(L)), // for debug anyway
1235    m_hash()
1236{
1237  assume( flags.m_rBaseRing == m_rBaseRing );
1238  if( L != NULL )
1239    Initialize(L);
1240}
1241
1242/// _p_LmDivisibleByNoComp for a | b*c
1243static inline BOOLEAN _p_LmDivisibleByNoComp(const poly a, const poly b, const poly c, const ring r)
1244{
1245  int i=r->VarL_Size - 1;
1246  unsigned long divmask = r->divmask;
1247  unsigned long la, lb;
1248
1249  if (r->VarL_LowIndex >= 0)
1250  {
1251    i += r->VarL_LowIndex;
1252    do
1253    {
1254      la = a->exp[i];
1255      lb = b->exp[i] + c->exp[i];
1256      if ((la > lb) ||
1257          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
1258      {
1259        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
1260        return FALSE;
1261      }
1262      i--;
1263    }
1264    while (i>=r->VarL_LowIndex);
1265  }
1266  else
1267  {
1268    do
1269    {
1270      la = a->exp[r->VarL_Offset[i]];
1271      lb = b->exp[r->VarL_Offset[i]] + c->exp[r->VarL_Offset[i]];
1272      if ((la > lb) ||
1273          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
1274      {
1275        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
1276        return FALSE;
1277      }
1278      i--;
1279    }
1280    while (i>=0);
1281  }
1282#ifdef HAVE_RINGS
1283  assume( !rField_is_Ring(r) ); // not implemented for rings...!
1284#endif
1285  return TRUE;
1286}
1287
1288bool CLeadingTerm::DivisibilityCheck(const poly product, const unsigned long not_sev, const ring r) const
1289{
1290  const poly p = m_lt;
1291
1292  assume( p_GetComp(p, r) == p_GetComp(product, r) );
1293
1294  const int k = m_label;
1295
1296//  assume( m_L->m[k] == p );
1297
1298  const unsigned long p_sev = m_sev;
1299
1300  assume( p_sev == p_GetShortExpVector(p, r) );     
1301
1302  return p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r);
1303
1304}
1305
1306/// as DivisibilityCheck(multiplier * t, ...) for monomial 'm'
1307/// and a module term 't'
1308bool CLeadingTerm::DivisibilityCheck(const poly m, const poly t, const unsigned long not_sev, const ring r) const
1309{
1310  const poly p = m_lt;
1311
1312  assume( p_GetComp(p, r) == p_GetComp(t, r) );
1313// assume( p_GetComp(m, r) == 0 );
1314
1315//  const int k = m_label;
1316//  assume( m_L->m[k] == p );
1317
1318  const unsigned long p_sev = m_sev;
1319  assume( p_sev == p_GetShortExpVector(p, r) );     
1320
1321  if (p_sev & not_sev)
1322    return false;
1323
1324  return _p_LmDivisibleByNoComp(p, m, t, r);
1325
1326//  return p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r);
1327
1328}
1329
1330
1331
1332/// TODO:
1333class CDivisorEnumerator: public SchreyerSyzygyComputationFlags
1334{
1335  private: 
1336    const CReducerFinder& m_reds;
1337    const poly m_product;
1338    const unsigned long m_not_sev;
1339    const long m_comp;
1340
1341    CReducerFinder::CReducersHash::const_iterator m_itr;
1342    CReducerFinder::TReducers::const_iterator m_current, m_finish;
1343
1344    bool m_active;
1345
1346  public:
1347    CDivisorEnumerator(const CReducerFinder& self, const poly product):
1348        SchreyerSyzygyComputationFlags(self),
1349        m_reds(self),
1350        m_product(product),
1351        m_not_sev(~p_GetShortExpVector(product, m_rBaseRing)),
1352        m_comp(p_GetComp(product, m_rBaseRing)),
1353        m_itr(), m_current(), m_finish(),
1354        m_active(false)
1355    {
1356      assume( m_comp >= 0 );
1357      assume( m_reds.m_L != NULL ); 
1358    }
1359
1360    inline bool Reset()
1361    {
1362      m_active = false;
1363     
1364      m_itr = m_reds.m_hash.find(m_comp);
1365
1366      if( m_itr == m_reds.m_hash.end() )
1367        return false;
1368
1369      assume( m_itr->first == m_comp );
1370
1371      m_current = (m_itr->second).begin();
1372      m_finish = (m_itr->second).end();
1373
1374      if (m_current == m_finish)
1375        return false;
1376
1377//      m_active = true;
1378      return true;     
1379    }
1380
1381    const CLeadingTerm& Current() const
1382    {
1383      assume( m_active );
1384      assume( m_current != m_finish );
1385
1386      return *(*m_current);
1387    }
1388 
1389    inline bool MoveNext()
1390    {
1391      assume( m_current != m_finish );
1392
1393      if( m_active )
1394        ++m_current;
1395      else
1396        m_active = true; // for Current()
1397     
1398      // looking for the next good entry
1399      for( ; m_current != m_finish; ++m_current )
1400      {
1401        assume( m_reds.m_L->m[Current().m_label] == Current().m_lt );
1402
1403        if( Current().DivisibilityCheck(m_product, m_not_sev, m_rBaseRing) )
1404        {
1405#ifndef NDEBUG
1406          if( __DEBUG__ )
1407          {
1408            Print("CDivisorEnumerator::MoveNext::est LS: q is divisible by LS[%d] !:((, diviser is: ", 1 + Current().m_label);
1409            dPrint(Current().m_lt, m_rBaseRing, m_rBaseRing, 1);
1410          }
1411#endif
1412//          m_active = true;
1413          return true;
1414        }
1415      }
1416
1417      // the end... :(
1418      assume( m_current == m_finish );
1419     
1420      m_active = false;
1421      return false;
1422    }
1423};
1424
1425
1426
1427bool CReducerFinder::IsDivisible(const poly product) const
1428{
1429  CDivisorEnumerator itr(*this, product);
1430  if( !itr.Reset() )
1431    return false;
1432
1433  return itr.MoveNext();
1434 
1435/* 
1436  const ring& r = m_rBaseRing;
1437 
1438  const long comp = p_GetComp(product, r);
1439  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
1440
1441  assume( comp >= 0 );
1442
1443  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
1444
1445  assume( m_L != NULL ); 
1446
1447  if( it == m_hash.end() )
1448    return false;
1449
1450  const TReducers& reducers = it->second;
1451
1452  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1453  {
1454    assume( m_L->m[(*vit)->m_label] == (*vit)->m_lt );
1455
1456    if( (*vit)->DivisibilityCheck(product, not_sev, r) )
1457    {
1458      if( __DEBUG__ )
1459      {
1460        Print("_FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", 1 + (*vit)->m_label);
1461        dPrint((*vit)->m_lt, r, r, 1);
1462      }
1463
1464      return true;
1465    }
1466  }
1467
1468  return false;
1469*/ 
1470}
1471
1472
1473#ifndef NDEBUG
1474void CReducerFinder::DebugPrint() const
1475{
1476  const ring& r = m_rBaseRing;
1477
1478  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++)
1479  {
1480    Print("Hash Key: %ld, Values: \n", it->first);
1481    const TReducers& reducers = it->second;
1482
1483    for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1484    {
1485      const poly p = (*vit)->m_lt;
1486
1487      assume( p_GetComp(p, r) == it->first );
1488
1489      const int k = (*vit)->m_label;
1490
1491      assume( m_L->m[k] == p );
1492
1493      const unsigned long p_sev = (*vit)->m_sev;
1494
1495      assume( p_sev == p_GetShortExpVector(p, r) );
1496
1497      Print("L[%d]: ", k); dPrint(p, r, r, 0); Print("SEV: %ld\n", p_sev);
1498    }
1499  }
1500}
1501#endif
1502
1503/// TODO:
1504class CDivisorEnumerator2: public SchreyerSyzygyComputationFlags
1505{
1506  private: 
1507    const CReducerFinder& m_reds;
1508    const poly m_multiplier, m_term;
1509    const unsigned long m_not_sev;
1510    const long m_comp;
1511
1512    CReducerFinder::CReducersHash::const_iterator m_itr;
1513    CReducerFinder::TReducers::const_iterator m_current, m_finish;
1514
1515    bool m_active;
1516
1517  public:
1518    CDivisorEnumerator2(const CReducerFinder& self, const poly m, const poly t):
1519        SchreyerSyzygyComputationFlags(self),
1520        m_reds(self),
1521        m_multiplier(m), m_term(t),
1522        m_not_sev(~p_GetShortExpVector(m, t, m_rBaseRing)),
1523        m_comp(p_GetComp(t, m_rBaseRing)),
1524        m_itr(), m_current(), m_finish(),
1525        m_active(false)
1526    {
1527      assume( m_comp >= 0 );
1528      assume( m_reds.m_L != NULL ); 
1529      assume( m_multiplier != NULL ); 
1530      assume( m_term != NULL );
1531//      assume( p_GetComp(m_multiplier, m_rBaseRing) == 0 );
1532    }
1533
1534    inline bool Reset()
1535    {
1536      m_active = false;
1537     
1538      m_itr = m_reds.m_hash.find(m_comp);
1539
1540      if( m_itr == m_reds.m_hash.end() )
1541        return false;
1542
1543      assume( m_itr->first == m_comp );
1544
1545      m_current = (m_itr->second).begin();
1546      m_finish = (m_itr->second).end();
1547
1548      if (m_current == m_finish)
1549        return false;
1550
1551      return true;     
1552    }
1553
1554    const CLeadingTerm& Current() const
1555    {
1556      assume( m_active );
1557      assume( m_current != m_finish );
1558
1559      return *(*m_current);
1560    }
1561 
1562    inline bool MoveNext()
1563    {
1564      assume( m_current != m_finish );
1565
1566      if( m_active )
1567        ++m_current;
1568      else 
1569        m_active = true;
1570   
1571     
1572      // looking for the next good entry
1573      for( ; m_current != m_finish; ++m_current )
1574      {
1575        assume( m_reds.m_L->m[Current().m_label] == Current().m_lt );
1576
1577        if( Current().DivisibilityCheck(m_multiplier, m_term, m_not_sev, m_rBaseRing) )
1578        {
1579#ifndef NDEBUG
1580          if( __DEBUG__ )
1581          {
1582            Print("CDivisorEnumerator::MoveNext::est LS: q is divisible by LS[%d] !:((, diviser is: ", 1 + Current().m_label);
1583            dPrint(Current().m_lt, m_rBaseRing, m_rBaseRing, 1);
1584          }
1585#endif
1586//          m_active = true;
1587          return true;
1588         
1589        }
1590      }
1591
1592      // the end... :(
1593      assume( m_current == m_finish );
1594     
1595      m_active = false;
1596      return false;
1597    }
1598};
1599   
1600poly CReducerFinder::FindReducer(const poly multiplier, const poly t,
1601                                 const poly syzterm,
1602                                 const CReducerFinder& syz_checker) const
1603{
1604  CDivisorEnumerator2 itr(*this, multiplier, t);
1605  if( !itr.Reset() )
1606    return NULL;
1607
1608  // don't care about the module component of multiplier (as it may be the syzygy term)
1609  // product = multiplier * t?
1610  const ring& r = m_rBaseRing;
1611
1612  assume( multiplier != NULL ); assume( t != NULL );
1613
1614  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
1615
1616  long c = 0;
1617
1618  if (syzterm != NULL)
1619    c = p_GetComp(syzterm, r) - 1;
1620
1621  assume( c >= 0 && c < IDELEMS(L) );
1622
1623  if (__DEBUG__ && (syzterm != NULL))
1624  {
1625    const poly m = L->m[c];
1626
1627    assume( m != NULL ); assume( pNext(m) == NULL );
1628
1629    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
1630
1631    poly pr = p_Mult_q( leadmonom(multiplier, r, false), leadmonom(t, r, false), r);
1632   
1633    assume( p_EqualPolys(lm, pr, r) );
1634
1635    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
1636    //  def @@product = leadmonomial(syzterm) * L[@@r];
1637
1638    p_Delete(&lm, r);   
1639    p_Delete(&pr, r);   
1640  }
1641   
1642  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
1643
1644  const poly q = p_New(r); pNext(q) = NULL;
1645
1646  if( __DEBUG__ )
1647    p_SetCoeff0(q, 0, r); // for printing q
1648
1649  while( itr.MoveNext() )
1650  {
1651    const poly p = itr.Current().m_lt;
1652    const int k  = itr.Current().m_label;
1653     
1654    p_ExpVectorSum(q, multiplier, t, r); // q == product == multiplier * t // TODO: do it once?
1655    p_ExpVectorDiff(q, q, p, r); // (LM(product) / LM(L[k]))
1656   
1657    p_SetComp(q, k + 1, r);
1658    p_Setm(q, r);
1659
1660    // cannot allow something like: a*gen(i) - a*gen(i)
1661    if (syzterm != NULL && (k == c))
1662      if (p_ExpVectorEqual(syzterm, q, r))
1663      {
1664#ifndef NDEBUG
1665        if( __DEBUG__ )
1666        {
1667          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1668          dPrint(syzterm, r, r, 1);
1669        }   
1670#endif
1671        continue;
1672      }
1673
1674    // while the complement (the fraction) is not reducible by leading syzygies
1675    if( to_check && syz_checker.IsDivisible(q) ) 
1676    {
1677#ifndef NDEBUG
1678      if( __DEBUG__ )
1679      {
1680        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
1681      }
1682#endif
1683      continue;
1684    }
1685
1686    number n = n_Mult( p_GetCoeff(multiplier, r), p_GetCoeff(t, r), r);
1687    p_SetCoeff0(q, n_Neg( n_Div(n, p_GetCoeff(p, r), r), r), r);
1688    n_Delete(&n, r);
1689   
1690    return q;
1691  }
1692   
1693/*
1694  const long comp = p_GetComp(t, r); assume( comp >= 0 );
1695  const unsigned long not_sev = ~p_GetShortExpVector(multiplier, t, r); // !
1696
1697//   for( int k = IDELEMS(L)-1; k>= 0; k-- )
1698//   {
1699//     const poly p = L->m[k];
1700//
1701//     if ( p_GetComp(p, r) != comp )
1702//       continue;
1703//
1704//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
1705
1706   // looking for an appropriate diviser p = L[k]...
1707  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
1708
1709  if( it == m_hash.end() )
1710    return NULL;
1711
1712  assume( m_L != NULL );
1713
1714  const TReducers& reducers = it->second;
1715
1716  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1717  {
1718
1719    const poly p = (*vit)->m_lt;
1720    const int k = (*vit)->m_label;
1721
1722    assume( L->m[k] == p );
1723
1724//    const unsigned long p_sev = (*vit)->m_sev;
1725//    assume( p_sev == p_GetShortExpVector(p, r) );     
1726
1727//    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
1728//      continue;     
1729
1730    if( !(*vit)->DivisibilityCheck(multiplier, t, not_sev, r) )
1731      continue;
1732   
1733   
1734//    if (p_sev & not_sev) continue;
1735//    if( !_p_LmDivisibleByNoComp(p, multiplier, t, r) ) continue;     
1736
1737
1738    p_ExpVectorSum(q, multiplier, t, r); // q == product == multiplier * t       
1739    p_ExpVectorDiff(q, q, p, r); // (LM(product) / LM(L[k]))
1740   
1741    p_SetComp(q, k + 1, r);
1742    p_Setm(q, r);
1743
1744    // cannot allow something like: a*gen(i) - a*gen(i)
1745    if (syzterm != NULL && (k == c))
1746      if (p_ExpVectorEqual(syzterm, q, r))
1747      {
1748        if( __DEBUG__ )
1749        {
1750          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1751          dPrint(syzterm, r, r, 1);
1752        }   
1753
1754        continue;
1755      }
1756
1757    // while the complement (the fraction) is not reducible by leading syzygies
1758    if( to_check && syz_checker.IsDivisible(q) )
1759    {
1760      if( __DEBUG__ )
1761      {
1762        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
1763      }
1764
1765      continue;
1766    }
1767
1768    number n = n_Mult( p_GetCoeff(multiplier, r), p_GetCoeff(t, r), r);
1769    p_SetCoeff0(q, n_Neg( n_Div(n, p_GetCoeff(p, r), r), r), r);
1770    n_Delete(&n, r);
1771   
1772    return q;
1773  }
1774*/ 
1775
1776  p_LmFree(q, r);
1777
1778  return NULL;
1779 
1780}
1781
1782
1783poly CReducerFinder::FindReducer(const poly product, const poly syzterm, const CReducerFinder& syz_checker) const
1784{
1785  CDivisorEnumerator itr(*this, product);
1786  if( !itr.Reset() )
1787    return NULL;
1788
1789
1790  const ring& r = m_rBaseRing;
1791
1792  assume( product != NULL );
1793
1794  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
1795
1796  long c = 0;
1797
1798  if (syzterm != NULL)
1799    c = p_GetComp(syzterm, r) - 1;
1800
1801  assume( c >= 0 && c < IDELEMS(L) );
1802
1803  if (__DEBUG__ && (syzterm != NULL))
1804  {
1805    const poly m = L->m[c];
1806
1807    assume( m != NULL ); assume( pNext(m) == NULL );
1808
1809    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
1810    assume( p_EqualPolys(lm, product, r) );
1811
1812    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
1813    //  def @@product = leadmonomial(syzterm) * L[@@r];
1814
1815    p_Delete(&lm, r);   
1816  }
1817
1818
1819  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
1820
1821  const poly q = p_New(r); pNext(q) = NULL;
1822
1823  if( __DEBUG__ )
1824    p_SetCoeff0(q, 0, r); // for printing q
1825
1826  while( itr.MoveNext() )
1827  {
1828    const poly p = itr.Current().m_lt;
1829    const int k  = itr.Current().m_label;
1830     
1831    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
1832    p_SetComp(q, k + 1, r);
1833    p_Setm(q, r);
1834
1835    // cannot allow something like: a*gen(i) - a*gen(i)
1836    if (syzterm != NULL && (k == c))
1837      if (p_ExpVectorEqual(syzterm, q, r))
1838      {
1839#ifndef NDEBUG
1840        if( __DEBUG__ )
1841        {
1842          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1843          dPrint(syzterm, r, r, 1);
1844        }   
1845#endif
1846        continue;
1847      }
1848
1849    // while the complement (the fraction) is not reducible by leading syzygies
1850    if( to_check && syz_checker.IsDivisible(q) ) 
1851    {
1852#ifndef NDEBUG
1853      if( __DEBUG__ )
1854      {
1855        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
1856      }
1857#endif
1858      continue;
1859    }
1860
1861    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
1862   
1863    return q;
1864  }
1865   
1866   
1867   
1868/*   
1869  const long comp = p_GetComp(product, r);
1870  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
1871
1872  assume( comp >= 0 );
1873
1874//   for( int k = IDELEMS(L)-1; k>= 0; k-- )
1875//   {
1876//     const poly p = L->m[k];
1877//
1878//     if ( p_GetComp(p, r) != comp )
1879//       continue;
1880//
1881//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
1882 
1883   // looking for an appropriate diviser p = L[k]...
1884  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
1885
1886  if( it == m_hash.end() )
1887    return NULL;
1888
1889  assume( m_L != NULL );
1890
1891  const TReducers& reducers = it->second;
1892
1893  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
1894
1895  const poly q = p_New(r); pNext(q) = NULL;
1896
1897  if( __DEBUG__ )
1898    p_SetCoeff0(q, 0, r); // for printing q
1899 
1900  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1901  {
1902    const poly p = (*vit)->m_lt;
1903
1904    assume( p_GetComp(p, r) == comp );
1905
1906    const int k = (*vit)->m_label;
1907
1908    assume( L->m[k] == p );
1909
1910    const unsigned long p_sev = (*vit)->m_sev;
1911
1912    assume( p_sev == p_GetShortExpVector(p, r) );     
1913
1914    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
1915      continue;     
1916
1917//     // ... which divides the product, looking for the _1st_ appropriate one!
1918//     if( !p_LmDivisibleByNoComp(p, product, r) ) // included inside  p_LmShortDivisibleBy!
1919//       continue;
1920
1921    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
1922    p_SetComp(q, k + 1, r);
1923    p_Setm(q, r);
1924
1925    // cannot allow something like: a*gen(i) - a*gen(i)
1926    if (syzterm != NULL && (k == c))
1927      if (p_ExpVectorEqual(syzterm, q, r))
1928      {
1929        if( __DEBUG__ )
1930        {
1931          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1932          dPrint(syzterm, r, r, 1);
1933        }   
1934
1935        continue;
1936      }
1937
1938    // while the complement (the fraction) is not reducible by leading syzygies
1939    if( to_check && syz_checker.IsDivisible(q) )
1940    {
1941      if( __DEBUG__ )
1942      {
1943        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
1944      }
1945     
1946      continue;
1947    }
1948
1949    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
1950    return q;
1951  }
1952*/
1953
1954  p_LmFree(q, r);
1955
1956  return NULL;
1957}
1958
1959
1960
1961CLCM::CLCM(const ideal& L, const SchreyerSyzygyComputationFlags& flags):
1962    SchreyerSyzygyComputationFlags(flags), std::vector<bool>(),
1963    m_compute(false), m_N(rVar(flags.m_rBaseRing))
1964{
1965  const ring& R = m_rBaseRing;
1966  assume( flags.m_rBaseRing == R );
1967  assume( R != NULL );
1968
1969  assume( L != NULL );
1970
1971  if( __TAILREDSYZ__ && !__HYBRIDNF__ && (L != NULL))
1972  {
1973    const int l = IDELEMS(L);
1974
1975    assume( l > 0 );
1976
1977    resize(l, false);
1978
1979    for( int k = l - 1; k >= 0; k-- )
1980    {
1981      const poly a = L->m[k]; assume( a != NULL );
1982
1983      for (unsigned int j = m_N; j > 0; j--)
1984        if ( !(*this)[j] )
1985          (*this)[j] = (p_GetExp(a, j, R) > 0);
1986    }
1987
1988    m_compute = true;   
1989  }
1990}
1991
1992
1993bool CLCM::Check(const poly m) const
1994{
1995  assume( m != NULL );
1996  if( m_compute && (m != NULL))
1997  { 
1998    const ring& R = m_rBaseRing;
1999
2000    assume( __TAILREDSYZ__ && !__HYBRIDNF__ );
2001
2002    for (unsigned int j = m_N; j > 0; j--)
2003      if ( (*this)[j] )
2004        if(p_GetExp(m, j, R) > 0)
2005          return true;
2006
2007    return false;
2008
2009  } else return true;
2010}
2011
2012
2013
2014
2015END_NAMESPACE               END_NAMESPACE_SINGULARXX
2016
2017
2018// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
Note: See TracBrowser for help on using the repository browser.