source: git/dyn_modules/syzextra/syzextra.cc @ 4eba3ad

jengelh-datetimespielwiese
Last change on this file since 4eba3ad was 4eba3ad, checked in by Oleksandr Motsak <motsak@…>, 9 years ago
introduced SchreyerSyzygyComputationFlags for storing and passing ring attributes within SchreyerSyzygyComputatio chg: renamed internal functions to avoid name conflicts chg: better namespace handling
  • Property mode set to 100644
File size: 21.6 KB
Line 
1// -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2/*****************************************************************************\
3 * Computer Algebra System SINGULAR   
4\*****************************************************************************/
5/** @file syzextra.cc
6 *
7 * Here we implement the Computation of Syzygies
8 *
9 * ABSTRACT: Computation of Syzygies due to Schreyer
10 *
11 * @author Oleksandr Motsak
12 *
13 **/
14/*****************************************************************************/
15
16// include header file
17#include <kernel/mod2.h>
18
19#include "syzextra.h"
20
21#include "DebugPrint.h"
22
23#include <omalloc/omalloc.h>
24
25#include <misc/intvec.h>
26#include <misc/options.h>
27
28#include <coeffs/coeffs.h>
29
30#include <polys/monomials/p_polys.h>
31#include <polys/monomials/ring.h>
32
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
225BEGIN_NAMESPACE(INTERNAL)
226
227ideal _ComputeLeadingSyzygyTerms(const ideal& id,
228                                const ring r,
229                                const SchreyerSyzygyComputationFlags attributes)
230{
231//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
232//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
233   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
234//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
235//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
236
237   assume(!__LEAD2SYZ__);
238
239  // 1. set of components S?
240  // 2. for each component c from S: set of indices of leading terms
241  // with this component?
242  // 3. short exp. vectors for each leading term?
243
244  const int size = IDELEMS(id);
245
246  if( size < 2 )
247  {
248    const ideal newid = idInit(1, 0); newid->m[0] = NULL; // zero ideal...       
249    return newid;
250  }
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 _Compute2LeadingSyzygyTerms(const ideal& id,
341                                 const ring r,
342                                 const SchreyerSyzygyComputationFlags attributes)
343{
344//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
345//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
346//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
347//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
348  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
349 
350
351  // 1. set of components S?
352  // 2. for each component c from S: set of indices of leading terms
353  // with this component?
354  // 3. short exp. vectors for each leading term?
355
356  const int size = IDELEMS(id);
357
358  if( size < 2 )
359  {
360    const ideal newid = idInit(1, 1); newid->m[0] = NULL; // zero module...   
361    return newid;
362  }
363
364
365    // TODO/NOTE: input is supposed to be sorted wrt "C,ds"!??
366
367    // components should come in groups: count elements in each group
368    // && estimate the real size!!!
369
370
371    // use just a vector instead???
372  ideal newid = idInit( (size * (size-1))/2, size); // maximal size: ideal case!
373
374  int k = 0;
375
376  for (int j = 0; j < size; j++)
377  {
378    const poly p = id->m[j];
379    assume( p != NULL );
380    const int  c = p_GetComp(p, r);
381
382    for (int i = j - 1; i >= 0; i--)
383    {
384      const poly pp = id->m[i];
385      assume( pp != NULL );
386      const int  cc = p_GetComp(pp, r);
387
388      if( c != cc )
389        continue;
390
391        // allocate memory & zero it out!
392      const poly m = p_Init(r); const poly mm = p_Init(r);
393
394
395        // m = LCM(p, pp) / p! mm = LCM(p, pp) / pp!
396        // TODO: optimize: knowing the ring structure: (C/lp)!
397
398      for (int v = rVar(r); v > 0; v--)
399      {
400        assume( v > 0 );
401        assume( v <= rVar(r) );
402
403        const short e1 = p_GetExp(p , v, r);
404        const short e2 = p_GetExp(pp, v, r);
405
406        if( e1 >= e2 )
407          p_SetExp(mm, v, e1 - e2, r); //            p_SetExp(m, v, 0, r);
408        else
409          p_SetExp(m, v, e2 - e1, r); //            p_SetExp(mm, v, 0, r);
410
411      }
412
413      assume( (j > i) && (i >= 0) );
414
415      p_SetComp(m, j + 1, r);
416      p_SetComp(mm, i + 1, r);
417
418      const number& lc1 = p_GetCoeff(p , r);
419      const number& lc2 = p_GetCoeff(pp, r);
420
421      number g = n_Lcm( lc1, lc2, r );
422
423      p_SetCoeff0(m ,       n_Div(g, lc1, r), r);
424      p_SetCoeff0(mm, n_Neg(n_Div(g, lc2, r), r), r);
425
426      n_Delete(&g, r);
427
428      p_Setm(m, r); // should not do anything!!!
429      p_Setm(mm, r); // should not do anything!!!
430
431      pNext(m) = mm; //        pNext(mm) = NULL;
432
433      newid->m[k++] = m;
434    }
435  }
436
437//   if( __DEBUG__ && FALSE )
438//   {
439//     PrintS("Compute2LeadingSyzygyTerms::Temp0: \n");
440//     dPrint(newid, r, r, 1);
441//   }
442
443  if( !__TAILREDSYZ__ )
444  {
445      // simplify(newid, 2 + 32)??
446      // sort(newid, "C,ds")[1]???
447    id_DelDiv(newid, r); // #define SIMPL_LMDIV 32
448
449//     if( __DEBUG__ && FALSE )
450//     {
451//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (deldiv): \n");
452//       dPrint(newid, r, r, 1);
453//     }
454  }
455  else
456  {
457      //      option(redSB); option(redTail);
458      //      TEST_OPT_REDSB
459      //      TEST_OPT_REDTAIL
460    assume( r == currRing );
461    BITSET _save_test = test;
462    test |= (Sy_bit(OPT_REDTAIL) | Sy_bit(OPT_REDSB));
463
464    intvec* w=new intvec(IDELEMS(newid));
465    ideal tmp = kStd(newid, currQuotient, isHomog, &w);
466    delete w;
467
468    test = _save_test;
469
470    id_Delete(&newid, r);
471    newid = tmp;
472
473//     if( __DEBUG__ && FALSE )
474//     {
475//       PrintS("Compute2LeadingSyzygyTerms::Temp1 (std): \n");
476//       dPrint(newid, r, r, 1);
477//     }
478
479  }
480
481  idSkipZeroes(newid);
482
483  Sort_c_ds(newid, r);
484 
485  return newid;
486}
487
488poly _FindReducer(poly product, poly syzterm,
489                 ideal L, ideal LS,
490                 const ring r,
491                 const SchreyerSyzygyComputationFlags attributes)
492{
493  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
494  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
495//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
496//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
497//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
498
499  assume( product != NULL );
500  assume( L != NULL );
501
502  int c = 0;
503
504  if (syzterm != NULL)
505    c = p_GetComp(syzterm, r) - 1;
506
507  assume( c >= 0 && c < IDELEMS(L) );
508
509  if (__SYZCHECK__ && syzterm != NULL)
510  {
511    const poly m = L->m[c];
512
513    assume( m != NULL ); assume( pNext(m) == NULL );
514
515    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
516    assume( p_EqualPolys(lm, product, r) );
517
518    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
519    //  def @@product = leadmonomial(syzterm) * L[@@r];
520
521    p_Delete(&lm, r);   
522  }
523
524  // looking for an appropriate diviser q = L[k]...
525  for( int k = IDELEMS(L)-1; k>= 0; k-- )
526  {
527    const poly p = L->m[k];   
528
529    // ... which divides the product, looking for the _1st_ appropriate one!
530    if( !p_LmDivisibleBy(p, product, r) )
531      continue;
532
533
534    const poly q = p_New(r);
535    pNext(q) = NULL;
536    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
537    p_SetComp(q, k + 1, r);
538    p_Setm(q, r);
539
540    // cannot allow something like: a*gen(i) - a*gen(i)
541    if (syzterm != NULL && (k == c))
542      if (p_ExpVectorEqual(syzterm, q, r))
543      {
544        if( __DEBUG__ )
545        {
546          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
547          dPrint(syzterm, r, r, 1);
548        }   
549
550        p_LmFree(q, r);
551        continue;
552      }
553
554    // while the complement (the fraction) is not reducible by leading syzygies
555    if( LS != NULL )
556    {
557      BOOLEAN ok = TRUE;
558
559      for(int kk = IDELEMS(LS)-1; kk>= 0; kk-- )
560      {
561        const poly pp = LS->m[kk];
562
563        if( p_LmDivisibleBy(pp, q, r) )
564        {
565
566          if( __DEBUG__ )
567          {
568            Print("_FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", kk+1);
569            dPrint(pp, r, r, 1);
570          }   
571
572          ok = FALSE; // q in <LS> :((
573          break;                 
574        }
575      }
576
577      if(!ok)
578      {
579        p_LmFree(q, r);
580        continue;
581      }
582    }
583
584    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
585    return q;
586
587  }
588
589
590  return NULL;
591}
592
593poly _ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
594                ideal L, ideal T, ideal LS,
595                const ring r, const SchreyerSyzygyComputationFlags attributes);
596
597poly _TraverseTail(poly multiplier, poly tail, 
598                  ideal L, ideal T, ideal LS,
599                  const ring r,
600                  const SchreyerSyzygyComputationFlags attributes)
601{
602//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
603//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
604//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
605//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
606//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
607
608  assume( multiplier != NULL );
609
610  assume( L != NULL );
611  assume( T != NULL );
612
613  poly s = NULL;
614
615  // iterate over the tail
616  for(poly p = tail; p != NULL; p = pNext(p))
617    s = p_Add_q(s, _ReduceTerm(multiplier, p, NULL, L, T, LS, r, attributes), r); 
618
619  return s;
620}
621
622
623poly _ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
624                ideal L, ideal T, ideal LS,
625                const ring r,
626                const SchreyerSyzygyComputationFlags attributes)
627{
628//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
629//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
630//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
631//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
632//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
633
634  assume( multiplier != NULL );
635  assume( term4reduction != NULL );
636
637
638  assume( L != NULL );
639  assume( T != NULL );
640
641  // assume(r == currRing); // ?
642
643  // simple implementation with FindReducer:
644  poly s = NULL;
645
646  if(1)
647  {
648    // NOTE: only LT(term4reduction) should be used in the following:
649    poly product = pp_Mult_mm(multiplier, term4reduction, r);
650    s = _FindReducer(product, syztermCheck, L, LS, r, attributes);
651    p_Delete(&product, r);
652  }
653
654  if( s == NULL ) // No Reducer?
655    return s;
656
657  poly b = leadmonom(s, r);
658
659  const int c = p_GetComp(s, r) - 1;
660  assume( c >= 0 && c < IDELEMS(T) );
661
662  const poly tail = T->m[c];
663
664  if( tail != NULL )
665    s = p_Add_q(s, _TraverseTail(b, tail, L, T, LS, r, attributes), r); 
666
667  return s;
668}
669
670
671poly _SchreyerSyzygyNF(poly syz_lead, poly syz_2,
672                      ideal L, ideal T, ideal LS,
673                      const ring r,
674                      const SchreyerSyzygyComputationFlags attributes)
675{
676//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
677//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
678//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
679//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
680//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
681
682  assume( syz_lead != NULL );
683  assume( syz_2 != NULL );
684
685  assume( L != NULL );
686  assume( T != NULL );
687
688  assume( IDELEMS(L) == IDELEMS(T) );
689
690  int  c = p_GetComp(syz_lead, r) - 1;
691
692  assume( c >= 0 && c < IDELEMS(T) );
693
694  poly p = leadmonom(syz_lead, r); // :( 
695  poly spoly = pp_Mult_qq(p, T->m[c], r);
696  p_Delete(&p, r);
697
698
699  c = p_GetComp(syz_2, r) - 1;
700  assume( c >= 0 && c < IDELEMS(T) );
701
702  p = leadmonom(syz_2, r); // :(
703  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
704  p_Delete(&p, r);
705
706  poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
707
708  while (spoly != NULL)
709  {
710    poly t = _FindReducer(spoly, NULL, L, LS, r, attributes);
711
712    p_LmDelete(&spoly, r);
713
714    if( t != NULL )
715    {
716      p = leadmonom(t, r); // :(
717      c = p_GetComp(t, r) - 1;
718
719      assume( c >= 0 && c < IDELEMS(T) );
720
721      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
722
723      p_Delete(&p, r);
724
725      tail = p_Add_q(tail, t, r);
726    }   
727  }
728
729  return tail;
730}
731
732
733void _ComputeSyzygy(const ideal L, const ideal T,
734                   ideal& LL, ideal& TT,
735                   const ring R,
736                   const SchreyerSyzygyComputationFlags attributes)
737{
738//  const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
739//  const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
740  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
741  const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
742  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
743
744  assume( R == currRing ); // For attributes :-/
745
746  assume( IDELEMS(L) == IDELEMS(T) );
747
748  if( __LEAD2SYZ__ )
749    LL = _Compute2LeadingSyzygyTerms(L, R, attributes); // 2 terms!
750  else
751    LL = _ComputeLeadingSyzygyTerms(L, R, attributes); // 1 term!
752
753  const int size = IDELEMS(LL);
754
755  TT = idInit(size, 0);
756
757  if( size == 1 && LL->m[0] == NULL )
758    return;
759
760
761  ideal LS = NULL;
762
763  if( __TAILREDSYZ__ )
764    LS = LL;
765
766  for( int k = size - 1; k >= 0; k-- )
767  {
768    const poly a = LL->m[k]; assume( a != NULL );
769
770    const int r = p_GetComp(a, R) - 1; 
771
772    assume( r >= 0 && r < IDELEMS(T) );
773    assume( r >= 0 && r < IDELEMS(L) );
774
775    poly aa = leadmonom(a, R); assume( aa != NULL); // :(   
776
777    poly a2 = pNext(a);   
778
779    if( a2 != NULL )
780    {
781      TT->m[k] = a2; pNext(a) = NULL;
782    }
783
784    if( ! __HYBRIDNF__ )
785    {
786      poly t = _TraverseTail(aa, T->m[r], L, T, LS, R, attributes);
787
788      if( a2 != NULL )
789      {
790        assume( __LEAD2SYZ__ );
791
792        const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
793
794        assume( r2 >= 0 && r2 < IDELEMS(T) );
795
796        TT->m[k] = p_Add_q(a2, p_Add_q(t, _TraverseTail(aa2, T->m[r2], L, T, LS, R, attributes), R), R);
797
798        p_Delete(&aa2, R);
799      } else
800        TT->m[k] = p_Add_q(t, _ReduceTerm(aa, L->m[r], a, L, T, LS, R, attributes), R);
801
802    } else
803    {
804      if( a2 == NULL )
805      {
806        aa = p_Mult_mm(aa, L->m[r], R);
807        a2 = _FindReducer(aa, a, L, LS, R, attributes); 
808      }
809      assume( a2 != NULL );
810
811      TT->m[k] = _SchreyerSyzygyNF(a, a2, L, T, LS, R, attributes); // will copy a2 :(
812
813      p_Delete(&a2, R);
814    }
815    p_Delete(&aa, R);   
816  }
817
818  TT->rank = id_RankFreeModule(TT, R);
819}
820
821END_NAMESPACE
822
823
824// TODO: wrapper layer, just for now... dissolve!
825
826void SchreyerSyzygyComputation::ComputeSyzygy()
827{
828  /// assumes m_syzLeads == m_syzTails == NULL!
829  FROM_NAMESPACE(INTERNAL, _ComputeSyzygy(m_idLeads, m_idTails, m_syzLeads, m_syzTails, m_rBaseRing, m_atttributes)); // TODO: just a wrapper for now :/
830}
831
832void SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms(bool bComputeSecondTerms)
833{
834  if( bComputeSecondTerms )
835    m_syzLeads = FROM_NAMESPACE(INTERNAL, _Compute2LeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
836  else
837    m_syzLeads = FROM_NAMESPACE(INTERNAL, _ComputeLeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
838 
839  // NOTE: set m_LS if tails are to be reduced!
840}
841
842poly SchreyerSyzygyComputation::FindReducer(poly product, poly syzterm) const
843{
844  return FROM_NAMESPACE(INTERNAL, _FindReducer(product, syzterm, m_idLeads, m_LS, m_rBaseRing, m_atttributes));
845}
846
847poly SchreyerSyzygyComputation::SchreyerSyzygyNF(poly syz_lead, poly syz_2) const
848{
849  return FROM_NAMESPACE(INTERNAL, _SchreyerSyzygyNF(syz_lead, syz_2, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
850}
851
852poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
853{
854  return FROM_NAMESPACE(INTERNAL, _TraverseTail(multiplier, tail, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
855}
856
857poly SchreyerSyzygyComputation::ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck) const
858{
859  return FROM_NAMESPACE(INTERNAL, _ReduceTerm(multiplier, term4reduction, syztermCheck, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
860}
861
862
863BEGIN_NAMESPACE_NONAME
864   
865static inline int atGetInt(idhdl rootRingHdl, const char* attribute, int def)
866{
867  return ((int)(long)(atGet(rootRingHdl, attribute, INT_CMD, (void*)def)));
868}
869
870END_NAMESPACE   
871
872SchreyerSyzygyComputationFlags::SchreyerSyzygyComputationFlags(idhdl rootRingHdl):
873#ifndef NDEBUG
874    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", TRUE) ),
875#else
876    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", FALSE) ),
877#endif
878    __SYZCHECK__( (BOOLEAN)atGetInt(rootRingHdl, "SYZCHECK", __DEBUG__) ),
879    __LEAD2SYZ__( (BOOLEAN)atGetInt(rootRingHdl, "LEAD2SYZ", 1) ), 
880    __TAILREDSYZ__( (BOOLEAN)atGetInt(rootRingHdl, "TAILREDSYZ", 1) ), 
881    __HYBRIDNF__( (BOOLEAN)atGetInt(rootRingHdl, "HYBRIDNF", 0) ) 
882{   
883  if( __DEBUG__ )
884  {
885    PrintS("SchreyerSyzygyComputationFlags: \n");
886    Print("   DEBUG     : \t%d\n", __DEBUG__);
887    Print("   SYZCHECK  : \t%d\n", __SYZCHECK__);
888    Print("   LEAD2SYZ  : \t%d\n", __LEAD2SYZ__);
889    Print("   TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
890  }
891
892  // TODO: just current setting!
893  assume( rootRingHdl == currRingHdl );
894  assume( rootRingHdl->typ == RING_CMD );
895  assume( rootRingHdl->data.uring == currRing );
896  // move the global ring here inside???
897}
898
899   
900
901
902END_NAMESPACE               END_NAMESPACE_SINGULARXX
903
904
905// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
Note: See TracBrowser for help on using the repository browser.