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

jengelh-datetimespielwiese
Last change on this file since dd24e5 was dd24e5, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
separated FindReducer into a separate class CReducerFinder (m_div) TODO: deal with m_LS via another CReducerFinder? NOTE: this change resulted in some speedups esp. for Hybrid approach (see schreyer.lib)
  • Property mode set to 100644
File size: 26.4 KB
Line 
1// -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2/*****************************************************************************\
3 * Computer Algebra System SINGULAR   
4\*****************************************************************************/
5/** @file syzextra.cc
6 *
7 * Here we implement the Computation of Syzygies
8 *
9 * ABSTRACT: Computation of Syzygies due to Schreyer
10 *
11 * @author Oleksandr Motsak
12 *
13 **/
14/*****************************************************************************/
15
16// include header file
17#include <kernel/mod2.h>
18
19#include "syzextra.h"
20
21#include "DebugPrint.h"
22
23#include <omalloc/omalloc.h>
24
25#include <misc/intvec.h>
26#include <misc/options.h>
27
28#include <coeffs/coeffs.h>
29
30#include <polys/monomials/p_polys.h>
31#include <polys/monomials/ring.h>
32
33#include <kernel/kstd1.h>
34#include <kernel/polys.h>
35#include <kernel/syz.h>
36#include <kernel/ideals.h>
37
38
39
40#include <Singular/tok.h>
41#include <Singular/ipid.h>
42#include <Singular/lists.h>
43#include <Singular/attrib.h>
44
45#include <Singular/ipid.h> 
46#include <Singular/ipshell.h> // For iiAddCproc
47
48#include <stdio.h>
49#include <stdlib.h>
50#include <string.h>
51
52// USING_NAMESPACE_SINGULARXX;
53USING_NAMESPACE( SINGULARXXNAME :: DEBUG )
54
55
56BEGIN_NAMESPACE_SINGULARXX     BEGIN_NAMESPACE(SYZEXTRA)
57
58
59BEGIN_NAMESPACE(SORT_c_ds)
60
61
62#ifdef _GNU_SOURCE
63static int cmp_c_ds(const void *p1, const void *p2, void *R)
64{
65#else
66static int cmp_c_ds(const void *p1, const void *p2)
67{
68  void *R = currRing; 
69#endif
70
71  const int YES = 1;
72  const int NO = -1;
73
74  const ring r =  (const ring) R; // TODO/NOTE: the structure is known: C, lp!!!
75
76  assume( r == currRing );
77
78  const poly a = *(const poly*)p1;
79  const poly b = *(const poly*)p2;
80
81  assume( a != NULL );
82  assume( b != NULL );
83
84  assume( p_LmTest(a, r) );
85  assume( p_LmTest(b, r) );
86
87
88  const signed long iCompDiff = p_GetComp(a, r) - p_GetComp(b, r);
89
90  // TODO: test this!!!!!!!!!!!!!!!!
91
92  //return -( compare (c, qsorts) )
93
94#ifndef NDEBUG
95  const int __DEBUG__ = 0;
96  if( __DEBUG__ )
97  {
98    PrintS("cmp_c_ds: a, b: \np1: "); dPrint(a, r, r, 2);
99    PrintS("b: "); dPrint(b, r, r, 2);
100    PrintLn();
101  }
102#endif
103
104
105  if( iCompDiff > 0 )
106    return YES;
107
108  if( iCompDiff < 0 )
109    return  NO;
110
111  assume( iCompDiff == 0 );
112
113  const signed long iDegDiff = p_Totaldegree(a, r) - p_Totaldegree(b, r);
114
115  if( iDegDiff > 0 )
116    return YES;
117
118  if( iDegDiff < 0 )
119    return  NO;
120
121  assume( iDegDiff == 0 );
122
123#ifndef NDEBUG
124  if( __DEBUG__ )
125  {
126    PrintS("cmp_c_ds: a & b have the same comp & deg! "); PrintLn();
127  }
128#endif
129
130  for (int v = rVar(r); v > 0; v--)
131  {
132    assume( v > 0 );
133    assume( v <= rVar(r) );
134
135    const signed int d = p_GetExp(a, v, r) - p_GetExp(b, v, r);
136
137    if( d > 0 )
138      return YES;
139
140    if( d < 0 )
141      return NO;
142
143    assume( d == 0 );
144  }
145
146  return 0; 
147}
148
149END_NAMESPACE
150/* namespace SORT_c_ds */
151
152/// return a new term: leading coeff * leading monomial of p
153/// with 0 leading component!
154poly leadmonom(const poly p, const ring r)
155{
156  poly m = NULL;
157
158  if( p != NULL )
159  {
160    assume( p != NULL );
161    assume( p_LmTest(p, r) );
162
163    m = p_LmInit(p, r);
164    p_SetCoeff0(m, n_Copy(p_GetCoeff(p, r), r), r);
165
166    p_SetComp(m, 0, r);
167    p_Setm(m, r);
168
169    assume( p_GetComp(m, r) == 0 );
170    assume( m != NULL );
171    assume( pNext(m) == NULL );
172    assume( p_LmTest(m, r) );
173  }
174
175  return m;
176}
177
178
179   
180poly p_Tail(const poly p, const ring r)
181{
182  if( p == NULL)
183    return NULL;
184  else
185    return p_Copy( pNext(p), r );
186}
187
188
189ideal id_Tail(const ideal id, const ring r)
190{
191  if( id == NULL)
192    return NULL;
193
194  const ideal newid = idInit(IDELEMS(id),id->rank);
195 
196  for (int i=IDELEMS(id) - 1; i >= 0; i--)
197    newid->m[i] = p_Tail( id->m[i], r );
198
199  newid->rank = id_RankFreeModule(newid, currRing);
200
201  return newid; 
202}
203
204
205
206void Sort_c_ds(const ideal id, const ring r)
207{
208  const int sizeNew = IDELEMS(id);
209
210#ifdef _GNU_SOURCE
211#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp, r)
212#else
213#define qsort_my(m, s, ss, r, cmp) qsort_r(m, s, ss, cmp)
214#endif
215 
216  if( sizeNew >= 2 )
217    qsort_my(id->m, sizeNew, sizeof(poly), r, FROM_NAMESPACE(SORT_c_ds, cmp_c_ds));
218
219#undef qsort_my
220
221  id->rank = id_RankFreeModule(id, r);
222}
223
224
225ideal SchreyerSyzygyComputation::Compute1LeadingSyzygyTerms()
226{
227  const ideal& id = m_idLeads;
228  const ring& r = m_rBaseRing;
229  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
230
231
232  //   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
233//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
234  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
235//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
236//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
237
238  assume(!__LEAD2SYZ__);
239
240  // 1. set of components S?
241  // 2. for each component c from S: set of indices of leading terms
242  // with this component?
243  // 3. short exp. vectors for each leading term?
244
245  const int size = IDELEMS(id);
246
247  if( size < 2 )
248  {
249    const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal...       
250    return newid;
251  }
252
253  // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!??
254
255  // components should come in groups: count elements in each group
256  // && estimate the real size!!!
257
258
259  // use just a vector instead???
260  const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
261
262  int k = 0;
263
264  for (int j = 0; j < size; j++)
265  {
266    const poly p = id->m[j];
267    assume( p != NULL );
268    const int  c = p_GetComp(p, r);
269
270    for (int i = j - 1; i >= 0; i--)
271    {
272      const poly pp = id->m[i];
273      assume( pp != NULL );
274      const int  cc = p_GetComp(pp, r);
275
276      if( c != cc )
277        continue;
278
279      const poly m = p_Init(r); // p_New???
280
281      // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)!
282      for (int v = rVar(r); v > 0; v--)
283      {
284        assume( v > 0 );
285        assume( v <= rVar(r) );
286
287        const short e1 = p_GetExp(p , v, r);
288        const short e2 = p_GetExp(pp, v, r);
289
290        if( e1 >= e2 )
291          p_SetExp(m, v, 0, r);
292        else
293          p_SetExp(m, v, e2 - e1, r);
294
295      }
296
297      assume( (j > i) && (i >= 0) );
298
299      p_SetComp(m, j + 1, r);
300      pNext(m) = NULL;
301      p_SetCoeff0(m, n_Init(1, r->cf), r); // for later...
302
303      p_Setm(m, r); // should not do anything!!!
304
305      newid->m[k++] = m;
306    }
307  }
308
309//   if( __DEBUG__ && FALSE )
310//   {
311//     PrintS("ComputeLeadingSyzygyTerms::Temp0: \n");
312//     dPrint(newid, r, r, 1);
313//   }
314
315  // the rest of newid is assumed to be zeroes...
316
317  // simplify(newid, 2 + 32)??
318  // sort(newid, "C,ds")[1]???
319  id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
320
321//   if( __DEBUG__ && FALSE )
322//   {
323//     PrintS("ComputeLeadingSyzygyTerms::Temp1: \n");
324//     dPrint(newid, r, r, 1);
325//   }
326
327  idSkipZeroes(newid); // #define SIMPL_NULL 2
328
329//   if( __DEBUG__ )
330//   {
331//     PrintS("ComputeLeadingSyzygyTerms::Output: \n");
332//     dPrint(newid, r, r, 1);
333//   }
334
335  Sort_c_ds(newid, r);
336
337  return newid;
338}
339
340ideal SchreyerSyzygyComputation::Compute2LeadingSyzygyTerms()
341{
342  const ideal& id = m_idLeads;
343  const ring& r = m_rBaseRing;
344  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
345 
346//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
347//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
348//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
349//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
350  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
351 
352
353  // 1. set of components S?
354  // 2. for each component c from S: set of indices of leading terms
355  // with this component?
356  // 3. short exp. vectors for each leading term?
357
358  const int size = IDELEMS(id);
359
360  if( size < 2 )
361  {
362    const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...   
363    return newid;
364  }
365
366
367  // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
368 
369  // components should come in groups: count elements in each group
370  // && estimate the real size!!!
371
372
373  // use just a vector instead???
374  ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
375
376  int k = 0;
377
378  for (int j = 0; j < size; j++)
379  {
380    const poly p = id->m[j];
381    assume( p != NULL );
382    const int  c = p_GetComp(p, r);
383
384    for (int i = j - 1; i >= 0; i--)
385    {
386      const poly pp = id->m[i];
387      assume( pp != NULL );
388      const int  cc = p_GetComp(pp, r);
389
390      if( c != cc )
391        continue;
392
393        // allocate memory & zero it out!
394      const poly m = p_Init(r); const poly mm = p_Init(r);
395
396
397        // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
398        // TODO: optimize: knowing the ring structure: (C/lp)!
399
400      for (int v = rVar(r); v > 0; v--)
401      {
402        assume( v > 0 );
403        assume( v <= rVar(r) );
404
405        const short e1 = p_GetExp(p , v, r);
406        const short e2 = p_GetExp(pp, v, r);
407
408        if( e1 >= e2 )
409          p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
410        else
411          p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
412
413      }
414
415      assume( (j > i) && (i >= 0) );
416
417      p_SetComp(m, j + 1, r);
418      p_SetComp(mm, i + 1, r);
419
420      const number& lc1 = p_GetCoeff(p , r);
421      const number& lc2 = p_GetCoeff(pp, r);
422
423      number g = n_Lcm( lc1, lc2, r );
424
425      p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
426      p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);
427
428      n_Delete(&g, r);
429
430      p_Setm(m, r); // should not do anything!!!
431      p_Setm(mm, r); // should not do anything!!!
432
433      pNext(m) = mm; //        pNext(mm) = NULL;
434
435      newid->m[k++] = m;
436    }
437  }
438
439//   if( __DEBUG__ && FALSE )
440//   {
441//     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
442//     dPrint(newid, r, r, 1);
443//   }
444
445  if( !__TAILREDSYZ__ )
446  {
447      // simplify(newid, 2 + 32)??
448      // sort(newid, "C,ds")[1]???
449    id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
450
451//     if( __DEBUG__ && FALSE )
452//     {
453//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
454//       dPrint(newid, r, r, 1);
455//     }
456  }
457  else
458  {
459      //      option(redSB); option(redTail);
460      //      TEST_OPT_REDSB
461      //      TEST_OPT_REDTAIL
462    assume( r == currRing );
463    BITSET _save_test = test;
464    test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB));
465
466    intvec* w=new intvec(IDELEMS(newid));
467    ideal tmp = kStd(newid, currQuotient, isHomog, &w);
468    delete w;
469
470    test = _save_test;
471
472    id_Delete(&newid, r);
473    newid = tmp;
474
475//     if( __DEBUG__ && FALSE )
476//     {
477//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
478//       dPrint(newid, r, r, 1);
479//     }
480
481  }
482
483  idSkipZeroes(newid);
484
485  Sort_c_ds(newid, r);
486 
487  return newid;
488}
489
490
491CReducerFinder::CLeadingTerm::CLeadingTerm(unsigned int _label,  const poly _lt, const ring R):
492    m_sev( p_GetShortExpVector(_lt, R) ),  m_label( _label ),  m_lt( _lt )
493    { }
494
495
496CReducerFinder::~CReducerFinder()
497{
498  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++ )
499  {
500    const TReducers& v = it->second;
501    for(TReducers::const_iterator vit = v.begin(); vit != v.end(); vit++ )
502      delete const_cast<CLeadingTerm*>(*vit);
503  }
504}
505                                   
506CReducerFinder::CReducerFinder(const SchreyerSyzygyComputation& data): m_data(data), m_hash()
507{
508
509  const ideal& L = data.m_idLeads;
510  const ring&  R = data.m_rBaseRing;
511//   const SchreyerSyzygyComputationFlags& attributes = data.m_atttributes;
512//
513//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
514//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
515//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
516//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
517
518
519  assume( L != NULL );
520  assume( R != NULL );
521  assume( R == currRing );
522
523  for( int k = IDELEMS(L) - 1; k >= 0; k-- )
524  {
525    const poly a = L->m[k]; assume( a != NULL );
526
527    // NOTE: label is k \in 0 ... |L|-1!!!
528    m_hash[p_GetComp(a, R)].push_back( new CLeadingTerm(k, a, R) );
529  }
530}
531
532
533CLCM::CLCM(const SchreyerSyzygyComputation& data): std::vector<bool>(), m_data(data), m_compute(false)
534{
535  const ideal& L = data.m_idLeads;
536  const ring&  R = data.m_rBaseRing;
537  const SchreyerSyzygyComputationFlags& attributes = data.m_atttributes;
538
539//  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
540//  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
541  const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
542  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
543
544
545  assume( L != NULL );
546  assume( R != NULL );
547  assume( R == currRing );
548
549  if( __TAILREDSYZ__ && !__HYBRIDNF__ )
550  {
551    const int l = IDELEMS(L);
552
553    resize(l, false);
554
555    const unsigned int N = rVar(R);
556
557    for( int k = l - 1; k >= 0; k-- )
558    {
559      const poly a = L->m[k]; assume( a != NULL );
560
561      for (unsigned int j = N; j > 0; j--)
562        if ( !(*this)[j] )
563          (*this)[j] = (p_GetExp(a, j, R) > 0);
564    }
565
566    m_compute = true;   
567  }
568}
569
570
571bool CLCM::Check(const poly m) const
572{
573  assume( m != NULL );
574  if( m_compute && (m != NULL))
575  { 
576    const ring& R = m_data.m_rBaseRing;
577    const SchreyerSyzygyComputationFlags& attributes = m_data.m_atttributes;
578
579  //  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
580  //  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
581    const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
582    const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
583
584    assume( R != NULL );
585    assume( R == currRing ); 
586
587    assume( __TAILREDSYZ__ && !__HYBRIDNF__ );
588   
589    const unsigned int N = rVar(R);
590
591    for (unsigned int j = N; j > 0; j--)
592      if ( (*this)[j] )
593        if(p_GetExp(m, j, R) > 0)
594          return true;
595
596    return false;
597   
598  } else return true;
599}
600
601void SchreyerSyzygyComputation::ComputeSyzygy()
602{
603//  FROM_NAMESPACE(INTERNAL, _ComputeSyzygy(m_idLeads, m_idTails, m_syzLeads, m_syzTails, m_rBaseRing, m_atttributes)); // TODO: just a wrapper for now :/
604
605  assume( m_idLeads != NULL );
606  assume( m_idTails != NULL );
607 
608  const ideal& L = m_idLeads;
609  const ideal& T = m_idTails;
610
611  ideal& TT = m_syzTails;
612  const ring& R = m_rBaseRing;
613  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
614
615//  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
616//  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
617  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
618  const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
619//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
620
621  assume( R == currRing ); // For attributes :-/
622
623  assume( IDELEMS(L) == IDELEMS(T) );
624
625  ComputeLeadingSyzygyTerms( __LEAD2SYZ__ ); // 2 terms OR 1 term!
626
627  ideal& LL = m_syzLeads;
628 
629  const int size = IDELEMS(LL);
630
631  TT = idInit(size, 0);
632
633  if( size == 1 && LL->m[0] == NULL )
634    return;
635
636
637  for( int k = size - 1; k >= 0; k-- )
638  {
639    const poly a = LL->m[k]; assume( a != NULL );
640
641    const int r = p_GetComp(a, R) - 1; 
642
643    assume( r >= 0 && r < IDELEMS(T) );
644    assume( r >= 0 && r < IDELEMS(L) );
645
646    poly aa = leadmonom(a, R); assume( aa != NULL); // :(   
647
648    poly a2 = pNext(a);   
649
650    // Splitting 2-terms Leading syzygy module
651    if( a2 != NULL )
652    {
653      TT->m[k] = a2; pNext(a) = NULL;
654    }
655
656    if( ! __HYBRIDNF__ )
657    {
658      poly t = TraverseTail(aa, T->m[r]);
659
660      if( a2 != NULL )
661      {
662        assume( __LEAD2SYZ__ );
663
664        const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
665
666        assume( r2 >= 0 && r2 < IDELEMS(T) );
667
668        TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2]), R), R);
669
670        p_Delete(&aa2, R);
671      } else
672        TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a), R);
673
674    } else
675    {
676      if( a2 == NULL )
677      {
678        aa = p_Mult_mm(aa, L->m[r], R); a2 = m_div.FindReducer(aa, a); 
679      }
680      assume( a2 != NULL );
681
682      TT->m[k] = SchreyerSyzygyNF(a, a2); // will copy a2 :(
683
684      p_Delete(&a2, R);
685    }
686    p_Delete(&aa, R);   
687  }
688
689  TT->rank = id_RankFreeModule(TT, R);
690}
691
692void SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms(bool bComputeSecondTerms)
693{
694  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
695
696  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
697  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
698
699  if( bComputeSecondTerms )
700  {
701    assume( __LEAD2SYZ__ );
702//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _Compute2LeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
703    m_syzLeads = Compute2LeadingSyzygyTerms();
704  }
705  else
706  {
707    assume( !__LEAD2SYZ__ );
708   
709    m_syzLeads = Compute1LeadingSyzygyTerms();
710  }
711//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _ComputeLeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
712 
713  // NOTE: set m_LS if tails are to be reduced!
714
715  if (__TAILREDSYZ__)
716    m_LS = m_syzLeads;
717
718  (void)( __LEAD2SYZ__ );
719}
720
721poly CReducerFinder::FindReducer(const poly product, const poly syzterm) const
722{
723//  return FROM_NAMESPACE(INTERNAL, _FindReducer(product, syzterm, m_idLeads, m_LS, m_rBaseRing, m_atttributes));
724//  poly _FindReducer(poly product, poly syzterm,
725//                     ideal L, ideal LS,
726//                     const ring r,
727//                     const SchreyerSyzygyComputationFlags attributes)
728
729  const ideal& L = m_data.m_idLeads;
730  const ideal& LS = m_data.m_LS;
731  const ring& r = m_data.m_rBaseRing;
732  const SchreyerSyzygyComputationFlags& attributes = m_data.m_atttributes;
733
734
735  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
736  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
737//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
738//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
739  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
740
741  assume( product != NULL );
742  assume( L != NULL );
743
744  long c = 0;
745
746  if (syzterm != NULL)
747    c = p_GetComp(syzterm, r) - 1;
748
749  assume( c >= 0 && c < IDELEMS(L) );
750 
751  if (__DEBUG__ && (syzterm != NULL))
752  {
753    const poly m = L->m[c];
754
755    assume( m != NULL ); assume( pNext(m) == NULL );
756
757    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
758    assume( p_EqualPolys(lm, product, r) );
759
760    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
761    //  def @@product = leadmonomial(syzterm) * L[@@r];
762
763    p_Delete(&lm, r);   
764  }
765
766  const long comp = p_GetComp(product, r);
767  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
768
769  assume( comp >= 0 );
770
771   // looking for an appropriate diviser p = L[k]...
772#if 1
773  CReducersHash::const_iterator it = m_hash.find(p_GetComp(product, r)); // same module component
774   
775  if( it == m_hash.end() )
776    return NULL;
777
778  const TReducers& reducers = it->second; 
779 
780  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
781  {
782     const poly p = (*vit)->m_lt;
783
784     assume( p_GetComp(p, r) == comp );
785     
786     const int k = (*vit)->m_label;
787
788     assume( L->m[k] == p );
789     
790     const unsigned long p_sev = (*vit)->m_sev;
791
792     assume( p_sev == p_GetShortExpVector(p, r) );     
793#else 
794   for( int k = IDELEMS(L)-1; k>= 0; k-- )
795   {
796     const poly p = L->m[k];
797     
798     if ( p_GetComp(p, r) != comp )
799       continue;
800       
801     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
802#endif
803     
804     if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
805       continue;     
806
807//     // ... which divides the product, looking for the _1st_ appropriate one!
808//     if( !p_LmDivisibleByNoComp(p, product, r) ) // included inside  p_LmShortDivisibleBy!
809//       continue;
810
811
812    const poly q = p_New(r);
813    pNext(q) = NULL;
814    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
815    p_SetComp(q, k + 1, r);
816    p_Setm(q, r);
817
818    // cannot allow something like: a*gen(i) - a*gen(i)
819    if (syzterm != NULL && (k == c))
820      if (p_ExpVectorEqual(syzterm, q, r))
821      {
822        if( __DEBUG__ )
823        {
824          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
825          dPrint(syzterm, r, r, 1);
826        }   
827
828        p_LmFree(q, r);
829        continue;
830      }
831
832    // while the complement (the fraction) is not reducible by leading syzygies
833    if( LS != NULL )
834    {
835      assume( __TAILREDSYZ__ );
836      BOOLEAN ok = TRUE;
837
838      // TODO: FindReducer in LS !!! there should be no divisors!
839      for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )
840      {
841        const poly pp = LS->m[kk];
842
843        if( p_LmDivisibleBy(pp, q, r) )
844        {
845
846          if( __DEBUG__ )
847          {
848            Print("_FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);
849            dPrint(pp, r, r, 1);
850          }   
851
852          ok = FALSE; // q in <LS> :((
853          break;                 
854        }
855      }
856
857      if(!ok)
858      {
859        p_LmFree(q, r);
860        continue;
861      }
862    }
863
864    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
865    return q;
866
867  }
868
869
870  return NULL;
871}
872
873
874poly SchreyerSyzygyComputation::SchreyerSyzygyNF(poly syz_lead, poly syz_2) const
875{
876// return FROM_NAMESPACE(INTERNAL, _SchreyerSyzygyNF(syz_lead, syz_2, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
877// poly _SchreyerSyzygyNF(poly syz_lead, poly syz_2,
878//                       ideal L, ideal T, ideal LS,
879//                       const ring r,
880//                       const SchreyerSyzygyComputationFlags attributes)
881// {
882
883  const ideal& L = m_idLeads;
884  const ideal& T = m_idTails;
885//  const ideal& LS = m_LS;
886  const ring& r = m_rBaseRing;
887
888//   const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
889//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
890//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
891//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
892//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
893//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
894
895  assume( syz_lead != NULL );
896  assume( syz_2 != NULL );
897
898  assume( L != NULL );
899  assume( T != NULL );
900
901  assume( IDELEMS(L) == IDELEMS(T) );
902
903  int  c = p_GetComp(syz_lead, r) - 1;
904
905  assume( c >= 0 && c < IDELEMS(T) );
906
907  poly p = leadmonom(syz_lead, r); // :( 
908  poly spoly = pp_Mult_qq(p, T->m[c], r);
909  p_Delete(&p, r);
910
911
912  c = p_GetComp(syz_2, r) - 1;
913  assume( c >= 0 && c < IDELEMS(T) );
914
915  p = leadmonom(syz_2, r); // :(
916  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
917  p_Delete(&p, r);
918
919  poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
920
921  while (spoly != NULL)
922  {
923    poly t = m_div.FindReducer(spoly, NULL);
924
925    p_LmDelete(&spoly, r);
926
927    if( t != NULL )
928    {
929      p = leadmonom(t, r); // :(
930      c = p_GetComp(t, r) - 1;
931
932      assume( c >= 0 && c < IDELEMS(T) );
933
934      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
935
936      p_Delete(&p, r);
937
938      tail = p_Add_q(tail, t, r);
939    }   
940  }
941
942  return tail;
943}
944
945
946poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
947{
948  const ideal& L = m_idLeads;
949  const ideal& T = m_idTails;
950//  const ideal& LS = m_LS;
951  const ring& r = m_rBaseRing;
952  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
953
954//  return FROM_NAMESPACE(INTERNAL, _TraverseTail(multiplier, tail, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
955// poly _TraverseTail(poly multiplier, poly tail,
956//                    ideal L, ideal T, ideal LS,
957//                    const ring r,
958//                    const SchreyerSyzygyComputationFlags attributes)
959// {
960
961
962//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
963//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
964//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
965//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
966   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
967
968  assume( multiplier != NULL );
969
970  assume( L != NULL );
971  assume( T != NULL );
972
973  poly s = NULL;
974
975  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
976    for(poly p = tail; p != NULL; p = pNext(p))   // iterate over the tail
977      s = p_Add_q(s, ReduceTerm(multiplier, p, NULL), r); 
978
979  return s;
980}
981
982
983
984
985poly SchreyerSyzygyComputation::ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck) const
986{
987  const ideal& L = m_idLeads;
988  const ideal& T = m_idTails;
989//  const ideal& LS = m_LS;
990  const ring& r = m_rBaseRing;
991  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
992
993//  return FROM_NAMESPACE(INTERNAL, _ReduceTerm(multiplier, term4reduction, syztermCheck, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
994// poly _ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
995//                 ideal L, ideal T, ideal LS,
996//                 const ring r,
997//                 const SchreyerSyzygyComputationFlags attributes)
998
999
1000
1001//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
1002//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
1003//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
1004//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
1005  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
1006
1007  assume( multiplier != NULL );
1008  assume( term4reduction != NULL );
1009
1010
1011  assume( L != NULL );
1012  assume( T != NULL );
1013
1014  // assume(r == currRing); // ?
1015
1016  // simple implementation with FindReducer:
1017  poly s = NULL;
1018
1019  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
1020  {
1021    // NOTE: only LT(term4reduction) should be used in the following:
1022    poly product = pp_Mult_mm(multiplier, term4reduction, r);
1023    s = m_div.FindReducer(product, syztermCheck);
1024    p_Delete(&product, r);
1025  }
1026
1027  if( s == NULL ) // No Reducer?
1028    return s;
1029
1030  poly b = leadmonom(s, r);
1031
1032  const int c = p_GetComp(s, r) - 1;
1033  assume( c >= 0 && c < IDELEMS(T) );
1034
1035  const poly tail = T->m[c];
1036
1037  if( tail != NULL )
1038    s = p_Add_q(s, TraverseTail(b, tail), r); 
1039
1040  return s;
1041}
1042
1043
1044
1045
1046
1047BEGIN_NAMESPACE_NONAME
1048   
1049static inline int atGetInt(idhdl rootRingHdl, const char* attribute, long def)
1050{
1051  return ((int)(long)(atGet(rootRingHdl, attribute, INT_CMD, (void*)def)));
1052}
1053
1054END_NAMESPACE   
1055
1056SchreyerSyzygyComputationFlags::SchreyerSyzygyComputationFlags(idhdl rootRingHdl):
1057#ifndef NDEBUG
1058    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", TRUE) ),
1059#else
1060    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", FALSE) ),
1061#endif
1062    __SYZCHECK__( (BOOLEAN)atGetInt(rootRingHdl, "SYZCHECK", __DEBUG__) ),
1063    __LEAD2SYZ__( (BOOLEAN)atGetInt(rootRingHdl, "LEAD2SYZ", 1) ), 
1064    __TAILREDSYZ__( (BOOLEAN)atGetInt(rootRingHdl, "TAILREDSYZ", 1) ), 
1065    __HYBRIDNF__( (BOOLEAN)atGetInt(rootRingHdl, "HYBRIDNF", 0) ) 
1066{   
1067  if( __DEBUG__ )
1068  {
1069    PrintS("SchreyerSyzygyComputationFlags: \n");
1070    Print("   DEBUG     : \t%d\n", __DEBUG__);
1071    Print("   SYZCHECK  : \t%d\n", __SYZCHECK__);
1072    Print("   LEAD2SYZ  : \t%d\n", __LEAD2SYZ__);
1073    Print("   TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
1074  }
1075
1076  // TODO: just current setting!
1077  assume( rootRingHdl == currRingHdl );
1078  assume( rootRingHdl->typ == RING_CMD );
1079  assume( rootRingHdl->data.uring == currRing );
1080  // move the global ring here inside???
1081}
1082
1083   
1084
1085
1086END_NAMESPACE               END_NAMESPACE_SINGULARXX
1087
1088
1089// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
Note: See TracBrowser for help on using the repository browser.