source: git/Singular/dyn_modules/syzextra/syzextra.cc @ 196580

spielwiese
Last change on this file since 196580 was 196580, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
Added SYZCHECK to syzextra + some more minor additions
  • Property mode set to 100644
File size: 58.2 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 * New implementations for the computation of syzygies and resolutions
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_NONAME
67
68static inline poly pp_Add_qq( const poly a, const poly b, const ring R)
69{
70  return p_Add_q( p_Copy(a, R), p_Copy(b, R), R );
71}
72
73static inline poly p_VectorProductLT( poly s,  const ideal& L, const ideal& T, const ring& R)
74{
75  assume( IDELEMS(L) == IDELEMS(T) );
76  poly vp = NULL; // resulting vector product
77
78  while( s != NULL )
79  {
80    const poly nxt = pNext(s);
81    pNext(s) = NULL;
82
83    if( !n_IsZero( p_GetCoeff(s, R), R) )
84    {
85      const int i = p_GetComp(s, R) - 1;
86      assume( i >= 0 ); assume( i < IDELEMS(L) );
87      p_SetComp(s, 0, R); p_SetmComp(s, R);
88
89      vp = p_Add_q( vp, pp_Mult_qq( s, L->m[i], R ), R); 
90      vp = p_Add_q( vp, pp_Mult_qq( s, T->m[i], R ), R); 
91    }
92
93    p_Delete(&s, R);
94
95    s = nxt;
96  };
97
98  assume( s == NULL );
99
100  return vp;
101}
102
103static inline int atGetInt(idhdl rootRingHdl, const char* attribute, long def)
104{
105  return ((int)(long)(atGet(rootRingHdl, attribute, INT_CMD, (void*)def)));
106}
107
108END_NAMESPACE
109
110BEGIN_NAMESPACE(SORT_c_ds)
111
112#if (defined(HAVE_QSORT_R) && (defined __APPLE__ || defined __MACH__ || defined __DARWIN__ || defined __FreeBSD__ || defined __BSD__ || defined OpenBSD3_1 || defined OpenBSD3_9))
113static int cmp_c_ds(void *R, const void *p1, const void *p2){
114#elif (defined(HAVE_QSORT_R) && (defined _GNU_SOURCE || defined __GNU__ || defined __linux__))
115static int cmp_c_ds(const void *p1, const void *p2, void *R){
116#else
117static int cmp_c_ds(const void *p1, const void *p2){ void *R = currRing;
118#endif
119  assume(R != NULL);
120  const int YES = 1;
121  const int NO = -1;
122
123  const ring r =  (const ring) R; // TODO/NOTE: the structure is known: C, lp!!!
124
125  assume( r == currRing ); // for now...
126
127  const poly a = *(const poly*)p1;
128  const poly b = *(const poly*)p2;
129
130  assume( a != NULL );
131  assume( b != NULL );
132
133  p_LmTest(a, r);
134  p_LmTest(b, r);
135
136
137  const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r);
138
139  // TODO: test this!!!!!!!!!!!!!!!!
140
141  //return -( compare (c, qsorts) )
142
143#ifndef SING_NDEBUG
144  const int __DEBUG__ = 0;
145  if( __DEBUG__ )
146  {
147    PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2);
148    PrintS("b: "); dPrint(b, r, r, 2);
149    PrintLn();
150  }
151#endif
152
153
154  if( iCompDiff > 0 )
155    return YES;
156
157  if( iCompDiff < 0 )
158    return  NO;
159
160  assume( iCompDiff == 0 );
161
162  const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r);
163
164  if( iDegDiff > 0 )
165    return YES;
166
167  if( iDegDiff < 0 )
168    return  NO;
169
170  assume( iDegDiff == 0 );
171
172#ifndef SING_NDEBUG
173  if( __DEBUG__ )
174  {
175    PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn();
176  }
177#endif
178
179  for (int v = rVar(r); v > 0; v--)
180  {
181    assume( v > 0 );
182    assume( v <= rVar(r) );
183
184    const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r);
185
186    if( d > 0 )
187      return YES;
188
189    if( d < 0 )
190      return NO;
191
192    assume( d == 0 );
193  }
194
195  return 0;
196}
197
198/*
199static int cmp_poly(const poly &a, const poly &b)
200{
201  const int YES = 1;
202  const int NO = -1;
203
204  const ring r =  (const ring) currRing; // TODO/NOTE: the structure is known: C, lp!!!
205
206  assume( r == currRing );
207
208  assume( a != NULL );
209  assume( b != NULL );
210
211  p_LmTest(a, r);
212  p_LmTest(b, r);
213  assume( p_GetComp(a, r) == 0 );
214  assume( p_GetComp(b, r) == 0 );
215
216#ifndef SING_NDEBUG
217  const int __DEBUG__ = 0;
218  if( __DEBUG__ )
219  {
220    PrintS("cmp_lex: a, b: \np1: "); dPrint(a, r, r, 2);
221    PrintS("b: "); dPrint(b, r, r, 2);
222    PrintLn();
223  }
224#endif
225
226  for (int v = rVar(r); v > 0; v--)
227  {
228    assume( v > 0 );
229    assume( v <= rVar(r) );
230
231    const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r);
232
233    if( d > 0 )
234      return YES;
235
236    if( d < 0 )
237      return NO;
238
239    assume( d == 0 );
240  }
241
242  return 0;
243}
244*/
245
246END_NAMESPACE
247/* namespace SORT_c_ds */
248
249/// writes a monomial (p),
250/// uses form x*gen(.) if ko != coloumn number of p
251static void writeLatexTerm(const poly t, const ring r, const bool bCurrSyz = true, const bool bLTonly = true)
252{
253  if( t == NULL )
254  {
255    PrintS(" 0 ");
256    return;
257  }
258
259  assume( r != NULL );
260  const coeffs C = r->cf; assume(C != NULL);
261
262  poly p = t;
263  BOOLEAN writePlus = FALSE;
264
265  do {
266  assume( p != NULL );
267
268  // write coef...
269  number& n = p_GetCoeff(p, r);
270
271  n_Normalize(n, C);
272
273  BOOLEAN writeMult = FALSE; ///< needs * before pp or module generator
274
275  BOOLEAN writeOne = FALSE; ///< need to write something after '-'!
276
277  if( n_IsZero(n, C) )
278  {
279    PrintS( writePlus? " + 0" : " 0 " ); writePlus = TRUE; writeMult = TRUE;
280//    return; // yes?
281  }
282
283  if (n_IsMOne(n, C))
284  {
285    PrintS(" - "); writeOne = TRUE; writePlus = FALSE;
286  }
287  else if (!n_IsOne(n, C))
288  {
289    if( writePlus && n_GreaterZero(n, C) )
290      PrintS(" + ");
291
292    StringSetS(""); n_WriteLong(n, C);
293    if (true)
294    {
295      char *s = StringEndS(); PrintS(s); omFree(s);
296    }
297
298
299    writeMult = TRUE;
300    writePlus = TRUE;
301  } else
302     writeOne = TRUE;
303
304  // missing '1' if no PP and gen...!?
305  // write monom...
306  const short N = rVar(r);
307
308  BOOLEAN wrotePP = FALSE; ///< needs * before module generator?
309
310  for (short i = 0; i < N; i++)
311  {
312    const long ee = p_GetExp(p, i+1, r);
313
314    if (ee!=0L)
315    {
316      if (writeMult)
317      {
318        PrintS(" ");
319        writeMult = FALSE;
320      } else
321      if( writePlus )
322        PrintS(" + ");
323
324      writePlus = FALSE;
325
326      if (ee != 1L)
327        Print(" %s^{%ld} ", rRingVar(i, r), ee);
328      else
329        Print(" %s ", rRingVar(i, r));
330
331      writeOne = FALSE;
332      wrotePP = TRUE;
333    }
334  }
335
336  writePlus = writePlus || wrotePP;
337  writeMult = writeMult || wrotePP;
338
339  // write module gen...
340  const long comp = p_GetComp(p, r);
341
342  if (comp > 0 )
343  {
344    if (writeMult)
345      PrintS("  ");
346     else
347      if( writePlus )
348        PrintS(" + ");
349
350    if (bCurrSyz)
351      Print(" \\\\GEN{%ld} ", comp);
352    else
353      Print(" \\\\GENP{%ld} ", comp);
354
355      writeOne = FALSE;
356  }
357
358  if ( writeOne )
359    PrintS( writePlus? " + 1 " : " 1 " );
360
361
362  pIter(p);
363
364  writePlus = TRUE;
365  } while( (!bLTonly) && (p != NULL) );
366
367}
368
369
370
371
372
373/// return a new term: leading coeff * leading monomial of p
374/// with 0 leading component!
375poly leadmonom(const poly p, const ring r, const bool bSetZeroComp)
376{
377  poly m = NULL;
378
379  if( p != NULL )
380  {
381    assume( p != NULL );
382    p_LmTest(p, r);
383
384    m = p_LmInit(p, r);
385    p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);
386
387    if( bSetZeroComp )
388      p_SetComp(m, 0, r);
389    p_Setm(m, r);
390
391
392    assume( m != NULL );
393    assume( pNext(m) == NULL );
394    p_LmTest(m, r);
395
396    if( bSetZeroComp )
397      assume( p_GetComp(m, r) == 0 );
398  }
399
400  return m;
401}
402
403
404
405poly p_Tail(const poly p, const ring r)
406{
407  if( p == NULL)
408    return NULL;
409  else
410    return p_Copy( pNext(p), r );
411}
412
413
414ideal id_Tail(const ideal id, const ring r)
415{
416  if( id == NULL)
417    return NULL;
418
419  const ideal newid = idInit(IDELEMS(id),id->rank);
420
421  for (int i=IDELEMS(id) - 1; i >= 0; i--)
422    newid->m[i] = p_Tail( id->m[i], r );
423
424  newid->rank = id_RankFreeModule(newid, currRing);
425
426  return newid;
427}
428
429
430
431void Sort_c_ds(const ideal id, const ring r)
432{
433  const int sizeNew = IDELEMS(id);
434
435#if ( (defined(HAVE_QSORT_R)) && (defined __APPLE__ || defined __MACH__ || defined __DARWIN__ || defined __FreeBSD__ || defined __BSD__ || defined OpenBSD3_1 || defined OpenBSD3_9) )
436#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, r, cmp)
437#elif ( (defined(HAVE_QSORT_R)) && (defined _GNU_SOURCE || defined __GNU__ || defined __linux__))
438#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r)
439#else
440#define qsort_my(m, s, ss, r, cmp) qsort(m, s, ss, cmp)
441#endif
442
443  if( sizeNew >= 2 )
444    qsort_my(id->m, sizeNew, sizeof(poly), r, FROM_NAMESPACE(SORT_c_ds, cmp_c_ds));
445
446#undef qsort_my
447
448  id->rank = id_RankFreeModule(id, r);
449}
450
451/// Clean up all the accumulated data
452void SchreyerSyzygyComputation::CleanUp()
453{
454  extern void id_Delete (ideal*, const ring);
455
456  id_Delete(const_cast<ideal*>(&m_idTails), m_rBaseRing); // TODO!!!
457
458  if( m_sum_bucket != NULL )
459  {
460    sBucketDestroy(&m_sum_bucket);
461    m_sum_bucket = NULL;
462  }
463
464  if( m_spoly_bucket != NULL )
465  {
466    kBucketDestroy(&m_spoly_bucket);
467    m_spoly_bucket = NULL;
468  }
469
470  for( TCache::iterator it = m_cache.begin(); it != m_cache.end(); it++ )
471  {
472    TP2PCache& T = it->second;
473
474    for(TP2PCache::iterator vit = T.begin(); vit != T.end(); vit++ )
475    {
476      p_Delete( (&(vit->second)), m_rBaseRing);
477      p_Delete( const_cast<poly*>(&(vit->first)), m_rBaseRing);
478    }
479  }
480}
481  /*
482  for( TTailTerms::const_iterator it = m_idTailTerms.begin(); it != m_idTailTerms.end(); it++ )
483  {
484    const TTail& v = *it;
485    for(TTail::const_iterator vit = v.begin(); vit != v.end(); vit++ )
486      delete const_cast<CTailTerm*>(*vit);
487  }
488  */
489
490
491
492int CReducerFinder::PreProcessTerm(const poly t, CReducerFinder& syzChecker) const
493{
494  assume( t != NULL );
495
496  if( __DEBUG__ & __TAILREDSYZ__ )
497    assume( !IsDivisible(t) ); // each input term should NOT be in <L>
498
499  const ring r = m_rBaseRing;
500
501
502  if( __TAILREDSYZ__ )
503    if( p_LmIsConstant(t, r) ) // most basic case of baing coprime with L, whatever that is...
504      return 1; // TODO: prove this...?
505
506  //   return false; // appears to be fine
507
508  const long comp = p_GetComp(t, r);
509
510  CReducersHash::const_iterator itr = m_hash.find(comp);
511
512  if ( itr == m_hash.end() )
513    return 2; // no such leading component!!!
514
515  assume( itr->first == comp );
516
517  const bool bIdealCase = (comp == 0);
518  const bool bSyzCheck = syzChecker.IsNonempty(); // need to check even in ideal case????? proof?  "&& !bIdealCase"
519
520  if( __TAILREDSYZ__ && (bIdealCase || bSyzCheck) )
521  {
522    const TReducers& v = itr->second;
523    const int N = rVar(r);
524    // TODO: extract exps of t beforehand?!
525    bool coprime = true;
526    for(TReducers::const_iterator vit = v.begin(); (vit != v.end()) && coprime; ++vit )
527    {
528      assume( m_L->m[(*vit)->m_label] == (*vit)->m_lt );
529
530      const poly p = (*vit)->m_lt;
531
532      assume( p_GetComp(p, r) == comp );
533
534      // TODO: check if coprime with Leads... if __TAILREDSYZ__ !
535      for( int var = N; var > 0; --var )
536        if( (p_GetExp(p, var, r) != 0) && (p_GetExp(t, var, r) != 0) )
537        {
538#ifndef SING_NDEBUG
539          if( __DEBUG__ | 0)
540          {
541            PrintS("CReducerFinder::PreProcessTerm, 't' is NOT co-prime with the following leading term: \n");
542            dPrint(p, r, r, 1);
543          }
544#endif
545          coprime = false; // t not coprime with p!
546          break;
547        }
548
549      if( bSyzCheck && coprime )
550      {
551        poly ss = p_LmInit(t, r);
552        p_SetCoeff0(ss, n_Init(1, r), r); // for delete & printout only!...
553        p_SetComp(ss, (*vit)->m_label + 1, r); // coeff?
554        p_Setm(ss, r);
555
556        coprime = ( syzChecker.IsDivisible(ss) );
557
558#ifndef SING_NDEBUG
559        if( __DEBUG__ && !coprime)
560        {
561          PrintS("CReducerFinder::PreProcessTerm, 't' is co-prime with p but may lead to NOT divisible syz.term: \n");
562          dPrint(ss, r, r, 1);
563        }
564#endif
565
566        p_LmDelete(&ss, r); // deletes coeff as well???
567      }
568
569    }
570
571#ifndef SING_NDEBUG
572    if( __DEBUG__ && coprime )
573      PrintS("CReducerFinder::PreProcessTerm, the following 't' is 'co-prime' with all of leading terms! \n");
574#endif
575
576    return coprime? 3: 0; // t was coprime with all of leading terms!!!
577
578  }
579  //   return true; // delete the term
580
581  return 0;
582}
583
584
585void SchreyerSyzygyComputation::SetUpTailTerms()
586{
587  const ideal idTails = m_idTails;
588  assume( idTails != NULL );
589  assume( idTails->m != NULL );
590  const ring r = m_rBaseRing;
591
592#ifndef SING_NDEBUG
593  if( __DEBUG__ | 0)
594  {
595    PrintS("SchreyerSyzygyComputation::SetUpTailTerms(): Tails: \n");
596    dPrint(idTails, r, r, 0);
597  }
598
599  unsigned long pp[4] = {0,0,0,0}; // count preprocessed terms...
600#endif
601
602  for( int p = IDELEMS(idTails) - 1; p >= 0; --p )
603    for( poly* tt = &(idTails->m[p]); (*tt) != NULL;  )
604    {
605      const poly t = *tt;
606      const int k = m_div.PreProcessTerm(t, m_checker); // 0..3
607      assume( 0 <= k && k <= 3 );
608
609#ifndef SING_NDEBUG
610      pp[k]++;
611#endif
612
613      if( k )
614      {
615#ifndef SING_NDEBUG
616        if( __DEBUG__)
617        {
618          Print("SchreyerSyzygyComputation::SetUpTailTerms(): PP (%d) the following TT: \n", k);
619          dPrint(t, r, r, 1);
620        }
621#endif
622
623        (*tt) = p_LmDeleteAndNext(t, r); // delete the lead and next...
624      }
625      else
626        tt = &pNext(t); // go next?
627
628    }
629
630#ifndef SING_NDEBUG
631  if( !__TREEOUTPUT__ )
632  if( TEST_OPT_PROT | 1)
633    Print("%%      **!!**      SchreyerSyzygyComputation::SetUpTailTerms()::PreProcessing(): X: {c: %lu, C: %lu, P: %lu} + %lu\n", pp[1], pp[2], pp[3], pp[0]);
634#endif
635
636#ifndef SING_NDEBUG
637  if( !__TREEOUTPUT__ )
638  if( __DEBUG__ | 0)
639  {
640    PrintS("SchreyerSyzygyComputation::SetUpTailTerms(): Preprocessed Tails: \n");
641    dPrint(idTails, r, r, 0);
642  }
643#endif
644
645}
646/*
647  m_idTailTerms.resize( IDELEMS(idTails) );
648
649  for( unsigned int p = IDELEMS(idTails) - 1; p >= 0; p -- )
650  {
651    TTail& v = m_idTailTerms[p];
652    poly t = idTails->m[p];
653    v.resize( pLength(t) );
654
655    unsigned int pp = 0;
656
657    while( t != NULL )
658    {
659      assume( t != NULL );
660      // TODO: compute L:t!
661//      ideal reducers;
662//      CReducerFinder m_reducers
663
664      CTailTerm* d = v[pp] = new CTailTerm();
665
666      ++pp; pIter(t);
667    }
668  }
669*/
670
671
672
673ideal SchreyerSyzygyComputation::Compute1LeadingSyzygyTerms()
674{
675  const ideal& id = m_idLeads;
676  const ring& r = m_rBaseRing;
677//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
678
679  assume(!__LEAD2SYZ__);
680
681  // 1. set of components S?
682  // 2. for each component c from S: set of indices of leading terms
683  // with this component?
684  // 3. short exp. vectors for each leading term?
685
686  const int size = IDELEMS(id);
687
688  if( size < 2 )
689  {
690    const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal...
691    return newid;
692  }
693
694  // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!??
695
696  // components should come in groups: count elements in each group
697  // && estimate the real size!!!
698
699
700  // use just a vector instead???
701  const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
702
703  int k = 0;
704
705  for (int j = 0; j < size; j++)
706  {
707    const poly p = id->m[j];
708    assume( p != NULL );
709    const int  c = p_GetComp(p, r);
710
711    for (int i = j - 1; i >= 0; i--)
712    {
713      const poly pp = id->m[i];
714      assume( pp != NULL );
715      const int  cc = p_GetComp(pp, r);
716
717      if( c != cc )
718        continue;
719
720      const poly m = p_Init(r); // p_New???
721
722      // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)!
723      for (int v = rVar(r); v > 0; v--)
724      {
725        assume( v > 0 );
726        assume( v <= rVar(r) );
727
728        const short e1 = p_GetExp(p , v, r);
729        const short e2 = p_GetExp(pp, v, r);
730
731        if( e1 >= e2 )
732          p_SetExp(m, v, 0, r);
733        else
734          p_SetExp(m, v, e2 - e1, r);
735
736      }
737
738      assume( (j > i) && (i >= 0) );
739
740      p_SetComp(m, j + 1, r);
741      pNext(m) = NULL;
742      p_SetCoeff0(m, n_Init(1, r->cf), r); // for later...
743
744      p_Setm(m, r); // should not do anything!!!
745
746      newid->m[k++] = m;
747    }
748  }
749
750//   if( __DEBUG__ && FALSE )
751//   {
752//     PrintS("ComputeLeadingSyzygyTerms::Temp0: \n");
753//     dPrint(newid, r, r, 1);
754//   }
755
756  // the rest of newid is assumed to be zeroes...
757
758  // simplify(newid, 2 + 32)??
759  // sort(newid, "C,ds")[1]???
760  id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
761
762//   if( __DEBUG__ && FALSE )
763//   {
764//     PrintS("ComputeLeadingSyzygyTerms::Temp1: \n");
765//     dPrint(newid, r, r, 1);
766//   }
767
768  idSkipZeroes(newid); // #define SIMPL_NULL 2
769
770//   if( __DEBUG__ )
771//   {
772//     PrintS("ComputeLeadingSyzygyTerms::Output: \n");
773//     dPrint(newid, r, r, 1);
774//   }
775
776  Sort_c_ds(newid, r);
777
778  return newid;
779}
780
781ideal SchreyerSyzygyComputation::Compute2LeadingSyzygyTerms()
782{
783  const ideal& id = m_idLeads;
784  const ring& r = m_rBaseRing;
785//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
786
787  // 1. set of components S?
788  // 2. for each component c from S: set of indices of leading terms
789  // with this component?
790  // 3. short exp. vectors for each leading term?
791
792  const int size = IDELEMS(id);
793
794  if( size < 2 )
795  {
796    const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...
797    return newid;
798  }
799
800
801  // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
802
803  // components should come in groups: count elements in each group
804  // && estimate the real size!!!
805
806
807  // use just a vector instead???
808  ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
809
810  int k = 0;
811
812  for (int j = 0; j < size; j++)
813  {
814    const poly p = id->m[j];
815    assume( p != NULL );
816    const int  c = p_GetComp(p, r);
817
818    for (int i = j - 1; i >= 0; i--)
819    {
820      const poly pp = id->m[i];
821      assume( pp != NULL );
822      const int  cc = p_GetComp(pp, r);
823
824      if( c != cc )
825        continue;
826
827        // allocate memory & zero it out!
828      const poly m = p_Init(r); const poly mm = p_Init(r);
829
830
831        // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
832        // TODO: optimize: knowing the ring structure: (C/lp)!
833
834      for (int v = rVar(r); v > 0; v--)
835      {
836        assume( v > 0 );
837        assume( v <= rVar(r) );
838
839        const short e1 = p_GetExp(p , v, r);
840        const short e2 = p_GetExp(pp, v, r);
841
842        if( e1 >= e2 )
843          p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
844        else
845          p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
846
847      }
848
849      assume( (j > i) && (i >= 0) );
850
851      p_SetComp(m, j + 1, r);
852      p_SetComp(mm, i + 1, r);
853
854      const number& lc1 = p_GetCoeff(p , r);
855      const number& lc2 = p_GetCoeff(pp, r);
856
857      number g = n_Lcm( lc1, lc2, r->cf );
858
859      p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
860      p_SetCoeff0(mm, n_InpNeg(n_Div(g, lc2, r), r), r);
861
862      n_Delete(&g, r);
863
864      p_Setm(m, r); // should not do anything!!!
865      p_Setm(mm, r); // should not do anything!!!
866
867      pNext(m) = mm; //        pNext(mm) = NULL;
868
869      newid->m[k++] = m;
870    }
871  }
872
873//   if( __DEBUG__ && FALSE )
874//   {
875//     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
876//     dPrint(newid, r, r, 1);
877//   }
878
879  if( !__TAILREDSYZ__ )
880  {
881      // simplify(newid, 2 + 32)??
882      // sort(newid, "C,ds")[1]???
883    id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
884
885//     if( __DEBUG__ && FALSE )
886//     {
887//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
888//       dPrint(newid, r, r, 1);
889//     }
890  }
891  else
892  {
893      //      option(redSB); option(redTail);
894      //      TEST_OPT_REDSB
895      //      TEST_OPT_REDTAIL
896    assume( r == currRing );
897
898    BITSET _save_test; SI_SAVE_OPT1(_save_test);
899    SI_RESTORE_OPT1(Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB) | _save_test);
900
901    intvec* w=new intvec(IDELEMS(newid));
902    ideal tmp = kStd(newid, currRing->qideal, isHomog, &w);
903    delete w;
904
905    SI_RESTORE_OPT1(_save_test)
906
907    id_Delete(&newid, r);
908    newid = tmp;
909
910//     if( __DEBUG__ && FALSE )
911//     {
912//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
913//       dPrint(newid, r, r, 1);
914//     }
915
916  }
917
918  idSkipZeroes(newid);
919
920  Sort_c_ds(newid, r);
921
922  return newid;
923}
924
925poly SchreyerSyzygyComputation::TraverseNF(const poly a, const poly a2) const
926{
927  const ideal& L = m_idLeads;
928  const ideal& T = m_idTails;
929
930  const ring& R = m_rBaseRing;
931
932  const int r = p_GetComp(a, R) - 1;
933
934  assume( r >= 0 && r < IDELEMS(T) );
935  assume( r >= 0 && r < IDELEMS(L) );
936
937  assume( a != NULL );
938
939  if( __TREEOUTPUT__ )
940  {
941     PrintS("{ \"nodelabel\": \""); writeLatexTerm(a, R); PrintS("\", \"children\": [");
942  }
943
944
945  poly aa = leadmonom(a, R); assume( aa != NULL); // :(
946
947  poly t = TraverseTail(aa, r);
948
949  if( a2 != NULL )
950  {
951    assume( __LEAD2SYZ__ );
952
953    if( __TREEOUTPUT__ )
954    {
955       PrintS("{ \"nodelabel\": \""); writeLatexTerm(a2, R); PrintS("\", \"children\": [");
956    }
957
958    // replace the following... ?
959    const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
960
961    assume( r2 >= 0 && r2 < IDELEMS(T) );
962
963    t = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, r2), R), R);
964
965    p_Delete(&aa2, R);
966
967    if( __TREEOUTPUT__ )
968       PrintS("]},");
969
970  } else
971    t = p_Add_q(t, ReduceTerm(aa, L->m[r], a), R);
972
973  p_Delete(&aa, R);
974
975  // TODO: print t???
976
977  if( __TREEOUTPUT__ )
978  {
979     PrintS("], \"noderesult\": \""); writeLatexTerm(t, R, true, false); PrintS("\" },");
980  }
981  return t;
982}
983
984
985void SchreyerSyzygyComputation::ComputeSyzygy()
986{
987  assume( m_idLeads != NULL );
988  assume( m_idTails != NULL );
989
990  const ideal& L = m_idLeads;
991  const ideal& T = m_idTails;
992
993  ideal& TT = m_syzTails;
994  const ring& R = m_rBaseRing;
995
996  if( m_sum_bucket == NULL )
997    m_sum_bucket = sBucketCreate(R);
998
999  if( m_spoly_bucket == NULL )
1000    m_spoly_bucket = kBucketCreate(R);
1001
1002
1003  assume( IDELEMS(L) == IDELEMS(T) );
1004#ifndef SING_NDEBUG
1005  int t, r;
1006#endif
1007
1008  if( __TREEOUTPUT__ )
1009  {
1010    Print("\n{ \"syzygylayer\": \"%d\", \"hybridnf\": \"%d\", \"diagrams\": \n[", __SYZNUMBER__, __HYBRIDNF__ );
1011  }
1012
1013  if( m_syzLeads == NULL )
1014  {
1015#ifndef SING_NDEBUG
1016    if( !__TREEOUTPUT__ )
1017    if( TEST_OPT_PROT | 1)
1018    {
1019      t = getTimer(); r = getRTimer();
1020      Print("%% %5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::ComputeLeadingSyzygyTerms: t: %d, r: %d\n", r, t, r);
1021    }
1022#endif
1023    ComputeLeadingSyzygyTerms( __LEAD2SYZ__ && !__IGNORETAILS__ ); // 2 terms OR 1 term!
1024#ifndef SING_NDEBUG
1025    if( !__TREEOUTPUT__ )
1026    if( TEST_OPT_PROT | 1)
1027    {
1028      t = getTimer() - t; r = getRTimer() - r;
1029      Print("%% %5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::ComputeLeadingSyzygyTerms: dt: %d, dr: %d\n", getRTimer(), t, r);
1030    }
1031#endif
1032
1033  }
1034
1035  assume( m_syzLeads != NULL );
1036  ideal& LL = m_syzLeads;
1037  const int size = IDELEMS(LL);
1038
1039  TT = idInit(size, 0);
1040
1041  if( size == 1 && LL->m[0] == NULL )
1042  {
1043     if( __TREEOUTPUT__ )
1044       PrintS("]},");
1045     return;
1046  }
1047
1048
1049  // use hybrid method?
1050  const bool method = (__HYBRIDNF__  == 1) || (__HYBRIDNF__ == 2 && __SYZNUMBER__ < 3);
1051
1052  if(  !__IGNORETAILS__)
1053  {
1054    if( T != NULL )
1055    {
1056#ifndef SING_NDEBUG
1057      if( !__TREEOUTPUT__ )
1058      if( TEST_OPT_PROT | 1 )
1059      {
1060        t = getTimer(); r = getRTimer();
1061        Print("%% %5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SetUpTailTerms(): t: %d, r: %d\n", r, t, r);
1062      }
1063#endif
1064
1065      SetUpTailTerms();
1066#ifndef SING_NDEBUG
1067      if( !__TREEOUTPUT__ )
1068      if( TEST_OPT_PROT | 1)
1069      {
1070        t = getTimer() - t; r = getRTimer() - r;
1071        Print("%% %5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SetUpTailTerms(): dt: %d, dr: %d\n", getRTimer(), t, r);
1072      }
1073#endif
1074    }
1075  }
1076
1077#ifndef SING_NDEBUG
1078  if( !__TREEOUTPUT__ )
1079  if( TEST_OPT_PROT | 1)
1080  {
1081    t = getTimer(); r = getRTimer();
1082    Print("%% %5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SyzygyLift: t: %d, r: %d\n", r, t, r);
1083  }
1084#endif
1085
1086  for( int k = size - 1; k >= 0; k-- )
1087  {
1088    const poly a = LL->m[k]; assume( a != NULL );
1089
1090    poly a2 = pNext(a);
1091
1092    // Splitting 2-terms Leading syzygy module
1093    if( a2 != NULL )
1094      pNext(a) = NULL;
1095
1096    if( __IGNORETAILS__ )
1097    {
1098      TT->m[k] = NULL;
1099
1100      assume( a2 != NULL );
1101
1102      if( a2 != NULL )
1103        p_Delete(&a2, R);
1104
1105      continue;
1106    }
1107
1108    //    TT->m[k] = a2;
1109
1110    if( method )
1111      TT->m[k] = SchreyerSyzygyNF(a, a2);
1112    else
1113      TT->m[k] = TraverseNF(a, a2); // TODO: WRONG
1114
1115    if( __SYZCHECK__ )     
1116    {
1117      // TODO: check the correctness (syzygy property): a + TT->m[k] should be a syzygy!!!
1118
1119      poly s = pp_Add_qq( a, TT->m[k], R); // current syzygy
1120
1121      poly vp = p_VectorProductLT(s, L, T, R); 
1122
1123      if( (__DEBUG__ | __TREEOUTPUT__) && (vp != NULL) )
1124      {
1125        Warn("SchreyerSyzygyComputation::ComputeSyzygy: failed syzygy property for syzygy [%d], non-zero image is as follows: ", k);       
1126        dPrint(vp, R, R, 0);       p_Delete(&vp, R);
1127       
1128        PrintS("SchreyerSyzygyComputation::ComputeSyzygy: Wrong syzygy is as follows: ");
1129        s = pp_Add_qq( a, TT->m[k], R);
1130        dPrint(s, R, R, 0); p_Delete(&s, R);
1131
1132        PrintS("SchreyerSyzygyComputation::ComputeSyzygy: Testing with the other method");
1133       
1134        if( !method )
1135          s = SchreyerSyzygyNF(a, a2);
1136        else
1137          s = TraverseNF(a, a2);
1138
1139        s = p_Add_q( a, s, R); // another syzygy
1140        PrintS("SchreyerSyzygyComputation::ComputeSyzygy: The other method gives the following  syzygy: ");
1141        dPrint(s, R, R, 0); 
1142       
1143        vp = p_VectorProductLT(s, L, T, R);
1144       
1145        if( vp == NULL )
1146        {
1147          PrintS("SchreyerSyzygyComputation::ComputeSyzygy: .... which is correct!!! ");
1148        } else
1149        {
1150          Warn("SchreyerSyzygyComputation::ComputeSyzygy: failed to compute syzygy tail[%d] with both methods!!! Non-zero image (2nd) is as follows: ", k);       
1151          dPrint(vp, R, R, 0);
1152        }
1153       
1154      } else
1155        assume( vp == NULL );
1156
1157      p_Delete(&vp, R);
1158    }
1159
1160  }
1161
1162#ifndef SING_NDEBUG
1163  if( !__TREEOUTPUT__ )
1164  if( TEST_OPT_PROT | 1)
1165  {
1166    t = getTimer() - t; r = getRTimer() - r;
1167    Print("%% %5d **!TIME4!** SchreyerSyzygyComputation::ComputeSyzygy::SyzygyLift: dt: %d, dr: %d\n", getRTimer(), t, r);
1168  }
1169#endif
1170
1171  TT->rank = id_RankFreeModule(TT, R);
1172
1173  if( __TREEOUTPUT__ )
1174    PrintS("\n]},");
1175}
1176
1177void SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms(bool bComputeSecondTerms)
1178{
1179//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
1180
1181//  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
1182//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
1183
1184  assume( m_syzLeads == NULL );
1185
1186  if( bComputeSecondTerms )
1187  {
1188    assume( __LEAD2SYZ__ );
1189//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _Compute2LeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
1190    m_syzLeads = Compute2LeadingSyzygyTerms();
1191  }
1192  else
1193  {
1194    assume( !__LEAD2SYZ__ );
1195
1196    m_syzLeads = Compute1LeadingSyzygyTerms();
1197  }
1198//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _ComputeLeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
1199
1200  // NOTE: set m_LS if tails are to be reduced!
1201  assume( m_syzLeads!= NULL );
1202
1203  if (__TAILREDSYZ__ && !__IGNORETAILS__ && (IDELEMS(m_syzLeads) > 0) && !((IDELEMS(m_syzLeads) == 1) && (m_syzLeads->m[0] == NULL)))
1204  {
1205    m_LS = m_syzLeads;
1206    m_checker.Initialize(m_syzLeads);
1207#ifndef SING_NDEBUG
1208    if( __DEBUG__ )
1209    {
1210      const ring& r = m_rBaseRing;
1211      PrintS("SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms: \n");
1212      PrintS("m_syzLeads: \n");
1213      dPrint(m_syzLeads, r, r, 1);
1214      PrintS("m_checker.Initialize(m_syzLeads) => \n");
1215      m_checker.DebugPrint();
1216    }
1217#endif
1218    assume( m_checker.IsNonempty() ); // TODO: this always fails... BUG????
1219  }
1220}
1221
1222#define NOPRODUCT 1
1223
1224poly SchreyerSyzygyComputation::SchreyerSyzygyNF(const poly syz_lead, poly syz_2) const
1225{
1226
1227  assume( !__IGNORETAILS__ );
1228
1229  const ideal& L = m_idLeads;
1230  const ideal& T = m_idTails;
1231  const ring& r = m_rBaseRing;
1232
1233  assume( syz_lead != NULL );
1234
1235  if( __TREEOUTPUT__ )
1236  {
1237//     PrintS("%%%% BEGIN LIFTHYBRID DIAGRAMM\n");
1238    PrintS("{   \"nodelabel\": \""); writeLatexTerm(syz_lead, r); PrintS("\", \"children\": [");
1239  }
1240
1241  if( syz_2 == NULL )
1242  {
1243    const int rr = p_GetComp(syz_lead, r) - 1;
1244
1245    assume( rr >= 0 && rr < IDELEMS(T) );
1246    assume( rr >= 0 && rr < IDELEMS(L) );
1247
1248#if NOPRODUCT
1249    syz_2 = m_div.FindReducer(syz_lead, L->m[rr], syz_lead, m_checker);
1250
1251    if( __TREEOUTPUT__ )
1252    {
1253      PrintS("{ \"nodelabel\": \""); writeLatexTerm(syz_2, r); PrintS("\" },");
1254    }
1255#else
1256    poly aa = leadmonom(syz_lead, r); assume( aa != NULL); // :(
1257    aa = p_Mult_mm(aa, L->m[rr], r);
1258
1259    if( __TREEOUTPUT__ )
1260    {
1261      PrintS("{ \"nodelabel\": \""); writeLatexTerm(syz_2, r); PrintS("\", \"edgelable\": \""); writeLatexTerm(aa, r, false); PrintS("\" },");
1262    }
1263
1264    syz_2 = m_div.FindReducer(aa, syz_lead, m_checker);
1265
1266    p_Delete(&aa, r);
1267#endif
1268
1269  }
1270
1271  assume( syz_2 != NULL ); // by construction of S-Polynomial
1272
1273  assume( L != NULL );
1274  assume( T != NULL );
1275
1276  assume( IDELEMS(L) == IDELEMS(T) );
1277
1278  int  c = p_GetComp(syz_lead, r) - 1;
1279
1280  assume( c >= 0 && c < IDELEMS(T) );
1281
1282  if( m_sum_bucket == NULL )
1283    m_sum_bucket = sBucketCreate(r);
1284
1285  if( m_spoly_bucket == NULL )
1286    m_spoly_bucket = kBucketCreate(r);
1287
1288  sBucket_pt& tail   = m_sum_bucket; assume( tail != NULL );
1289  kBucket_pt& bucket = m_spoly_bucket; assume( bucket != NULL );
1290  kbTest(bucket);
1291
1292
1293//  kBucketInit(bucket, NULL, 0); // not needed!?
1294
1295  poly p = leadmonom(syz_lead, r); // :(
1296//  poly spoly = pp_Mult_qq(p, T->m[c], r);
1297  kBucket_Plus_mm_Mult_pp(bucket, p, T->m[c], 0); // TODO: store pLength(T->m[c]) separately!?
1298  p_Delete(&p, r);
1299
1300  kbTest(bucket);
1301
1302  c = p_GetComp(syz_2, r) - 1;
1303  assume( c >= 0 && c < IDELEMS(T) );
1304
1305  p = leadmonom(syz_2, r); // :(
1306//  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
1307  kBucket_Plus_mm_Mult_pp(bucket, p, T->m[c], 0); // pLength(T->m[c])?!
1308  kbTest(bucket);
1309  p_Delete(&p, r);
1310
1311  // TODO: use bucket!?
1312//  const bool bUsePolynomial = TEST_OPT_NOT_BUCKETS; //  || (pLength(spoly) < MIN_LENGTH_BUCKET);
1313//  CPolynomialSummator tail(r, bUsePolynomial);
1314  sBucket_Add_p(tail, syz_2, 1); // tail.AddAndDelete(syz_2, 1);
1315
1316  kbTest(bucket);
1317  for( poly spoly = kBucketExtractLm(bucket); spoly != NULL; p_LmDelete(&spoly, r), spoly = kBucketExtractLm(bucket))
1318  {
1319    kbTest(bucket);
1320    poly t = m_div.FindReducer(spoly, NULL, m_checker);
1321
1322    if( t != NULL )
1323    {
1324      p = leadmonom(t, r); // :(
1325      c = p_GetComp(t, r) - 1;
1326
1327      assume( c >= 0 && c < IDELEMS(T) );
1328
1329      if( __TREEOUTPUT__ )
1330      {
1331        PrintS("{ \"nodelabel\": \""); writeLatexTerm(t, r); PrintS("\", \"edgelable\": \""); writeLatexTerm(spoly, r, false); PrintS("\" },");
1332      }
1333
1334      kBucket_Plus_mm_Mult_pp(bucket, p, T->m[c], 0); // pLength(T->m[c])?
1335//      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
1336
1337      p_Delete(&p, r);
1338
1339      sBucket_Add_p(tail, t, 1); // tail.AddAndDelete(t, 1);
1340    } // otherwise discard that leading term altogether!
1341    kbTest(bucket);
1342  }
1343
1344  // now bucket must be empty!
1345  kbTest(bucket);
1346  assume( kBucketClear(bucket) == NULL );
1347
1348  poly result; int len;
1349  sBucketClearAdd(tail, &result, &len); // TODO: use Merge with sBucket???
1350  assume( pLength(result) == len );
1351
1352  if( __TREEOUTPUT__ )
1353  {
1354    PrintS("]},");
1355  }
1356
1357
1358  return result;
1359}
1360
1361// namespace     {
1362
1363// };
1364
1365bool my_p_LmCmp (poly a, poly b, const ring r) { return p_LmCmp(a, b, r) == -1; } // TODO: change to simple lex. memory compare!
1366
1367// NOTE: need p_Copy?????? for image + multiplier!!???
1368// NOTE: better store complete syz. terms!!?
1369poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, const int tail) const
1370{
1371  const ring& r = m_rBaseRing;
1372
1373  assume(m_idTails !=  NULL && m_idTails->m != NULL);
1374  assume( tail >= 0 && tail < IDELEMS(m_idTails) );
1375
1376/*  return ComputeImage(multiplier, tail); */
1377
1378  // TODO: store (multiplier, tail) -.-^-.-^-.--> !
1379  TCache::iterator top_itr = m_cache.find(tail);
1380
1381  if ( top_itr != m_cache.end() )
1382  {
1383     assume( top_itr->first == tail );
1384
1385     TP2PCache& T = top_itr->second;
1386
1387     TP2PCache::iterator itr = T.find(multiplier);
1388
1389     if( itr != T.end() ) // Yey - Reuse!!!
1390     {
1391       assume( p_LmEqual(itr->first, multiplier, r) );
1392
1393       if( itr->second == NULL )
1394         return (NULL);
1395
1396       poly p = p_Copy(itr->second, r); // no copy???
1397
1398       if( __TREEOUTPUT__ )
1399       {
1400//         PrintS("{ \"nodelabel\": \""); writeLatexTerm(multiplier, r, false);
1401//         Print(" \\\\cdot \\\\GEN{%d}\", \"children\": [ ", tail + 1);
1402          Print("{ \"lookedupnode\": \"1\", \"nodelabel\": \"");
1403          writeLatexTerm(p, r, true, false);
1404       }
1405
1406       if( !n_Equal( pGetCoeff(multiplier), pGetCoeff(itr->first), r) ) // normalize coeffs!?
1407       {
1408         number n = n_Div( pGetCoeff(multiplier), pGetCoeff(itr->first), r);
1409         p = p_Mult_nn(p, n, r);
1410         n_Delete(&n, r);
1411
1412         if( __TREEOUTPUT__ )
1413           Print("\", \"RESCALEDRESULT\": \"1\" },");
1414       } else
1415       {
1416         if( __TREEOUTPUT__ )
1417           Print("\", \"RESCALEDRESULT\": \"0\" },");
1418       }
1419
1420       return p;
1421     }
1422
1423     if( __TREEOUTPUT__ )
1424     {
1425//        PrintS("{ \"nodelabel\": \""); writeLatexTerm(multiplier, r, false); Print(" \\\\cdot \\\\GEN{%d}\", \"children\": [", tail + 1);
1426     }
1427
1428     const poly p = ComputeImage(multiplier, tail);
1429     T.insert( TP2PCache::value_type(p_Copy(multiplier, r), p) ); //     T[ multiplier ] = p;
1430
1431//     if( p == NULL )
1432//        return (NULL);
1433
1434     if( __TREEOUTPUT__ )
1435     {
1436//        PrintS("], \"storedResult\": \""); writeLatexTerm(p, r, true, false); PrintS("\" },");
1437     }
1438
1439     return p_Copy(p, r);
1440  }
1441
1442  CCacheCompare o(r); TP2PCache T(o);
1443
1444  if( __TREEOUTPUT__ )
1445  {
1446//     PrintS("{ \"nodelabel\": \""); writeLatexTerm(multiplier, r, false); Print(" \\\\cdot \\\\GEN{%d}\", \"children\": [", tail + 1);
1447  }
1448
1449
1450  const poly p = ComputeImage(multiplier, tail);
1451
1452  T.insert( TP2PCache::value_type(p_Copy(multiplier, r), p) );
1453
1454  m_cache.insert( TCache::value_type(tail, T) );
1455
1456//  if( p == NULL )
1457//    return (NULL);
1458
1459  if( __TREEOUTPUT__ )
1460  {
1461//     PrintS("], \"storedResult\": \""); writeLatexTerm(p, r, true, false); PrintS("\" },");
1462  }
1463
1464  return p_Copy(p, r);
1465}
1466
1467poly SchreyerSyzygyComputation::ComputeImage(poly multiplier, const int tail) const
1468{
1469  const poly t = m_idTails->m[tail]; // !!!
1470
1471  if(t != NULL)
1472    return TraverseTail(multiplier, t);
1473
1474  return NULL;
1475}
1476
1477
1478poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
1479{
1480  assume( !__IGNORETAILS__ );
1481
1482  const ideal& L = m_idLeads;
1483  const ideal& T = m_idTails;
1484  const ring& r = m_rBaseRing;
1485
1486  assume( multiplier != NULL );
1487
1488  assume( L != NULL );
1489  assume( T != NULL );
1490
1491
1492  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
1493  {
1494//    const bool bUsePolynomial = TEST_OPT_NOT_BUCKETS; //  || (pLength(tail) < MIN_LENGTH_BUCKET);
1495    if( m_sum_bucket == NULL )
1496      m_sum_bucket = sBucketCreate(r);
1497
1498    sBucket_pt& sum   = m_sum_bucket; assume( sum != NULL );
1499    poly s; int len;
1500
1501//    CPolynomialSummator sum(r, bUsePolynomial);
1502//    poly s = NULL;
1503    for(poly p = tail; p != NULL; p = pNext(p))   // iterate over the tail
1504    {
1505      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1506      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1507      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1508      const poly rt = ReduceTerm(multiplier, p, NULL); // TODO: also return/store length?
1509      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1510      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1511      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1512
1513      const int lp = pLength(rt);
1514      if( rt != NULL && lp != 0 )
1515        sBucket_Add_p(sum, rt, lp);
1516    }
1517
1518    sBucketClearAdd(sum, &s, &len);
1519    assume( pLength(s) == len );
1520    return s;
1521  }
1522
1523  return NULL;
1524
1525}
1526
1527
1528
1529
1530poly SchreyerSyzygyComputation::ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck) const
1531{
1532  assume( !__IGNORETAILS__ );
1533
1534  const ideal& L = m_idLeads;
1535  const ideal& T = m_idTails;
1536  const ring& r = m_rBaseRing;
1537
1538  assume( multiplier != NULL );
1539  assume( term4reduction != NULL );
1540
1541
1542  assume( L != NULL );
1543  assume( T != NULL );
1544
1545  // simple implementation with FindReducer:
1546  poly s = NULL;
1547
1548  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
1549  {
1550#if NOPRODUCT
1551    s = m_div.FindReducer(multiplier, term4reduction, syztermCheck, m_checker);
1552
1553    if( s == NULL ) // No Reducer?
1554     return s;
1555
1556    if( __TREEOUTPUT__ )
1557    {
1558      PrintS("{ \"nodelabel\": \""); writeLatexTerm(s, r);
1559    }
1560
1561#else
1562    // NOTE: only LT(term4reduction) should be used in the following:
1563    poly product = pp_Mult_mm(multiplier, term4reduction, r);
1564    s = m_div.FindReducer(product, syztermCheck, m_checker);
1565
1566    if( s == NULL ) // No Reducer?
1567     return s;
1568
1569    if( __TREEOUTPUT__ )
1570    {
1571      PrintS("{ \"nodelabel\": \""); writeLatexTerm(s, r); PrintS("\", \"edgelable\": \""); writeLatexTerm(product, r, false);
1572    }
1573
1574    p_Delete(&product, r);
1575#endif
1576  }
1577
1578  if( s == NULL ) // No Reducer?
1579    return s;
1580
1581
1582  poly b = leadmonom(s, r);
1583
1584  const int c = p_GetComp(s, r) - 1;
1585  assume( c >= 0 && c < IDELEMS(T) );
1586
1587
1588  if( __TREEOUTPUT__ )
1589     PrintS("\", \"children\": [");
1590
1591  const poly t = TraverseTail(b, c); // T->m[c];
1592
1593  if( t != NULL )
1594    s = p_Add_q(s, t, r);
1595
1596  if( __TREEOUTPUT__ )
1597  {
1598
1599    PrintS("], \"noderesult\": \"");
1600    writeLatexTerm(s, r, true, false);
1601    PrintS("\" },"); 
1602  }
1603
1604  return s;
1605}
1606
1607SchreyerSyzygyComputationFlags::SchreyerSyzygyComputationFlags(idhdl rootRingHdl):
1608    __DEBUG__( atGetInt(rootRingHdl,"DEBUG", 0) ),
1609    __LEAD2SYZ__( atGetInt(rootRingHdl, "LEAD2SYZ", 0) ),
1610    __TAILREDSYZ__( atGetInt(rootRingHdl, "TAILREDSYZ", 1) ),
1611    __HYBRIDNF__( atGetInt(rootRingHdl, "HYBRIDNF", 1) ),
1612    __IGNORETAILS__( atGetInt(rootRingHdl, "IGNORETAILS", 0) ),
1613    __SYZNUMBER__( atGetInt(rootRingHdl, "SYZNUMBER", 0) ),
1614    __TREEOUTPUT__( atGetInt(rootRingHdl, "TREEOUTPUT", 0) ),
1615    __SYZCHECK__( atGetInt(rootRingHdl, "SYZCHECK", 0) ),
1616    m_rBaseRing( rootRingHdl->data.uring )
1617{
1618#ifndef SING_NDEBUG
1619  if( __DEBUG__ )
1620  {
1621    PrintS("SchreyerSyzygyComputationFlags: \n");
1622    Print("        DEBUG: \t%d\n", __DEBUG__);
1623//    Print("   SYZCHECK  : \t%d\n", __SYZCHECK__);
1624    Print("     LEAD2SYZ: \t%d\n", __LEAD2SYZ__);
1625    Print("   TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
1626    Print("  IGNORETAILS: \t%d\n", __IGNORETAILS__);
1627    Print("   TREEOUTPUT: \t%d\n", __TREEOUTPUT__);
1628    Print("     SYZCHECK: \t%d\n", __SYZCHECK__);
1629  }
1630#endif
1631
1632  // TODO: just current setting!
1633  assume( rootRingHdl == currRingHdl );
1634  assume( rootRingHdl->typ == RING_CMD );
1635  assume( m_rBaseRing == currRing );
1636  // move the global ring here inside???
1637}
1638
1639
1640
1641CLeadingTerm::CLeadingTerm(unsigned int _label,  const poly _lt, const ring R):
1642    m_sev( p_GetShortExpVector(_lt, R) ),  m_label( _label ),  m_lt( _lt )
1643{ }
1644
1645
1646CReducerFinder::~CReducerFinder()
1647{
1648  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++ )
1649  {
1650    const TReducers& v = it->second;
1651    for(TReducers::const_iterator vit = v.begin(); vit != v.end(); vit++ )
1652      delete const_cast<CLeadingTerm*>(*vit);
1653  }
1654}
1655
1656
1657void CReducerFinder::Initialize(const ideal L)
1658{
1659  assume( m_L == NULL || m_L == L );
1660  if( m_L == NULL )
1661    m_L = L;
1662
1663  assume( m_L == L );
1664
1665  if( L != NULL )
1666  {
1667    const ring& R = m_rBaseRing;
1668    assume( R != NULL );
1669
1670    for( int k = IDELEMS(L) - 1; k >= 0; k-- )
1671    {
1672      const poly a = L->m[k]; // assume( a != NULL );
1673
1674      // NOTE: label is k \in 0 ... |L|-1!!!
1675      if( a != NULL )
1676        m_hash[p_GetComp(a, R)].push_back( new CLeadingTerm(k, a, R) );
1677    }
1678  }
1679}
1680
1681CReducerFinder::CReducerFinder(const ideal L, const SchreyerSyzygyComputationFlags& flags):
1682    SchreyerSyzygyComputationFlags(flags),
1683    m_L(const_cast<ideal>(L)), // for debug anyway
1684    m_hash()
1685{
1686  assume( flags.m_rBaseRing == m_rBaseRing );
1687  if( L != NULL )
1688    Initialize(L);
1689}
1690
1691/// _p_LmDivisibleByNoComp for a | b*c
1692static inline BOOLEAN _p_LmDivisibleByNoComp(const poly a, const poly b, const poly c, const ring r)
1693{
1694  int i=r->VarL_Size - 1;
1695  unsigned long divmask = r->divmask;
1696  unsigned long la, lb;
1697
1698  if (r->VarL_LowIndex >= 0)
1699  {
1700    i += r->VarL_LowIndex;
1701    do
1702    {
1703      la = a->exp[i];
1704      lb = b->exp[i] + c->exp[i];
1705      if ((la > lb) ||
1706          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
1707      {
1708        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
1709        return FALSE;
1710      }
1711      i--;
1712    }
1713    while (i>=r->VarL_LowIndex);
1714  }
1715  else
1716  {
1717    do
1718    {
1719      la = a->exp[r->VarL_Offset[i]];
1720      lb = b->exp[r->VarL_Offset[i]] + c->exp[r->VarL_Offset[i]];
1721      if ((la > lb) ||
1722          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
1723      {
1724        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
1725        return FALSE;
1726      }
1727      i--;
1728    }
1729    while (i>=0);
1730  }
1731#ifdef HAVE_RINGS
1732  assume( !rField_is_Ring(r) ); // not implemented for rings...!
1733#endif
1734  return TRUE;
1735}
1736
1737bool CLeadingTerm::DivisibilityCheck(const poly product, const unsigned long not_sev, const ring r) const
1738{
1739  const poly p = m_lt;
1740
1741  assume( p_GetComp(p, r) == p_GetComp(product, r) );
1742
1743//  const int k = m_label;
1744//  assume( m_L->m[k] == p );
1745
1746  const unsigned long p_sev = m_sev;
1747
1748  assume( p_sev == p_GetShortExpVector(p, r) );
1749
1750  return p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r);
1751
1752}
1753
1754/// as DivisibilityCheck(multiplier * t, ...) for monomial 'm'
1755/// and a module term 't'
1756bool CLeadingTerm::DivisibilityCheck(const poly m, const poly t, const unsigned long not_sev, const ring r) const
1757{
1758  const poly p = m_lt;
1759
1760  assume( p_GetComp(p, r) == p_GetComp(t, r) );
1761// assume( p_GetComp(m, r) == 0 );
1762
1763//  const int k = m_label;
1764//  assume( m_L->m[k] == p );
1765
1766  const unsigned long p_sev = m_sev;
1767  assume( p_sev == p_GetShortExpVector(p, r) );
1768
1769  if (p_sev & not_sev)
1770    return false;
1771
1772  return _p_LmDivisibleByNoComp(p, m, t, r);
1773
1774//  return p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r);
1775
1776}
1777
1778
1779
1780/// TODO:
1781class CDivisorEnumerator: public SchreyerSyzygyComputationFlags
1782{
1783  private:
1784    const CReducerFinder& m_reds;
1785    const poly m_product;
1786    const unsigned long m_not_sev;
1787    const long m_comp;
1788
1789    CReducerFinder::CReducersHash::const_iterator m_itr;
1790    CReducerFinder::TReducers::const_iterator m_current, m_finish;
1791
1792    bool m_active;
1793
1794  public:
1795    CDivisorEnumerator(const CReducerFinder& self, const poly product):
1796        SchreyerSyzygyComputationFlags(self),
1797        m_reds(self),
1798        m_product(product),
1799        m_not_sev(~p_GetShortExpVector(product, m_rBaseRing)),
1800        m_comp(p_GetComp(product, m_rBaseRing)),
1801        m_itr(), m_current(), m_finish(),
1802        m_active(false)
1803    {
1804      assume( m_comp >= 0 );
1805      assume( m_reds.m_L != NULL );
1806    }
1807
1808    inline bool Reset()
1809    {
1810      m_active = false;
1811
1812      m_itr = m_reds.m_hash.find(m_comp);
1813
1814      if( m_itr == m_reds.m_hash.end() )
1815        return false;
1816
1817      assume( m_itr->first == m_comp );
1818
1819      m_current = (m_itr->second).begin();
1820      m_finish = (m_itr->second).end();
1821
1822      if (m_current == m_finish)
1823        return false;
1824
1825//      m_active = true;
1826      return true;
1827    }
1828
1829    const CLeadingTerm& Current() const
1830    {
1831      assume( m_active );
1832      assume( m_current != m_finish );
1833
1834      return *(*m_current);
1835    }
1836
1837    inline bool MoveNext()
1838    {
1839      assume( m_current != m_finish );
1840
1841      if( m_active )
1842        ++m_current;
1843      else
1844        m_active = true; // for Current()
1845
1846      // looking for the next good entry
1847      for( ; m_current != m_finish; ++m_current )
1848      {
1849        assume( m_reds.m_L->m[Current().m_label] == Current().m_lt );
1850
1851        if( Current().DivisibilityCheck(m_product, m_not_sev, m_rBaseRing) )
1852        {
1853#ifndef SING_NDEBUG
1854          if( __DEBUG__ )
1855          {
1856            Print("CDivisorEnumerator::MoveNext::est LS: q is divisible by LS[%d] !:((, diviser is: ", 1 + Current().m_label);
1857            dPrint(Current().m_lt, m_rBaseRing, m_rBaseRing, 1);
1858          }
1859#endif
1860//          m_active = true;
1861          return true;
1862        }
1863      }
1864
1865      // the end... :(
1866      assume( m_current == m_finish );
1867
1868      m_active = false;
1869      return false;
1870    }
1871};
1872
1873
1874
1875bool CReducerFinder::IsDivisible(const poly product) const
1876{
1877  CDivisorEnumerator itr(*this, product);
1878  if( !itr.Reset() )
1879    return false;
1880
1881  return itr.MoveNext();
1882
1883/*
1884  const ring& r = m_rBaseRing;
1885
1886  const long comp = p_GetComp(product, r);
1887  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
1888
1889  assume( comp >= 0 );
1890
1891  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
1892
1893  assume( m_L != NULL );
1894
1895  if( it == m_hash.end() )
1896    return false;
1897  // assume comp!
1898
1899  const TReducers& reducers = it->second;
1900
1901  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1902  {
1903    assume( m_L->m[(*vit)->m_label] == (*vit)->m_lt );
1904
1905    if( (*vit)->DivisibilityCheck(product, not_sev, r) )
1906    {
1907      if( __DEBUG__ )
1908      {
1909        Print("_FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", 1 + (*vit)->m_label);
1910        dPrint((*vit)->m_lt, r, r, 1);
1911      }
1912
1913      return true;
1914    }
1915  }
1916
1917  return false;
1918*/
1919}
1920
1921
1922#ifndef SING_NDEBUG
1923void CReducerFinder::DebugPrint() const
1924{
1925  const ring& r = m_rBaseRing;
1926
1927  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++)
1928  {
1929    Print("Hash Key: %ld, Values: \n", it->first);
1930    const TReducers& reducers = it->second;
1931
1932    for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1933    {
1934      const poly p = (*vit)->m_lt;
1935
1936      assume( p_GetComp(p, r) == it->first );
1937
1938      const int k = (*vit)->m_label;
1939
1940      assume( m_L->m[k] == p );
1941
1942      const unsigned long p_sev = (*vit)->m_sev;
1943
1944      assume( p_sev == p_GetShortExpVector(p, r) );
1945
1946      Print("L[%d]: ", k); dPrint(p, r, r, 0); Print("SEV: %ld\n", p_sev);
1947    }
1948  }
1949}
1950#endif
1951
1952/// TODO:
1953class CDivisorEnumerator2: public SchreyerSyzygyComputationFlags
1954{
1955  private:
1956    const CReducerFinder& m_reds;
1957    const poly m_multiplier, m_term;
1958    const unsigned long m_not_sev;
1959    const long m_comp;
1960
1961    CReducerFinder::CReducersHash::const_iterator m_itr;
1962    CReducerFinder::TReducers::const_iterator m_current, m_finish;
1963
1964    bool m_active;
1965
1966  public:
1967    CDivisorEnumerator2(const CReducerFinder& self, const poly m, const poly t):
1968        SchreyerSyzygyComputationFlags(self),
1969        m_reds(self),
1970        m_multiplier(m), m_term(t),
1971        m_not_sev(~p_GetShortExpVector(m, t, m_rBaseRing)),
1972        m_comp(p_GetComp(t, m_rBaseRing)),
1973        m_itr(), m_current(), m_finish(),
1974        m_active(false)
1975    {
1976      assume( m_comp >= 0 );
1977      assume( m_reds.m_L != NULL );
1978      assume( m_multiplier != NULL );
1979      assume( m_term != NULL );
1980//      assume( p_GetComp(m_multiplier, m_rBaseRing) == 0 );
1981    }
1982
1983    inline bool Reset()
1984    {
1985      m_active = false;
1986
1987      m_itr = m_reds.m_hash.find(m_comp);
1988
1989      if( m_itr == m_reds.m_hash.end() )
1990        return false;
1991
1992      assume( m_itr->first == m_comp );
1993
1994      m_current = (m_itr->second).begin();
1995      m_finish = (m_itr->second).end();
1996
1997      if (m_current == m_finish)
1998        return false;
1999
2000      return true;
2001    }
2002
2003    const CLeadingTerm& Current() const
2004    {
2005      assume( m_active );
2006      assume( m_current != m_finish );
2007
2008      return *(*m_current);
2009    }
2010
2011    inline bool MoveNext()
2012    {
2013      assume( m_current != m_finish );
2014
2015      if( m_active )
2016        ++m_current;
2017      else
2018        m_active = true;
2019
2020
2021      // looking for the next good entry
2022      for( ; m_current != m_finish; ++m_current )
2023      {
2024        assume( m_reds.m_L->m[Current().m_label] == Current().m_lt );
2025
2026        if( Current().DivisibilityCheck(m_multiplier, m_term, m_not_sev, m_rBaseRing) )
2027        {
2028#ifndef SING_NDEBUG
2029          if( __DEBUG__ )
2030          {
2031            Print("CDivisorEnumerator::MoveNext::est LS: q is divisible by LS[%d] !:((, diviser is: ", 1 + Current().m_label);
2032            dPrint(Current().m_lt, m_rBaseRing, m_rBaseRing, 1);
2033          }
2034#endif
2035//          m_active = true;
2036          return true;
2037
2038        }
2039      }
2040
2041      // the end... :(
2042      assume( m_current == m_finish );
2043
2044      m_active = false;
2045      return false;
2046    }
2047};
2048
2049poly CReducerFinder::FindReducer(const poly multiplier, const poly t,
2050                                 const poly syzterm,
2051                                 const CReducerFinder& syz_checker) const
2052{
2053  CDivisorEnumerator2 itr(*this, multiplier, t);
2054  if( !itr.Reset() )
2055    return NULL;
2056
2057  // don't care about the module component of multiplier (as it may be the syzygy term)
2058  // product = multiplier * t?
2059  const ring& r = m_rBaseRing;
2060
2061  assume( multiplier != NULL ); assume( t != NULL );
2062
2063  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
2064
2065  long c = 0;
2066
2067  if (syzterm != NULL)
2068    c = p_GetComp(syzterm, r) - 1;
2069
2070  assume( c >= 0 && c < IDELEMS(L) );
2071
2072  if (__DEBUG__ && (syzterm != NULL))
2073  {
2074    const poly m = L->m[c];
2075
2076    assume( m != NULL ); assume( pNext(m) == NULL );
2077
2078    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
2079
2080    poly pr = p_Mult_q( leadmonom(multiplier, r, false), leadmonom(t, r, false), r);
2081
2082    assume( p_EqualPolys(lm, pr, r) );
2083
2084    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
2085    //  def @@product = leadmonomial(syzterm) * L[@@r];
2086
2087    p_Delete(&lm, r);
2088    p_Delete(&pr, r);
2089  }
2090
2091  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
2092
2093  const poly q = p_New(r); pNext(q) = NULL;
2094
2095  if( __DEBUG__ )
2096    p_SetCoeff0(q, 0, r); // for printing q
2097
2098  while( itr.MoveNext() )
2099  {
2100    const poly p = itr.Current().m_lt;
2101    const int k  = itr.Current().m_label;
2102
2103    p_ExpVectorSum(q, multiplier, t, r); // q == product == multiplier * t // TODO: do it once?
2104    p_ExpVectorDiff(q, q, p, r); // (LM(product) / LM(L[k]))
2105
2106    p_SetComp(q, k + 1, r);
2107    p_Setm(q, r);
2108
2109    // cannot allow something like: a*gen(i) - a*gen(i)
2110    if (syzterm != NULL && (k == c))
2111      if (p_ExpVectorEqual(syzterm, q, r))
2112      {
2113#ifndef SING_NDEBUG
2114        if( __DEBUG__ )
2115        {
2116          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
2117          dPrint(syzterm, r, r, 1);
2118        }
2119#endif
2120        continue;
2121      }
2122
2123    // while the complement (the fraction) is not reducible by leading syzygies
2124    if( to_check && syz_checker.IsDivisible(q) )
2125    {
2126#ifndef SING_NDEBUG
2127      if( __DEBUG__ )
2128      {
2129        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
2130      }
2131#endif
2132      continue;
2133    }
2134
2135    number n = n_Mult( p_GetCoeff(multiplier, r), p_GetCoeff(t, r), r);
2136    p_SetCoeff0(q, n_InpNeg( n_Div(n, p_GetCoeff(p, r), r), r), r);
2137    n_Delete(&n, r);
2138
2139    return q;
2140  }
2141
2142/*
2143  const long comp = p_GetComp(t, r); assume( comp >= 0 );
2144  const unsigned long not_sev = ~p_GetShortExpVector(multiplier, t, r); // !
2145
2146//   for( int k = IDELEMS(L)-1; k>= 0; k-- )
2147//   {
2148//     const poly p = L->m[k];
2149//
2150//     if ( p_GetComp(p, r) != comp )
2151//       continue;
2152//
2153//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
2154
2155   // looking for an appropriate diviser p = L[k]...
2156  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
2157
2158  if( it == m_hash.end() )
2159    return NULL;
2160
2161  // assume comp!
2162
2163  assume( m_L != NULL );
2164
2165  const TReducers& reducers = it->second;
2166
2167  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
2168  {
2169
2170    const poly p = (*vit)->m_lt;
2171    const int k = (*vit)->m_label;
2172
2173    assume( L->m[k] == p );
2174
2175//    const unsigned long p_sev = (*vit)->m_sev;
2176//    assume( p_sev == p_GetShortExpVector(p, r) );
2177
2178//    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
2179//      continue;
2180
2181    if( !(*vit)->DivisibilityCheck(multiplier, t, not_sev, r) )
2182      continue;
2183
2184
2185//    if (p_sev & not_sev) continue;
2186//    if( !_p_LmDivisibleByNoComp(p, multiplier, t, r) ) continue;
2187
2188
2189    p_ExpVectorSum(q, multiplier, t, r); // q == product == multiplier * t
2190    p_ExpVectorDiff(q, q, p, r); // (LM(product) / LM(L[k]))
2191
2192    p_SetComp(q, k + 1, r);
2193    p_Setm(q, r);
2194
2195    // cannot allow something like: a*gen(i) - a*gen(i)
2196    if (syzterm != NULL && (k == c))
2197      if (p_ExpVectorEqual(syzterm, q, r))
2198      {
2199        if( __DEBUG__ )
2200        {
2201          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
2202          dPrint(syzterm, r, r, 1);
2203        }
2204
2205        continue;
2206      }
2207
2208    // while the complement (the fraction) is not reducible by leading syzygies
2209    if( to_check && syz_checker.IsDivisible(q) )
2210    {
2211      if( __DEBUG__ )
2212      {
2213        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
2214      }
2215
2216      continue;
2217    }
2218
2219    number n = n_Mult( p_GetCoeff(multiplier, r), p_GetCoeff(t, r), r);
2220    p_SetCoeff0(q, n_InpNeg( n_Div(n, p_GetCoeff(p, r), r), r), r);
2221    n_Delete(&n, r);
2222
2223    return q;
2224  }
2225*/
2226
2227  p_LmFree(q, r);
2228
2229  return NULL;
2230
2231}
2232
2233
2234poly CReducerFinder::FindReducer(const poly product, const poly syzterm, const CReducerFinder& syz_checker) const
2235{
2236  CDivisorEnumerator itr(*this, product);
2237  if( !itr.Reset() )
2238    return NULL;
2239
2240
2241  const ring& r = m_rBaseRing;
2242
2243  assume( product != NULL );
2244
2245  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
2246
2247  long c = 0;
2248
2249  if (syzterm != NULL)
2250    c = p_GetComp(syzterm, r) - 1;
2251
2252  assume( c >= 0 && c < IDELEMS(L) );
2253
2254  if (__DEBUG__ && (syzterm != NULL))
2255  {
2256    const poly m = L->m[c];
2257
2258    assume( m != NULL ); assume( pNext(m) == NULL );
2259
2260    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
2261    assume( p_EqualPolys(lm, product, r) );
2262
2263    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
2264    //  def @@product = leadmonomial(syzterm) * L[@@r];
2265
2266    p_Delete(&lm, r);
2267  }
2268
2269
2270  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
2271
2272  const poly q = p_New(r); pNext(q) = NULL;
2273
2274  if( __DEBUG__ )
2275    p_SetCoeff0(q, 0, r); // for printing q
2276
2277  while( itr.MoveNext() )
2278  {
2279    const poly p = itr.Current().m_lt;
2280    const int k  = itr.Current().m_label;
2281
2282    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
2283    p_SetComp(q, k + 1, r);
2284    p_Setm(q, r);
2285
2286    // cannot allow something like: a*gen(i) - a*gen(i)
2287    if (syzterm != NULL && (k == c))
2288      if (p_ExpVectorEqual(syzterm, q, r))
2289      {
2290#ifndef SING_NDEBUG
2291        if( __DEBUG__ )
2292        {
2293          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
2294          dPrint(syzterm, r, r, 1);
2295        }
2296#endif
2297        continue;
2298      }
2299
2300    // while the complement (the fraction) is not reducible by leading syzygies
2301    if( to_check && syz_checker.IsDivisible(q) )
2302    {
2303#ifndef SING_NDEBUG
2304      if( __DEBUG__ )
2305      {
2306        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
2307      }
2308#endif
2309      continue;
2310    }
2311
2312    p_SetCoeff0(q, n_InpNeg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
2313
2314    return q;
2315  }
2316
2317
2318
2319/*
2320  const long comp = p_GetComp(product, r);
2321  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
2322
2323  assume( comp >= 0 );
2324
2325//   for( int k = IDELEMS(L)-1; k>= 0; k-- )
2326//   {
2327//     const poly p = L->m[k];
2328//
2329//     if ( p_GetComp(p, r) != comp )
2330//       continue;
2331//
2332//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
2333
2334   // looking for an appropriate diviser p = L[k]...
2335  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
2336
2337  if( it == m_hash.end() )
2338    return NULL;
2339
2340  assume( m_L != NULL );
2341
2342  const TReducers& reducers = it->second;
2343
2344  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
2345
2346  const poly q = p_New(r); pNext(q) = NULL;
2347
2348  if( __DEBUG__ )
2349    p_SetCoeff0(q, 0, r); // for printing q
2350
2351  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
2352  {
2353    const poly p = (*vit)->m_lt;
2354
2355    assume( p_GetComp(p, r) == comp );
2356
2357    const int k = (*vit)->m_label;
2358
2359    assume( L->m[k] == p );
2360
2361    const unsigned long p_sev = (*vit)->m_sev;
2362
2363    assume( p_sev == p_GetShortExpVector(p, r) );
2364
2365    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
2366      continue;
2367
2368//     // ... which divides the product, looking for the _1st_ appropriate one!
2369//     if( !p_LmDivisibleByNoComp(p, product, r) ) // included inside  p_LmShortDivisibleBy!
2370//       continue;
2371
2372    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
2373    p_SetComp(q, k + 1, r);
2374    p_Setm(q, r);
2375
2376    // cannot allow something like: a*gen(i) - a*gen(i)
2377    if (syzterm != NULL && (k == c))
2378      if (p_ExpVectorEqual(syzterm, q, r))
2379      {
2380        if( __DEBUG__ )
2381        {
2382          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
2383          dPrint(syzterm, r, r, 1);
2384        }
2385
2386        continue;
2387      }
2388
2389    // while the complement (the fraction) is not reducible by leading syzygies
2390    if( to_check && syz_checker.IsDivisible(q) )
2391    {
2392      if( __DEBUG__ )
2393      {
2394        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
2395      }
2396
2397      continue;
2398    }
2399
2400    p_SetCoeff0(q, n_InpNeg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
2401    return q;
2402  }
2403*/
2404
2405  p_LmFree(q, r);
2406
2407  return NULL;
2408}
2409
2410
2411
2412CLCM::CLCM(const ideal& L, const SchreyerSyzygyComputationFlags& flags):
2413    SchreyerSyzygyComputationFlags(flags), std::vector<bool>(),
2414    m_compute(false), m_N(rVar(flags.m_rBaseRing))
2415{
2416  const ring& R = m_rBaseRing;
2417  assume( flags.m_rBaseRing == R );
2418  assume( R != NULL );
2419
2420  assume( L != NULL );
2421
2422  if( __TAILREDSYZ__ && !__HYBRIDNF__ && (L != NULL))
2423  {
2424    const int l = IDELEMS(L);
2425
2426    assume( l > 0 );
2427
2428    resize(l, false);
2429
2430    for( int k = l - 1; k >= 0; k-- )
2431    {
2432      const poly a = L->m[k]; assume( a != NULL );
2433
2434      for (unsigned int j = m_N; j > 0; j--)
2435        if ( !(*this)[j] )
2436          (*this)[j] = (p_GetExp(a, j, R) > 0);
2437    }
2438
2439    m_compute = true;
2440  }
2441}
2442
2443
2444bool CLCM::Check(const poly m) const
2445{
2446  assume( m != NULL );
2447  if( m_compute && (m != NULL))
2448  {
2449    const ring& R = m_rBaseRing;
2450
2451    assume( __TAILREDSYZ__ && !__HYBRIDNF__ );
2452
2453    for (unsigned int j = m_N; j > 0; j--)
2454      if ( (*this)[j] )
2455        if(p_GetExp(m, j, R) > 0)
2456          return true;
2457
2458    return false;
2459
2460  } else return true;
2461}
2462
2463CCacheCompare::CCacheCompare(): m_ring(currRing) {}
2464
2465
2466template class std::vector<bool>;
2467template class std::vector<CLeadingTerm const*>;
2468template class std::map< CReducerFinder::TComponentKey, CReducerFinder::TReducers >;
2469
2470template class std::map<TCacheKey, TCacheValue, struct CCacheCompare>;
2471template class std::map<int, TP2PCache>;
2472
2473
2474
2475END_NAMESPACE               END_NAMESPACE_SINGULARXX
2476
2477
2478// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
Note: See TracBrowser for help on using the repository browser.