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

jengelh-datetimespielwiese
Last change on this file since e98c64 was e98c64, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
adding IGNORETAILS option for computing pure Schreyer Frames (test with 'testAGRFrame')
  • 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//   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  if( m_syzLeads == NULL )
514    ComputeLeadingSyzygyTerms( __LEAD2SYZ__ ); // 2 terms OR 1 term!
515
516  assume( m_syzLeads != NULL );
517
518  ideal& LL = m_syzLeads;
519
520 
521  const int size = IDELEMS(LL);
522
523  TT = idInit(size, 0);
524
525  if( size == 1 && LL->m[0] == NULL )
526    return;
527
528
529  for( int k = size - 1; k >= 0; k-- )
530  {
531    const poly a = LL->m[k]; assume( a != NULL );
532
533    poly a2 = pNext(a);   
534
535    // Splitting 2-terms Leading syzygy module
536    if( a2 != NULL )
537      pNext(a) = NULL;
538
539    if( __IGNORETAILS__ )
540    {
541      TT->m[k] = NULL;
542
543      if( a2 != NULL )
544        p_Delete(&a2, R);
545
546      continue;
547    }
548   
549    TT->m[k] = a2;
550
551    const int r = p_GetComp(a, R) - 1; 
552
553    assume( r >= 0 && r < IDELEMS(T) );
554    assume( r >= 0 && r < IDELEMS(L) );         
555
556    poly aa = leadmonom(a, R); assume( aa != NULL); // :(
557
558    if( ! __HYBRIDNF__ )
559    {
560      poly t = TraverseTail(aa, T->m[r]);
561
562      if( a2 != NULL )
563      {
564        assume( __LEAD2SYZ__ );
565
566        const int r2 = p_GetComp(a2, R) - 1; poly aa2 = leadmonom(a2, R); // :(
567
568        assume( r2 >= 0 && r2 < IDELEMS(T) );
569
570        TT->m[k] = p_Add_q(a2, p_Add_q(t, TraverseTail(aa2, T->m[r2]), R), R);
571
572        p_Delete(&aa2, R);
573      } else
574        TT->m[k] = p_Add_q(t, ReduceTerm(aa, L->m[r], a), R);
575
576    } else
577    {
578      if( a2 == NULL )
579      {
580        aa = p_Mult_mm(aa, L->m[r], R); a2 = m_div.FindReducer(aa, a, m_checker); 
581      }
582      assume( a2 != NULL );
583
584      TT->m[k] = SchreyerSyzygyNF(a, a2); // will copy a2 :(
585
586      p_Delete(&a2, R);
587    }
588
589    p_Delete(&aa, R);   
590  }
591
592  TT->rank = id_RankFreeModule(TT, R);
593}
594
595void SchreyerSyzygyComputation::ComputeLeadingSyzygyTerms(bool bComputeSecondTerms)
596{
597//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
598
599//  const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
600//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
601
602  assume( m_syzLeads == NULL );
603
604  if( bComputeSecondTerms )
605  {
606    assume( __LEAD2SYZ__ );
607//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _Compute2LeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
608    m_syzLeads = Compute2LeadingSyzygyTerms();
609  }
610  else
611  {
612    assume( !__LEAD2SYZ__ );
613   
614    m_syzLeads = Compute1LeadingSyzygyTerms();
615  }
616//    m_syzLeads = FROM_NAMESPACE(INTERNAL, _ComputeLeadingSyzygyTerms(m_idLeads, m_rBaseRing, m_atttributes));
617 
618  // NOTE: set m_LS if tails are to be reduced!
619  assume( m_syzLeads!= NULL );
620
621  if (__TAILREDSYZ__ && (IDELEMS(m_syzLeads) > 0) && !__IGNORETAILS__)
622  {
623    m_LS = m_syzLeads;
624    m_checker.Initialize(m_syzLeads);
625    assume( m_checker.IsNonempty() ); // TODO: this always fails... BUG????
626  }
627}
628
629poly SchreyerSyzygyComputation::SchreyerSyzygyNF(poly syz_lead, poly syz_2) const
630{
631// return FROM_NAMESPACE(INTERNAL, _SchreyerSyzygyNF(syz_lead, syz_2, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
632// poly _SchreyerSyzygyNF(poly syz_lead, poly syz_2,
633//                       ideal L, ideal T, ideal LS,
634//                       const ring r,
635//                       const SchreyerSyzygyComputationFlags attributes)
636// {
637
638  assume( !__IGNORETAILS__ );
639
640  const ideal& L = m_idLeads;
641  const ideal& T = m_idTails;
642  const ring& r = m_rBaseRing;
643
644//   const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
645//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
646//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
647//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
648//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
649//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
650
651  assume( syz_lead != NULL );
652  assume( syz_2 != NULL );
653
654  assume( L != NULL );
655  assume( T != NULL );
656
657  assume( IDELEMS(L) == IDELEMS(T) );
658
659  int  c = p_GetComp(syz_lead, r) - 1;
660
661  assume( c >= 0 && c < IDELEMS(T) );
662
663  poly p = leadmonom(syz_lead, r); // :( 
664  poly spoly = pp_Mult_qq(p, T->m[c], r);
665  p_Delete(&p, r);
666
667
668  c = p_GetComp(syz_2, r) - 1;
669  assume( c >= 0 && c < IDELEMS(T) );
670
671  p = leadmonom(syz_2, r); // :(
672  spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
673  p_Delete(&p, r);
674
675  poly tail = p_Copy(syz_2, r); // TODO: use bucket!?
676
677  while (spoly != NULL)
678  {
679    poly t = m_div.FindReducer(spoly, NULL, m_checker);
680
681    p_LmDelete(&spoly, r);
682
683    if( t != NULL )
684    {
685      p = leadmonom(t, r); // :(
686      c = p_GetComp(t, r) - 1;
687
688      assume( c >= 0 && c < IDELEMS(T) );
689
690      spoly = p_Add_q(spoly, pp_Mult_qq(p, T->m[c], r), r);
691
692      p_Delete(&p, r);
693
694      tail = p_Add_q(tail, t, r);
695    }   
696  }
697
698  return tail;
699}
700
701
702poly SchreyerSyzygyComputation::TraverseTail(poly multiplier, poly tail) const
703{
704  assume( !__IGNORETAILS__ );
705
706  const ideal& L = m_idLeads;
707  const ideal& T = m_idTails;
708//  const ideal& LS = m_LS;
709  const ring& r = m_rBaseRing;
710//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
711
712//  return FROM_NAMESPACE(INTERNAL, _TraverseTail(multiplier, tail, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
713// poly _TraverseTail(poly multiplier, poly tail,
714//                    ideal L, ideal T, ideal LS,
715//                    const ring r,
716//                    const SchreyerSyzygyComputationFlags attributes)
717// {
718
719
720//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
721//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
722//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
723//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
724//   const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
725
726  assume( multiplier != NULL );
727
728  assume( L != NULL );
729  assume( T != NULL );
730
731  poly s = NULL;
732
733  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
734    for(poly p = tail; p != NULL; p = pNext(p))   // iterate over the tail
735      s = p_Add_q(s, ReduceTerm(multiplier, p, NULL), r); 
736
737  return s;
738}
739
740
741
742
743poly SchreyerSyzygyComputation::ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck) const
744{
745  assume( !__IGNORETAILS__ );
746
747  const ideal& L = m_idLeads;
748  const ideal& T = m_idTails;
749//  const ideal& LS = m_LS;
750  const ring& r = m_rBaseRing;
751//  const SchreyerSyzygyComputationFlags& attributes = m_atttributes;
752
753//  return FROM_NAMESPACE(INTERNAL, _ReduceTerm(multiplier, term4reduction, syztermCheck, m_idLeads, m_idTails, m_LS, m_rBaseRing, m_atttributes));
754// poly _ReduceTerm(poly multiplier, poly term4reduction, poly syztermCheck,
755//                 ideal L, ideal T, ideal LS,
756//                 const ring r,
757//                 const SchreyerSyzygyComputationFlags attributes)
758
759
760
761//   const BOOLEAN __DEBUG__      = attributes.__DEBUG__;
762//   const BOOLEAN __SYZCHECK__   = attributes.__SYZCHECK__;
763//   const BOOLEAN __LEAD2SYZ__   = attributes.__LEAD2SYZ__;
764//   const BOOLEAN __HYBRIDNF__   = attributes.__HYBRIDNF__;
765//  const BOOLEAN __TAILREDSYZ__ = attributes.__TAILREDSYZ__;
766
767  assume( multiplier != NULL );
768  assume( term4reduction != NULL );
769
770
771  assume( L != NULL );
772  assume( T != NULL );
773
774  // assume(r == currRing); // ?
775
776  // simple implementation with FindReducer:
777  poly s = NULL;
778
779  if( (!__TAILREDSYZ__) || m_lcm.Check(multiplier) )
780  {
781    // NOTE: only LT(term4reduction) should be used in the following:
782    poly product = pp_Mult_mm(multiplier, term4reduction, r);
783    s = m_div.FindReducer(product, syztermCheck, m_checker);
784    p_Delete(&product, r);
785  }
786
787  if( s == NULL ) // No Reducer?
788    return s;
789
790  poly b = leadmonom(s, r);
791
792  const int c = p_GetComp(s, r) - 1;
793  assume( c >= 0 && c < IDELEMS(T) );
794
795  const poly tail = T->m[c];
796
797  if( tail != NULL )
798    s = p_Add_q(s, TraverseTail(b, tail), r); 
799
800  return s;
801}
802
803
804
805
806
807BEGIN_NAMESPACE_NONAME
808   
809static inline int atGetInt(idhdl rootRingHdl, const char* attribute, long def)
810{
811  return ((int)(long)(atGet(rootRingHdl, attribute, INT_CMD, (void*)def)));
812}
813
814END_NAMESPACE   
815
816SchreyerSyzygyComputationFlags::SchreyerSyzygyComputationFlags(idhdl rootRingHdl):
817#ifndef NDEBUG
818    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", TRUE) ),
819#else
820    __DEBUG__( (BOOLEAN)atGetInt(rootRingHdl,"DEBUG", FALSE) ),
821#endif
822//    __SYZCHECK__( (BOOLEAN)atGetInt(rootRingHdl, "SYZCHECK", __DEBUG__) ),
823    __LEAD2SYZ__( (BOOLEAN)atGetInt(rootRingHdl, "LEAD2SYZ", 1) ), 
824    __TAILREDSYZ__( (BOOLEAN)atGetInt(rootRingHdl, "TAILREDSYZ", 1) ), 
825    __HYBRIDNF__( (BOOLEAN)atGetInt(rootRingHdl, "HYBRIDNF", 0) ),
826    __IGNORETAILS__( (BOOLEAN)atGetInt(rootRingHdl, "IGNORETAILS", 0) ),
827    m_rBaseRing( rootRingHdl->data.uring )
828{   
829  if( __DEBUG__ )
830  {
831    PrintS("SchreyerSyzygyComputationFlags: \n");
832    Print("        DEBUG: \t%d\n", __DEBUG__);
833//    Print("   SYZCHECK  : \t%d\n", __SYZCHECK__);
834    Print("     LEAD2SYZ: \t%d\n", __LEAD2SYZ__);
835    Print("   TAILREDSYZ: \t%d\n", __TAILREDSYZ__);
836    Print("  IGNORETAILS: \t%d\n", __IGNORETAILS__);
837   
838  }
839
840  // TODO: just current setting!
841  assume( rootRingHdl == currRingHdl );
842  assume( rootRingHdl->typ == RING_CMD );
843  assume( m_rBaseRing == currRing );
844  // move the global ring here inside???
845}
846
847   
848
849CReducerFinder::CLeadingTerm::CLeadingTerm(unsigned int _label,  const poly _lt, const ring R):
850    m_sev( p_GetShortExpVector(_lt, R) ),  m_label( _label ),  m_lt( _lt )
851{ }
852
853
854CReducerFinder::~CReducerFinder()
855{
856  for( CReducersHash::const_iterator it = m_hash.begin(); it != m_hash.end(); it++ )
857  {
858    const TReducers& v = it->second;
859    for(TReducers::const_iterator vit = v.begin(); vit != v.end(); vit++ )
860      delete const_cast<CLeadingTerm*>(*vit);
861  }
862}
863
864
865void CReducerFinder::Initialize(const ideal L)
866{
867  assume( m_L == NULL || m_L == L );
868  if( m_L == NULL )
869    m_L = L;
870
871  assume( m_L == L );
872 
873  if( L != NULL )
874  {
875    const ring& R = m_rBaseRing;
876    assume( R != NULL );
877   
878    for( int k = IDELEMS(L) - 1; k >= 0; k-- )
879    {
880      const poly a = L->m[k]; // assume( a != NULL );
881
882      // NOTE: label is k \in 0 ... |L|-1!!!
883      if( a != NULL )
884        m_hash[p_GetComp(a, R)].push_back( new CLeadingTerm(k, a, R) );
885    }
886  }
887}
888
889CReducerFinder::CReducerFinder(const ideal L, const SchreyerSyzygyComputationFlags& flags):
890    SchreyerSyzygyComputationFlags(flags),
891    m_L(const_cast<ideal>(L)), // for debug anyway
892    m_hash()
893{
894  assume( flags.m_rBaseRing == m_rBaseRing );
895  if( L != NULL )
896    Initialize(L);
897}
898
899
900bool CReducerFinder::IsDivisible(const poly product) const
901{
902  const ring& r = m_rBaseRing;
903 
904  const long comp = p_GetComp(product, r);
905  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
906
907  assume( comp >= 0 );
908
909  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
910
911  if( it == m_hash.end() )
912    return false;
913
914  assume( m_L != NULL ); 
915
916  const TReducers& reducers = it->second;
917
918  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
919  {
920    const poly p = (*vit)->m_lt;
921
922    assume( p_GetComp(p, r) == comp );
923
924    const int k = (*vit)->m_label;
925
926    assume( m_L->m[k] == p );
927
928    const unsigned long p_sev = (*vit)->m_sev;
929
930    assume( p_sev == p_GetShortExpVector(p, r) );     
931
932    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
933      continue;
934
935    if( __DEBUG__ )
936    {
937      Print("_FindReducer::Test LS: q is divisible by LS[%d] !:((, diviser is: ", k+1);
938      dPrint(p, r, r, 1);
939    }
940
941    return true;
942  }
943
944  return false;
945}
946
947
948poly CReducerFinder::FindReducer(const poly product, const poly syzterm, const CReducerFinder& syz_checker) const
949{
950  const ring& r = m_rBaseRing;
951
952  assume( product != NULL );
953
954  const ideal& L = m_L; assume( L != NULL ); // for debug/testing only!
955
956  long c = 0;
957
958  if (syzterm != NULL)
959    c = p_GetComp(syzterm, r) - 1;
960
961  assume( c >= 0 && c < IDELEMS(L) );
962
963  if (__DEBUG__ && (syzterm != NULL))
964  {
965    const poly m = L->m[c];
966
967    assume( m != NULL ); assume( pNext(m) == NULL );
968
969    poly lm = p_Mult_mm(leadmonom(syzterm, r), m, r);
970    assume( p_EqualPolys(lm, product, r) );
971
972    //  def @@c = leadcomp(syzterm); int @@r = int(@@c);
973    //  def @@product = leadmonomial(syzterm) * L[@@r];
974
975    p_Delete(&lm, r);   
976  }
977
978  const long comp = p_GetComp(product, r);
979  const unsigned long not_sev = ~p_GetShortExpVector(product, r);
980
981  assume( comp >= 0 );
982
983//   for( int k = IDELEMS(L)-1; k>= 0; k-- )
984//   {
985//     const poly p = L->m[k];
986//
987//     if ( p_GetComp(p, r) != comp )
988//       continue;
989//
990//     const unsigned long p_sev = p_GetShortExpVector(p, r); // to be stored in m_hash!!!
991 
992   // looking for an appropriate diviser p = L[k]...
993  CReducersHash::const_iterator it = m_hash.find(comp); // same module component
994
995  if( it == m_hash.end() )
996    return NULL;
997
998  assume( m_L != NULL );
999
1000  const TReducers& reducers = it->second;
1001
1002  const BOOLEAN to_check = (syz_checker.IsNonempty()); // __TAILREDSYZ__ &&
1003
1004  const poly q = p_New(r); pNext(q) = NULL;
1005
1006  if( __DEBUG__ )
1007    p_SetCoeff0(q, 0, r); // for printing q
1008 
1009  for(TReducers::const_iterator vit = reducers.begin(); vit != reducers.end(); vit++ )
1010  {
1011    const poly p = (*vit)->m_lt;
1012
1013    assume( p_GetComp(p, r) == comp );
1014
1015    const int k = (*vit)->m_label;
1016
1017    assume( L->m[k] == p );
1018
1019    const unsigned long p_sev = (*vit)->m_sev;
1020
1021    assume( p_sev == p_GetShortExpVector(p, r) );     
1022
1023    if( !p_LmShortDivisibleByNoComp(p, p_sev, product, not_sev, r) )
1024      continue;     
1025
1026//     // ... which divides the product, looking for the _1st_ appropriate one!
1027//     if( !p_LmDivisibleByNoComp(p, product, r) ) // included inside  p_LmShortDivisibleBy!
1028//       continue;
1029
1030    p_ExpVectorDiff(q, product, p, r); // (LM(product) / LM(L[k]))
1031    p_SetComp(q, k + 1, r);
1032    p_Setm(q, r);
1033
1034    // cannot allow something like: a*gen(i) - a*gen(i)
1035    if (syzterm != NULL && (k == c))
1036      if (p_ExpVectorEqual(syzterm, q, r))
1037      {
1038        if( __DEBUG__ )
1039        {
1040          Print("_FindReducer::Test SYZTERM: q == syzterm !:((, syzterm is: ");
1041          dPrint(syzterm, r, r, 1);
1042        }   
1043
1044        continue;
1045      }
1046
1047    // while the complement (the fraction) is not reducible by leading syzygies
1048    if( to_check && syz_checker.IsDivisible(q) ) 
1049    {
1050      if( __DEBUG__ )
1051      {
1052        PrintS("_FindReducer::Test LS: q is divisible by LS[?] !:((: ");
1053      }
1054     
1055      continue;
1056    }
1057
1058    p_SetCoeff0(q, n_Neg( n_Div( p_GetCoeff(product, r), p_GetCoeff(p, r), r), r), r);
1059    return q;
1060  }
1061
1062  p_LmFree(q, r);
1063
1064  return NULL;
1065}
1066
1067
1068
1069CLCM::CLCM(const ideal& L, const SchreyerSyzygyComputationFlags& flags):
1070    SchreyerSyzygyComputationFlags(flags), std::vector<bool>(),
1071    m_compute(false), m_N(rVar(flags.m_rBaseRing))
1072{
1073  const ring& R = m_rBaseRing;
1074  assume( flags.m_rBaseRing == R );
1075  assume( R != NULL );
1076
1077  assume( L != NULL );
1078
1079  if( __TAILREDSYZ__ && !__HYBRIDNF__ && (L != NULL))
1080  {
1081    const int l = IDELEMS(L);
1082
1083    assume( l > 0 );
1084
1085    resize(l, false);
1086
1087    for( int k = l - 1; k >= 0; k-- )
1088    {
1089      const poly a = L->m[k]; assume( a != NULL );
1090
1091      for (unsigned int j = m_N; j > 0; j--)
1092        if ( !(*this)[j] )
1093          (*this)[j] = (p_GetExp(a, j, R) > 0);
1094    }
1095
1096    m_compute = true;   
1097  }
1098}
1099
1100
1101bool CLCM::Check(const poly m) const
1102{
1103  assume( m != NULL );
1104  if( m_compute && (m != NULL))
1105  { 
1106    const ring& R = m_rBaseRing;
1107
1108    assume( __TAILREDSYZ__ && !__HYBRIDNF__ );
1109
1110    for (unsigned int j = m_N; j > 0; j--)
1111      if ( (*this)[j] )
1112        if(p_GetExp(m, j, R) > 0)
1113          return true;
1114
1115    return false;
1116
1117  } else return true;
1118}
1119
1120
1121
1122
1123END_NAMESPACE               END_NAMESPACE_SINGULARXX
1124
1125
1126// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
Note: See TracBrowser for help on using the repository browser.