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

spielwiese
Last change on this file since 026171 was 026171, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
introduced a bitmask of `leading` variables NOTE: this is the 2nd Schreyer's optimization idea chg: further removal of unused variables
  • Property mode set to 100644
File size: 24.0 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
254  // TODO/NOTE: input is supposed to be (reverse-) sorted wrt "(c,ds)"!??
255
256  // components should come in groups: count elements in each group
257  // && estimate the real size!!!
258
259
260  // use just a vector instead???
261  const ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
262
263  int k = 0;
264
265  for (int j = 0; j < size; j++)
266  {
267    const poly p = id->m[j];
268    assume( p != NULL );
269    const int  c = p_GetComp(p, r);
270
271    for (int i = j - 1; i >= 0; i--)
272    {
273      const poly pp = id->m[i];
274      assume( pp != NULL );
275      const int  cc = p_GetComp(pp, r);
276
277      if( c != cc )
278        continue;
279
280      const poly m = p_Init(r); // p_New???
281
282      // m = LCM(p, pp) / p! // TODO: optimize: knowing the ring structure: (C/lp)!
283      for (int v = rVar(r); v > 0; v--)
284      {
285        assume( v > 0 );
286        assume( v <= rVar(r) );
287
288        const short e1 = p_GetExp(p , v, r);
289        const short e2 = p_GetExp(pp, v, r);
290
291        if( e1 >= e2 )
292          p_SetExp(m, v, 0, r);
293        else
294          p_SetExp(m, v, e2 - e1, r);
295
296      }
297
298      assume( (j > i) && (i >= 0) );
299
300      p_SetComp(m, j + 1, r);
301      pNext(m) = NULL;
302      p_SetCoeff0(m, n_Init(1, r->cf), r); // for later...
303
304      p_Setm(m, r); // should not do anything!!!
305
306      newid->m[k++] = m;
307    }
308  }
309
310//   if( __DEBUG__ && FALSE )
311//   {
312//     PrintS("ComputeLeadingSyzygyTerms::Temp0: \n");
313//     dPrint(newid, r, r, 1);
314//   }
315
316  // the rest of newid is assumed to be zeroes...
317
318  // simplify(newid, 2 + 32)??
319  // sort(newid, "C,ds")[1]???
320  id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
321
322//   if( __DEBUG__ && FALSE )
323//   {
324//     PrintS("ComputeLeadingSyzygyTerms::Temp1: \n");
325//     dPrint(newid, r, r, 1);
326//   }
327
328  idSkipZeroes(newid); // #define SIMPL_NULL 2
329
330//   if( __DEBUG__ )
331//   {
332//     PrintS("ComputeLeadingSyzygyTerms::Output: \n");
333//     dPrint(newid, r, r, 1);
334//   }
335
336  Sort_c_ds(newid, r);
337
338  return newid;
339}
340
341ideal SchreyerSyzygyComputation::Compute2LeadingSyzygyTerms()
342{
343  const ideal& id = m_idLeads;
344  const ring& r = m_rBaseRing;
345  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
346 
347//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
348//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
349//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
350//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
351  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
352 
353
354  // 1. set of components S?
355  // 2. for each component c from S: set of indices of leading terms
356  // with this component?
357  // 3. short exp. vectors for each leading term?
358
359  const int size = IDELEMS(id);
360
361  if( size < 2 )
362  {
363    const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...   
364    return newid;
365  }
366
367
368  // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
369 
370  // components should come in groups: count elements in each group
371  // && estimate the real size!!!
372
373
374  // use just a vector instead???
375  ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
376
377  int k = 0;
378
379  for (int j = 0; j < size; j++)
380  {
381    const poly p = id->m[j];
382    assume( p != NULL );
383    const int  c = p_GetComp(p, r);
384
385    for (int i = j - 1; i >= 0; i--)
386    {
387      const poly pp = id->m[i];
388      assume( pp != NULL );
389      const int  cc = p_GetComp(pp, r);
390
391      if( c != cc )
392        continue;
393
394        // allocate memory & zero it out!
395      const poly m = p_Init(r); const poly mm = p_Init(r);
396
397
398        // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
399        // TODO: optimize: knowing the ring structure: (C/lp)!
400
401      for (int v = rVar(r); v > 0; v--)
402      {
403        assume( v > 0 );
404        assume( v <= rVar(r) );
405
406        const short e1 = p_GetExp(p , v, r);
407        const short e2 = p_GetExp(pp, v, r);
408
409        if( e1 >= e2 )
410          p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
411        else
412          p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
413
414      }
415
416      assume( (j > i) && (i >= 0) );
417
418      p_SetComp(m, j + 1, r);
419      p_SetComp(mm, i + 1, r);
420
421      const number& lc1 = p_GetCoeff(p , r);
422      const number& lc2 = p_GetCoeff(pp, r);
423
424      number g = n_Lcm( lc1, lc2, r );
425
426      p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
427      p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);
428
429      n_Delete(&g, r);
430
431      p_Setm(m, r); // should not do anything!!!
432      p_Setm(mm, r); // should not do anything!!!
433
434      pNext(m) = mm; //        pNext(mm) = NULL;
435
436      newid->m[k++] = m;
437    }
438  }
439
440//   if( __DEBUG__ && FALSE )
441//   {
442//     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
443//     dPrint(newid, r, r, 1);
444//   }
445
446  if( !__TAILREDSYZ__ )
447  {
448      // simplify(newid, 2 + 32)??
449      // sort(newid, "C,ds")[1]???
450    id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
451
452//     if( __DEBUG__ && FALSE )
453//     {
454//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
455//       dPrint(newid, r, r, 1);
456//     }
457  }
458  else
459  {
460      //      option(redSB); option(redTail);
461      //      TEST_OPT_REDSB
462      //      TEST_OPT_REDTAIL
463    assume( r == currRing );
464    BITSET _save_test = test;
465    test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB));
466
467    intvec* w=new intvec(IDELEMS(newid));
468    ideal tmp = kStd(newid, currQuotient, isHomog, &w);
469    delete w;
470
471    test = _save_test;
472
473    id_Delete(&newid, r);
474    newid = tmp;
475
476//     if( __DEBUG__ && FALSE )
477//     {
478//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
479//       dPrint(newid, r, r, 1);
480//     }
481
482  }
483
484  idSkipZeroes(newid);
485
486  Sort_c_ds(newid, r);
487 
488  return newid;
489}
490
491CLCM::CLCM(const SchreyerSyzygyComputation& data): std::vector<bool>(), m_data(data), m_compute(false)
492{
493  const ideal& L = data.m_idLeads;
494  const ring&  R = data.m_rBaseRing;
495  const SchreyerSyzygyComputationFlags& attributes = data.m_atttributes;
496
497//  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
498//  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
499  const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
500  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
501
502
503  assume( L != NULL );
504  assume( R != NULL );
505  assume( R == currRing );
506
507  if( __TAILREDSYZ__ && !__HYBRIDNF__ )
508  {
509    const int l = IDELEMS(L);
510
511    resize(l, false);
512
513    const unsigned int N = rVar(R);
514
515    for( int k = l - 1; k >= 0; k-- )
516    {
517      const poly a = L->m[k]; assume( a != NULL );
518
519      for (unsigned int j = N; j > 0; j--)
520        if ( !(*this)[j] )
521          (*this)[j] = (p_GetExp(a, j, R) > 0);
522    }
523
524    m_compute = true;   
525  }
526}
527
528
529bool CLCM::Check(const poly m) const
530{
531  assume( m != NULL );
532  if( m_compute && (m != NULL))
533  { 
534    const ring& R = m_data.m_rBaseRing;
535    const SchreyerSyzygyComputationFlags& attributes = m_data.m_atttributes;
536
537  //  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
538  //  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
539    const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
540    const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
541
542    assume( R != NULL );
543    assume( R == currRing ); 
544
545    assume( __TAILREDSYZ__ && !__HYBRIDNF__ );
546   
547    const unsigned int N = rVar(R);
548
549    for (unsigned int j = N; j > 0; j--)
550      if ( (*this)[j] )
551        if(p_GetExp(m, j, R) > 0)
552          return true;
553
554    return false;
555   
556  } else return true;
557}
558
559void SchreyerSyzygyComputation::ComputeSyzygy()
560{
561//  FROM_NAMESPACE(INTERNAL, _ComputeSyzygy(m_idLeads, m_idTails, m_syzLeads, m_syzTails, m_rBaseRing, m_atttributes)); // TODO: just a wrapper for now :/
562
563  assume( m_idLeads != NULL );
564  assume( m_idTails != NULL );
565 
566  const ideal& L = m_idLeads;
567  const ideal& T = m_idTails;
568
569  ideal& TT = m_syzTails;
570  const ring& R = m_rBaseRing;
571  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
572
573//  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
574//  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
575  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
576  const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
577//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
578
579  assume( R == currRing ); // For attributes :-/
580
581  assume( IDELEMS(L) == IDELEMS(T) );
582
583  ComputeLeadingSyzygyTerms( __LEAD2SYZ__ ); // 2 terms OR 1 term!
584
585  ideal& LL = m_syzLeads;
586 
587  const int size = IDELEMS(LL);
588
589  TT = idInit(size, 0);
590
591  if( size == 1 && LL->m[0] == NULL )
592    return;
593
594
595  for( int k = size - 1; k >= 0; k-- )
596  {
597    const poly a = LL->m[k]; assume( a != NULL );
598
599    const int r = p_GetComp(a, R) - 1; 
600
601    assume( r >= 0 && r < IDELEMS(T) );
602    assume( r >= 0 && r < IDELEMS(L) );
603
604    poly aa = leadmonom(a, R); assume( aa != NULL); // :(   
605
606    poly a2 = pNext(a);   
607
608    if( a2 != NULL )
609    {
610      TT->m[k] = a2; pNext(a) = NULL;
611    }
612
613    if( ! __HYBRIDNF__ )
614    {
615      poly t = TraverseTail(aa, T->m[r]);
616
617      if( a2 != NULL )
618      {
619        assume( __LEAD2SYZ__ );
620
621        const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
622
623        assume( r2 >= 0 && r2 < IDELEMS(T) );
624
625        TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2]), R), R);
626
627        p_Delete(&aa2, R);
628      } else
629        TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a), R);
630
631    } else
632    {
633      if( a2 == NULL )
634      {
635        aa = p_Mult_mm(aa, L->m[r], R);
636        a2 = FindReducer(aa, a); 
637      }
638      assume( a2 != NULL );
639
640      TT->m[k] = SchreyerSyzygyNF(a, a2); // will copy a2 :(
641
642      p_Delete(&a2, R);
643    }
644    p_Delete(&aa, R);   
645  }
646
647  TT->rank = id_RankFreeModule(TT, R);
648}
649
650void SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms(bool bComputeSecondTerms)
651{
652  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
653
654  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
655  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
656
657  if( bComputeSecondTerms )
658  {
659    assume( __LEAD2SYZ__ );
660//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _Compute2LeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
661    m_syzLeads = Compute2LeadingSyzygyTerms();
662  }
663  else
664  {
665    assume( !__LEAD2SYZ__ );
666   
667    m_syzLeads = Compute1LeadingSyzygyTerms();
668  }
669//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _ComputeLeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
670 
671  // NOTE: set m_LS if tails are to be reduced!
672
673  if (__TAILREDSYZ__)
674    m_LS = m_syzLeads;
675
676  (void)( __LEAD2SYZ__ );
677}
678
679poly SchreyerSyzygyComputation::FindReducer(poly product, poly syzterm) const
680{
681//  return FROM_NAMESPACE(INTERNAL, _FindReducer(product, syzterm, m_idLeads, m_LS, m_rBaseRing, m_atttributes));
682//  poly _FindReducer(poly product, poly syzterm,
683//                     ideal L, ideal LS,
684//                     const ring r,
685//                     const SchreyerSyzygyComputationFlags attributes)
686
687  const ideal& L = m_idLeads;
688  const ideal& LS = m_LS;
689  const ring& r = m_rBaseRing;
690  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
691
692
693  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
694  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
695//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
696//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
697  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
698
699  assume( product != NULL );
700  assume( L != NULL );
701
702  int c = 0;
703
704  if (syzterm != NULL)
705    c = p_GetComp(syzterm, r) - 1;
706
707  assume( c >= 0 && c < IDELEMS(L) );
708
709  if (__SYZCHECK__ && syzterm != NULL)
710  {
711    const poly m = L->m[c];
712
713    assume( m != NULL ); assume( pNext(m) == NULL );
714
715    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
716    assume( p_EqualPolys(lm, product, r) );
717
718    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
719    //  def @@product = leadmonomial(syzterm) * L[@@r];
720
721    p_Delete(&lm, r);   
722  }
723
724  // looking for an appropriate diviser q = L[k]...
725  for( int k = IDELEMS(L)-1; k>= 0; k-- )
726  {
727    const poly p = L->m[k];   
728
729    // ... which divides the product, looking for the _1st_ appropriate one!
730    if( !p_LmDivisibleBy(p, product, r) )
731      continue;
732
733
734    const poly q = p_New(r);
735    pNext(q) = NULL;
736    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
737    p_SetComp(q, k + 1, r);
738    p_Setm(q, r);
739
740    // cannot allow something like: a*gen(i) - a*gen(i)
741    if (syzterm != NULL && (k == c))
742      if (p_ExpVectorEqual(syzterm, q, r))
743      {
744        if( __DEBUG__ )
745        {
746          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
747          dPrint(syzterm, r, r, 1);
748        }   
749
750        p_LmFree(q, r);
751        continue;
752      }
753
754    // while the complement (the fraction) is not reducible by leading syzygies
755    if( LS != NULL )
756    {
757      assume( __TAILREDSYZ__ );
758      BOOLEAN ok = TRUE;
759
760      for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )
761      {
762        const poly pp = LS->m[kk];
763
764        if( p_LmDivisibleBy(pp, q, r) )
765        {
766
767          if( __DEBUG__ )
768          {
769            Print("_FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);
770            dPrint(pp, r, r, 1);
771          }   
772
773          ok = FALSE; // q in <LS> :((
774          break;                 
775        }
776      }
777
778      if(!ok)
779      {
780        p_LmFree(q, r);
781        continue;
782      }
783    }
784
785    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
786    return q;
787
788  }
789
790
791  return NULL;
792}
793
794
795poly SchreyerSyzygyComputation::SchreyerSyzygyNF(poly syz_lead, poly syz_2) const
796{
797// return FROM_NAMESPACE(INTERNAL, _SchreyerSyzygyNF(syz_lead, syz_2, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
798// poly _SchreyerSyzygyNF(poly syz_lead, poly syz_2,
799//                       ideal L, ideal T, ideal LS,
800//                       const ring r,
801//                       const SchreyerSyzygyComputationFlags attributes)
802// {
803
804  const ideal& L = m_idLeads;
805  const ideal& T = m_idTails;
806//  const ideal& LS = m_LS;
807  const ring& r = m_rBaseRing;
808
809//   const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
810//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
811//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
812//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
813//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
814//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
815
816  assume( syz_lead != NULL );
817  assume( syz_2 != NULL );
818
819  assume( L != NULL );
820  assume( T != NULL );
821
822  assume( IDELEMS(L) == IDELEMS(T) );
823
824  int  c = p_GetComp(syz_lead, r) - 1;
825
826  assume( c >= 0 && c < IDELEMS(T) );
827
828  poly p = leadmonom(syz_lead, r); // :( 
829  poly spoly = pp_Mult_qq(p, T->m[c], r);
830  p_Delete(&p, r);
831
832
833  c = p_GetComp(syz_2, r) - 1;
834  assume( c >= 0 && c < IDELEMS(T) );
835
836  p = leadmonom(syz_2, r); // :(
837  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
838  p_Delete(&p, r);
839
840  poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
841
842  while (spoly != NULL)
843  {
844    poly t = FindReducer(spoly, NULL);
845
846    p_LmDelete(&spoly, r);
847
848    if( t != NULL )
849    {
850      p = leadmonom(t, r); // :(
851      c = p_GetComp(t, r) - 1;
852
853      assume( c >= 0 && c < IDELEMS(T) );
854
855      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
856
857      p_Delete(&p, r);
858
859      tail = p_Add_q(tail, t, r);
860    }   
861  }
862
863  return tail;
864}
865
866
867poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
868{
869  const ideal& L = m_idLeads;
870  const ideal& T = m_idTails;
871//  const ideal& LS = m_LS;
872  const ring& r = m_rBaseRing;
873  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
874
875//  return FROM_NAMESPACE(INTERNAL, _TraverseTail(multiplier, tail, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
876// poly _TraverseTail(poly multiplier, poly tail,
877//                    ideal L, ideal T, ideal LS,
878//                    const ring r,
879//                    const SchreyerSyzygyComputationFlags attributes)
880// {
881
882
883//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
884//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
885//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
886//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
887   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
888
889  assume( multiplier != NULL );
890
891  assume( L != NULL );
892  assume( T != NULL );
893
894  poly s = NULL;
895
896  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
897    for(poly p = tail; p != NULL; p = pNext(p))   // iterate over the tail
898      s = p_Add_q(s, ReduceTerm(multiplier, p, NULL), r); 
899
900  return s;
901}
902
903
904
905
906poly SchreyerSyzygyComputation::ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck) const
907{
908  const ideal& L = m_idLeads;
909  const ideal& T = m_idTails;
910//  const ideal& LS = m_LS;
911  const ring& r = m_rBaseRing;
912  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
913
914//  return FROM_NAMESPACE(INTERNAL, _ReduceTerm(multiplier, term4reduction, syztermCheck, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
915// poly _ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
916//                 ideal L, ideal T, ideal LS,
917//                 const ring r,
918//                 const SchreyerSyzygyComputationFlags attributes)
919
920
921
922//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
923//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
924//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
925//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
926  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
927
928  assume( multiplier != NULL );
929  assume( term4reduction != NULL );
930
931
932  assume( L != NULL );
933  assume( T != NULL );
934
935  // assume(r == currRing); // ?
936
937  // simple implementation with FindReducer:
938  poly s = NULL;
939
940  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
941  {
942    // NOTE: only LT(term4reduction) should be used in the following:
943    poly product = pp_Mult_mm(multiplier, term4reduction, r);
944    s = FindReducer(product, syztermCheck);
945    p_Delete(&product, r);
946  }
947
948  if( s == NULL ) // No Reducer?
949    return s;
950
951  poly b = leadmonom(s, r);
952
953  const int c = p_GetComp(s, r) - 1;
954  assume( c >= 0 && c < IDELEMS(T) );
955
956  const poly tail = T->m[c];
957
958  if( tail != NULL )
959    s = p_Add_q(s, TraverseTail(b, tail), r); 
960
961  return s;
962}
963
964
965
966
967
968BEGIN_NAMESPACE_NONAME
969   
970static inline int atGetInt(idhdl rootRingHdl, const char* attribute, long def)
971{
972  return ((int)(long)(atGet(rootRingHdl, attribute, INT_CMD, (void*)def)));
973}
974
975END_NAMESPACE   
976
977SchreyerSyzygyComputationFlags::SchreyerSyzygyComputationFlags(idhdl rootRingHdl):
978#ifndef NDEBUG
979    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", TRUE) ),
980#else
981    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", FALSE) ),
982#endif
983    __SYZCHECK__( (BOOLEAN)atGetInt(rootRingHdl, "SYZCHECK", __DEBUG__) ),
984    __LEAD2SYZ__( (BOOLEAN)atGetInt(rootRingHdl, "LEAD2SYZ", 1) ), 
985    __TAILREDSYZ__( (BOOLEAN)atGetInt(rootRingHdl, "TAILREDSYZ", 1) ), 
986    __HYBRIDNF__( (BOOLEAN)atGetInt(rootRingHdl, "HYBRIDNF", 0) ) 
987{   
988  if( __DEBUG__ )
989  {
990    PrintS("SchreyerSyzygyComputationFlags: \n");
991    Print("   DEBUG     : \t%d\n", __DEBUG__);
992    Print("   SYZCHECK  : \t%d\n", __SYZCHECK__);
993    Print("   LEAD2SYZ  : \t%d\n", __LEAD2SYZ__);
994    Print("   TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
995  }
996
997  // TODO: just current setting!
998  assume( rootRingHdl == currRingHdl );
999  assume( rootRingHdl->typ == RING_CMD );
1000  assume( rootRingHdl->data.uring == currRing );
1001  // move the global ring here inside???
1002}
1003
1004   
1005
1006
1007END_NAMESPACE               END_NAMESPACE_SINGULARXX
1008
1009
1010// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
Note: See TracBrowser for help on using the repository browser.