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

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