source: git/dyn_modules/syzextra/syzextra.cc @ 5cecde

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