source: git/dyn_modules/syzextra/syzextra.cc @ 31a08c2

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